Logo coherent WaveBurst  
Library Reference Guide
Logo
ClusterToFrameNew01_dtmin_approx_cent.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 //
20 // Convert netcluster principal components to NN (Neural Network) frames (NDIMxNDIM pixels)
21 // Author : Gabriele Vedovato
22 // Note : run with ACLiC -> root -l 'ClusterToFrameNew01_dtmin_approx_cent.C+("ifile.root","ofile.root")'
23 //
24 
25 #define XIFO 4
26 
27 #pragma GCC system_header
28 
29 #include "cwb.hh"
30 #include "config.hh"
31 #include "network.hh"
32 #include "wavearray.hh"
33 #include "TString.h"
34 #include "TObjArray.h"
35 #include "TObjString.h"
36 #include "TPaletteAxis.h"
37 #include "TRandom.h"
38 #include "TComplex.h"
39 #include "TMath.h"
40 #include "mdc.hh"
41 #include "watplot.hh"
42 #include <vector>
43 #include <algorithm>
44 #include "TTreeFormula.h"
45 TH2F* hist2D=NULL;
46 TCanvas* canvas;
47 #define REF 0.1
48 #define NPIX 20
49 #define NDIM 8
50 #define NBIN 20
51 #define MAXL 10
52 #define MINL 4
53 #define maxLayers 1025
54 
55 void PlotFrame(std::vector<double>* nn_frame, int cid);
56 
57 double ConvertToFrame(netcluster* pwc, int cid, int nifo, char type, int irate, double &dT,double &dF,double &Fc,int &ninf,std::vector<double>* nn_frame, std::vector<double>* dt);
58 
59 //#define WATPLOT // if defined the event monster plots are produced
60 //#define PLOT_FRAME // frame display
61 //#define distribution
62 
64 #ifdef WATPLOT
65  watplot WTS(const_cast<char*>("wts"));
66 #endif
67 
68  // open input root file
69  TFile* ifile = new TFile(ifName);
70  if(ifile==NULL) {
71  cout << "ClusterToFrame - Error : error opening file : " << ifName.Data() << endl;
72  gSystem->Exit(1);
73  }
74 
75  // get waveburst tree
76  TTree* nn_itree = (TTree *) ifile->Get("waveburst");
77  if(nn_itree==NULL) {
78  cout << "ClusterToFrame - Error : tree waveburst not found !!!" << endl;
79  gSystem->Exit(1);
80  }
81 
82  netcluster* pwc = new netcluster();
83  nn_itree->SetBranchAddress("netcluster",&pwc);
84  int Ntot=nn_itree->GetEntries();
85 
86  // get number of detectors
87  int nIFO = nn_itree->GetUserInfo()->GetSize();
88  if(nIFO==0) {
89  cout<<"ClusterToFrame - Error : wrong user detector list in header tree"<<endl;
90  gSystem->Exit(1);
91  }
92 
93 // open output root file
94  TFile ofile(ofName,"RECREATE");
95 // create output tree
96  TTree* nn_otree = (TTree*)nn_itree->CloneTree(0);
97  //nn_otree->SetName("nn");
98  //nn_otree->SetTitle("nn");
99 
100 // add the NDIM*NDIM pixels vector (nn_frame) for NN analysis and the time vector (dt) which contains the time-duration of each selected pixel
101  std::vector<double>* nn_frame = new vector<double>;
102  std::vector<double>* dt = new vector<double>;
103  nn_frame->resize(NDIM*NDIM);
104  dt->resize(NPIX);
105 
106  int ind=0;
107  double dT=0.;
108  double dF=0.;
109  double Fc=0.;
110  double deltat[NPIX];
111 
112  for (int iu=0; iu<NPIX; iu++) deltat[iu]=0.;
113  int ninf=0;
114 
115 //defining brannces for adding other information
116  nn_otree->Branch("nnframe", &nn_frame, NDIM*NDIM, 0);
117  nn_otree->Branch("dt_corepixels",&dt,NPIX,0);
118  nn_otree->Branch("index",&ind,"ind/I");
119  nn_otree->Branch("t_duration",&dT,"duration/D");
120  nn_otree->Branch("f_duration",&dF,"f_duration/D");
121  nn_otree->Branch("central_f",&Fc,"central_f/D");
122 
123 // read nn_itree
124  int nn_isize = nn_itree->GetEntries();
125  cout << "nn_itree size : " << nn_isize << endl;
126  double SNRf=0.;
127 
128  for(int m=0;m<nn_isize;m++) {
129  if(m%100==0) cout << "ClusterToFrame -> " << m << "/" << nn_isize << endl;
130 
131  nn_itree->GetEntry(m);
132 
133 #ifdef WATPLOT // monster event display
134  WTS.canvas->cd();
135  char fname[1024];
136  sprintf(fname, "monster_event_display_%d.png",m);
137  cout << "write " << fname << endl;
138  WTS.plot(pwc, 1, nIFO, 'L', 0, const_cast<char*>("COLZ")); // use ID=1
139  WTS.canvas->Print(fname);
140  WTS.clear();
141 #endif
142 
143 // using the real conversion-function (ConvertToFrame)
144 
145  ind=m;
146  double SNRfi=0.;
147  int ID = 1; // ID=1 because each netcluster contains only 1 event
148 
149  for (int kk=0; kk<NDIM*NDIM; kk++) (*nn_frame)[kk]=0.;
150  for (int kk=0; kk<NPIX; kk++) (*dt)[kk]=0.;
151  SNRfi=ConvertToFrame(pwc, m, nIFO, 'L', 0, dT, dF, Fc,ninf, nn_frame,dt);
152  SNRf=SNRf+SNRfi; //evaluating the avarage weight of the latest pixel
153 
154 #ifdef PLOT_FRAME // frame display
155  PlotFrame(nn_frame, m);
156 #endif
157 
158 //filling the output tree
159  nn_otree->Fill();
160  }
161 
162 //write and close the output file
163  cout<<"nn_size: "<<nn_isize<<endl;
164  nn_otree->Write();
165  ofile.Close();
166  cout<<ofName<<endl;
167 
168  ifile->Close();
169  cout<<"number events with SNRfi<"<<REF<<": "<<ninf<<endl;
170  cout<<"media importanza ultimo bin: "<<SNRf/nn_isize<<endl;
171  gSystem->Exit(0);
172 }
173 
174 double
175 ConvertToFrame(netcluster* pwc, int n_ev, int nifo, char type, int irate,double &dT, double &dF, double &Fc,int &ninf, std::vector<double>* nn_frame, std::vector<double>* dt) {
176  int cid=1;
177  double FRAME[NDIM][NDIM];
178  if(type != 'L' && type != 'N') return 0.;
179 
180  double RATE = pwc->rate; // original rate>4096
181 
182  std::vector<int>* vint = &(pwc->cList[cid-1]); // pixel list
183 
184  int V = vint->size(); // cluster size
185  if(!V) return 0.;
186 
187  bool optrate=false;
188  if(irate==1) { // extract optimal rate
189  optrate=true;
190  std::vector<int>* pv = pwc->cRate.size() ? &(pwc->cRate[cid-1]) : NULL;
191  irate = pv!=NULL ? (*pv)[0] : 0;
192  }
193 
194  //Layers: number of pixels pixels contained in the frequency range (64-2048) in a particular level of decomposition
195 
196  // int minLayers=1000;
197  // int maxLayers=0;
198  double minTime=1e20;
199  double maxTime=0.;
200  double minFreq=1e20;
201  double maxFreq=0.;
202 
203  //in freq sta prendendo i centri dei pixels mentre nei tempi prendi i limite inf
204  int num_core=0;
205 
206 //loop over the pixels
207  for(int j=0; j<V; j++) { // loop over the pixels
208  netpixel* pix = pwc->getPixel(cid,j);
209  if(!pix->core) continue; // selection of the only core-pixels
210 
211  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
212 
213 // if(pix->layers<minLayers) minLayers=pix->layers;
214 // if(pix->layers>maxLayers) maxLayers=pix->layers;
215 
216 // int minLayers=pow(2,MINL)+1;
217 // int maxLayers=pow(2,MAXL)+1;
218 
219 
220  //deltat[j]=((pix->layers-1)/RATE);
221  int max_ind=0;
222  double rate = RATE/(pix->layers-1);
223  double time=((double)pix->time/rate)/pix->layers;
224  if(time<minTime) minTime=time;
225  if(time>maxTime) {maxTime=time;
226  //max_ind=j;
227  //max_amp=pow(pix->getdata(),2);
228  }
229  num_core=num_core+1;
230  //pix->frequency: nindice in freq
231  //rate/2=df
232  //siccome in freq inizio con pixels mezzi così(pix->frequency*rate/2.) arrivo alla freq centrale-> layers sono quindi sempre dispari
233  //pix->time=itime*(pix->layer)+ifreq {ifreq e itime idici interi di array}
234  //time=itime*dt[s]
235  //freq=ifreq*df[Hz]
236 
237  double freq = pix->frequency*rate/2.;
238  if(freq<minFreq) minFreq=freq;
239  if(freq>maxFreq) maxFreq=freq;
240  }
241  int const ncorepix=num_core;
242  double deltat[ncorepix];
243  for (int u=0;u<ncorepix;u++) deltat[u]=0.;
244  //cout<<"nPIX: "<<NPIX<<" size: "<<dt->size()<<endl;
245  //std::sort(dt->begin(),dt->begin()+NPIX);
246  // for (int u=0;u<NPIX;u++) cout<<"dt: "<<(*dt)[u+n_ev*NPIX]<<endl;
247  // for (int u=0;u<NPIX;u++) cout<<"evento numero "<<n_ev<<" pixel numero: "<<u<<" dt="<<deltat[u]<<endl;
248  //int minRate=RATE/(maxLayers-1);
249  // int maxRate=RATE/(minLayers-1);
250  int minRate=RATE/(pow(2,MAXL));
251  int maxRate=RATE/(pow(2,MINL));
252  dT=maxTime-minTime;
253  // cout<<" dT: "<<dT<<endl;
254  dF=maxFreq-minFreq;
255  // cout<<" dF: "<<dF<<endl;
256  Fc=(maxFreq+minFreq)/2;
257  // cout<<" Fc: "<<Fc<<endl;
258  double mindt = 1./maxRate;
259  double maxdt = 1./minRate;
260  double mindf = minRate/2.;
261  double maxdf = maxRate/2.;
262 
263  //cout << "minRate : " << minRate << "\t\t\t maxRate : " << maxRate << endl;
264  //cout << "minTime : " << minTime << "\t\t\t maxTime : " << maxTime << endl;
265  //cout << "minFreq : " << minFreq << "\t\t\t maxFreq : " << maxFreq << endl;
266  //cout << "mindt : " << mindt << "\t\t\t maxdt : " << maxdt << endl;
267  //cout << "mindf : " << mindf << "\t\t\t maxdf : " << maxdf << endl;
268 
269  int iminTime = int(minTime);
270  int imaxTime = int(maxTime+1);
271  int nTime = (imaxTime-iminTime)*maxRate;
272 
273  if(hist2D) { delete hist2D; hist2D=NULL; }
274  hist2D=new TH2F("WTF", "WTF", nTime, iminTime, imaxTime, 2*maxLayers-2, 0, RATE/2);
275  hist2D->SetXTitle("time, sec");
276  hist2D->SetYTitle("frequency, Hz");
277 
278  double dFreq = (maxFreq-minFreq)/10.>2*maxdf ? (maxFreq-minFreq)/10. : 2*maxdf ;
279  double mFreq = minFreq-dFreq<0 ? 0 : minFreq-dFreq;
280  double MFreq = maxFreq+dFreq>RATE/2 ? RATE/2 : maxFreq+dFreq;
281  hist2D->GetYaxis()->SetRangeUser(mFreq, MFreq);
282 
283  double dTime = (maxTime-minTime)/10.>2*maxdt ? (maxTime-minTime)/10. : 2*maxdt ;
284  double mTime = minTime-dTime<iminTime ? iminTime : minTime-dTime;
285  double MTime = maxTime+dTime>imaxTime ? imaxTime : maxTime+dTime;
286  hist2D->GetXaxis()->SetRangeUser(mTime,MTime);
287  int ncore=0;
288  int npix=0;
289  double Null2=0;
290  double Likelihood2=0;
291  double Null=0;
292  double Likelihood=0;
293  double amp=0.;
294  double ttime=0.;
295  TTree* tS=new TTree("Time_SNR","Time_SNR");
296  tS->Branch("amplitude",&amp,"amplitude/D");
297  tS->Branch("central_time",&ttime,"central_time/D");
298  for(int n=0; n<V; n++) {
299  netpixel* pix = pwc->getPixel(cid,n);
300  if(!pix->core) continue;
301  double rate = RATE/(pix->layers-1);
302  deltat[ncore]=1./rate;
303  ncore=ncore+1;
304  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
305  npix++;
306  double sSNR=0;
307  double wSNR=0;
308  for(int m=0; m<nifo; m++) {
309  sSNR += pow(pix->getdata('S',m),2); // snr whitened reconstructed signal 00
310  sSNR += pow(pix->getdata('P',m),2); // snr whitened reconstructed signal 90
311  wSNR += pow(pix->getdata('W',m),2); // snr whitened at the detector 00
312  wSNR += pow(pix->getdata('U',m),2); // snr whitened at the detector 90
313  }
314  double null = wSNR-sSNR;
315  amp=sSNR;
316  //cout<<"pixel n: "<<n<<" sSNR: "<<sSNR;
317  int iRATE = RATE/(pix->layers-1);
318  int M=maxRate/iRATE;
319  int K=2*(maxLayers-1)/(pix->layers-1);
320  double itime=((double)pix->time/(double)iRATE)/pix->layers;
321  ttime=itime;
322  //cout<<" ttime: "<<ttime<<endl;
323  tS->Fill();
324  int i=(itime-iminTime)*maxRate;
325  int j=pix->frequency*K;
326  if(iRATE!=irate && irate!=0) continue;
327  Null+=null;
328  Likelihood+=sSNR;
329  int L=0;int R=1;while (R < iRATE) {R*=2;L++;}
330  }
331  //cout<<"ncore: "<<ncore<<endl;
332  //cout<<"minTime: "<<minTime<<" maxTime "<<maxTime<<endl;
333  double deltaBIN=(maxTime-minTime)/(NBIN-1);
334  //cout<<"deltaBIN: "<<deltaBIN<<endl;
335  TH1D* h=new TH1D("SNR_distribution","SNR_distribution",NBIN,1,NBIN+1);
336  //double timec[NPIX];
337  double timec2[ncorepix];
338  // for (int u=0; u<NPIX;u++) timec[u]=0.;
339  char cuts[1024];
340  double time_bin[NBIN];
341  for (int nbin=0;nbin<NBIN;nbin++){
342  double mintime_bin=maxTime-(nbin)*deltaBIN;
343  if (nbin==0||nbin==NBIN-1) mintime_bin=mintime_bin-0.001;
344  time_bin[nbin]=mintime_bin;
345  // cout<<"numero bin "<<nbin<<" mintime_bin "<<mintime_bin<<endl;
346  sprintf(cuts,"central_time>%5.5f",mintime_bin);
347  TTreeFormula* treeformula=new TTreeFormula("cuts",cuts,tS);
348  int err = treeformula->Compile(cuts);
349  if(err) {
350  cout << "DisplayROC.C - wrong input cuts " << cuts << endl;
351  //return -1;
352  }
353  double ybin=0.;
354  tS->SetBranchAddress("central_time",&ttime);
355  tS->SetBranchAddress("amplitude",&amp);
356  for(int n=0;n<tS->GetEntries();n++) { // loop over the detected events
357  tS->GetEntry(n); // load event #n
358 // timec[n]=ttime;
359  timec2[n]=ttime;
360  // cout<<"ttime: "<<ttime<<" amp: "<<amp<<" di evento nume: "<<n+1<<" su "<<tS->GetEntries()<<endl;
361  if(treeformula!=NULL) {
362  if(treeformula->EvalInstance()==0) { // skip entry if it does not satisfy selection cuts
363  // cout << "Skip entry number : " << n << endl;
364  continue;
365  }
366  ybin=ybin+amp;
367  // cout<<"n "<<n<<" min time bin: "<<mintime_bin<<" central time: "<<ttime<<" amp: "<<amp<<" xBIN: "<<mintime_bin<<" yBIN: "<<ybin<<endl;
368  }
369  }
370  h->Fill(nbin+1,ybin);
371  ybin=0.;
372  }
373 std::sort(timec2,timec2+ncorepix);
374 //std::sort(timec,timec+NPIX);
375 int index[ncorepix];
376 TMath::Sort(ncorepix,timec2,index);
377  for (int u=0;u<NPIX;u++){
378  if(u<ncorepix) {
379  //dt->push_back(deltat[index[u]]);
380  (*dt)[u] = deltat[index[u]];
381  //cout<<"dt: "<<deltat[index[u]]<<" central_time "<<timec2[u]<<endl;
382  }
383  //else dt->push_back(0.);
384  else (*dt)[u] = 0.;
385 // cout<<""<<endl;
386 // cout<<"deltat "<<deltat[u]<<" dt: "<<(*dt)[u]<<endl;
387  }
388 #ifdef distribution
389 TCanvas* h_canv=new TCanvas("SNR","SNR",0,0,500,500);
390 h->Draw();
391 char h_char[512];
392 sprintf(h_char,"SNRisto/event%i.png",n_ev);
393 h_canv->SaveAs(h_char);
394 #endif
395 double SNRfi0=0.;
396 double SNRfi=0.;
397 int js=1;
398 while (h->GetBinContent(js)==0&&js<NBIN) js=js+1;
399 SNRfi0=(h->GetBinContent(js)/h->GetBinContent(NBIN));
400 int indtemp=0;
401 double ampTOT=0.;
402 double ampFin=0.;
403 double amp2[ncorepix];
404 for (int u=0;u<ncorepix;u++) amp2[u]=0.;
405 double t_limit=0.;
406  for(int iu=0;iu<ncorepix;iu++){
407  tS->SetBranchAddress("central_time",&ttime);
408  tS->SetBranchAddress("amplitude",&amp);
409  char cuts2[1024];
410  sprintf(cuts2,"central_time>%5.5f",timec2[ncorepix-1-iu]-0.001);
411  TTreeFormula* treeformula2=new TTreeFormula("cuts2",cuts2,tS);
412  int err2 = treeformula2->Compile(cuts2);
413  if(err2) {
414  cout << "DisplayROC.C - wrong input cuts " << cuts2 << endl;
415  //return -1;
416  break;
417  }
418  ampTOT=0.;
419  for(int n=0;n<tS->GetEntries();n++) { // loop over the detected events
420  tS->GetEntry(n); // load event #n
421  //cout<<"tempo limite analizzato "<<timec2[ncorepix-1-iu]<<" ttime "<<ttime<<" evento nume: "<<n+1<<" su "<<tS->GetEntries()<<endl;
422  if(treeformula2!=NULL) {
423  if(treeformula2->EvalInstance()==0) { // skip entry if it does not satisfy selection cuts
424  // cout << "Skip entry number : " << n << endl;
425  continue;
426  }
427  ampTOT=ampTOT+amp;
428  //cout<<"timec "<<timec2[ncorepix-1-iu]<< " ampTOT "<<ampTOT<<" amp "<<amp<<endl;
429  }
430  }
431  amp2[iu]=ampTOT;
432  //if(ampTOT>REF&&iu!=0) t_limit= timec2[ncorepix-iu];//così t_limit è l'ultimo da togliere
433  if((ampTOT/h->GetBinContent(NBIN))>REF) {
434  if(iu>0) ampFin=amp2[iu-1];
435  //cout<<"fial amplitude"<<ampFin<<endl;
436  t_limit= timec2[ncorepix-1-iu];//così t_limit è il primo da lasciare
437  //if(iu>0) cout<<"ultimo tempo da buttare:"<<timec2[ncorepix-iu];
438  break;
439  }
440  // if((ampTOT/h->GetBinContent(NBIN))>REF&&iu==0) t_limit=0.;
441 
442  }
443 //cout<<" primo timo da tenere:"<<t_limit<<endl;
444 //SNRfi=ampFin/(h->GetBinContent(NBIN-1));
445 SNRfi=ampFin/(h->GetBinContent(NBIN));
446 //cout<<"SNRfi0: "<<SNRfi0<<" SNRfi: "<<SNRfi<<endl;
447 
448  for(int n=0; n<V; n++) {
449  netpixel* pix = pwc->getPixel(cid,n);
450  if(!pix->core) continue;
451  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
452  double rate = RATE/(pix->layers-1);
453  double sSNR2=0;
454  double wSNR2=0;
455  //cout<<"limit time: "<<t_limit<<" time from direct pixel: "<<((double)pix->time/rate)/pix->layers<<endl;
456  //cout<<"likelihhod: "<<Likelihood2<<endl;
457  if (SNRfi>0.&&SNRfi<=REF&&(((double)pix->time/rate)/pix->layers)>t_limit) continue;
458  for(int m=0; m<nifo; m++) {
459  sSNR2 += pow(pix->getdata('S',m),2); // snr whitened reconstructed signal 00
460  sSNR2 += pow(pix->getdata('P',m),2); // snr whitened reconstructed signal 90
461  wSNR2 += pow(pix->getdata('W',m),2); // snr whitened at the detector 00
462  wSNR2 += pow(pix->getdata('U',m),2); // snr whitened at the detector 90
463  }
464  double null2 = wSNR2-sSNR2;
465  int iRATE = RATE/(pix->layers-1);
466  int M=maxRate/iRATE;
467  int K=2*(maxLayers-1)/(pix->layers-1);
468  double itime=((double)pix->time/(double)iRATE)/pix->layers;
469  ttime=itime;
470  int i=(itime-iminTime)*maxRate;
471  int j=pix->frequency*K;
472  if(iRATE!=irate && irate!=0) continue;
473  Null2+=null2;
474  Likelihood2+=sSNR2;
475  //cout<<"likelihood dopo aggiornamento: "<<Likelihood2<<endl;
476  int L=0;int R=1;while (R < iRATE) {R*=2;L++;}
477  for(int m=0;m<M;m++) {
478  for(int k=0;k<K;k++) {
479  double SNRsum=hist2D->GetBinContent(i+1+m,j+1+k-K/2);
480  if(type=='L') hist2D->SetBinContent(i+1+m,j+1+k-K/2,sSNR2+SNRsum);
481  else hist2D->SetBinContent(i+1+m,j+1+k-K/2,null2+SNRsum);
482  }
483  }
484  }
485 
486  int imin=hist2D->GetNbinsX();
487  int imax=0;
488  int jmin=hist2D->GetNbinsY();
489  int jmax=0;
490 
491  for (int i=0;i<=hist2D->GetNbinsX();i++) {
492  for (int j=0;j<=hist2D->GetNbinsY();j++) {
493  double X = hist2D->GetXaxis()->GetBinCenter(i)+hist2D->GetXaxis()->GetBinWidth(i)/2.;
494  double Y = hist2D->GetYaxis()->GetBinCenter(j)+hist2D->GetYaxis()->GetBinWidth(j)/2.;
495  double Z = hist2D->GetBinContent(i,j);
496  if(Z) {
497  if(i<imin) imin=i;
498  if(i>imax) imax=i;
499  if(j<jmin) jmin=j;
500  if(j>jmax) jmax=j;
501  }
502  //if(Z) cout << i << " " << j << " " << X << " " << Y << " " << Z << endl;
503  }
504  }
505 
506  //cout << "imin " << imin << " imax " << imax << endl;
507  //cout << "jmin " << jmin << " jmax " << jmax << endl;
508 
509 
510  int di = imax-imin+1;
511  int dj = jmax-jmin+1;
512 
513  //cout << "di " << di << " dj " << dj << endl;
514 
515  int ri = NDIM-di%NDIM;
516  int rj = NDIM-dj%NDIM;
517  if (di%NDIM==0) ri=0;
518  if (dj%NDIM==0) rj=0;
519 /* int lri = ri/2;
520  int rri = ri-lri;
521  //cout << "lri " << lri << " rri " << rri << endl;
522 
523  imin-=lri;
524  imax+=rri;
525 
526  int lrj = rj/2;
527  int rrj = rj-lrj;
528  //cout << "lrj " << lrj << " rrj " << rrj << endl;
529 
530  jmin-=lrj;
531  jmax+=rrj;
532 */
533 //cout << "di: " << di << " ri: " << ri << endl;
534 // cout << "dj: " << dj << " rj: " << rj << endl;
535 // cout << "imin " << imin << " imax " << imax << endl;
536 // cout << "jmin " << jmin << " jmax " << jmax << endl;
537  bool ripi=false;
538  bool ripj=false;
539 //cout<<"resto diverso da zero"<<endl;
540 
541  int lri = ri/2;
542  int rri = ri-lri;
543  int lrj = rj/2;
544  int rrj = rj-lrj;
545 
546 
547  if (ri<4 || di<=4){
548  di = di+ri;
549  imin-=lri;
550  imax+=rri;
551  //imin-=ri;
552  }
553 
554  else {
555  imin+=di%NDIM;
556  di = di-di%NDIM;
557  ripi=true;
558  }
559 
560  if (rj<4 || dj<=4){
561  dj = dj+rj;
562  jmin-=lrj;
563  jmax+=rrj;
564 
565  //jmax+=rj;
566  }
567  else {
568  jmax-=dj%NDIM;
569  dj = dj-dj%NDIM;
570  ripj=true;
571  }
572 
573 // cout << "imin " << imin << " imax " << imax << endl;
574 // cout << "jmin " << jmin << " jmax " << jmax << endl;
575 //cout<<" def limiti"<<endl;
576  int xdim = (imax-imin+1)/NDIM;
577  int ydim = (jmax-jmin+1)/NDIM;
578  //cout << "xdim " << xdim << " ydim " << ydim << endl;
579 //cout<<" def xdim,ydim"<<endl;
580 
581  for(int i=0;i<NDIM;i++) for(int j=0;j<NDIM;j++) FRAME[i][j]=0;
582 
583  for (int i=0;i<=hist2D->GetNbinsX();i++) {
584 //cout<<"i, j "<<i<<" "<<" i min "<<imin<<xmax<<" xdim "<<xdim<<" jmin "<<jmin<<" ydim "<<ydim<<endl;
585  for (int j=0;j<=hist2D->GetNbinsY();j++) {
586 
587  bool ripi2=false;
588  bool ripj2=false;
589 
590  double X = hist2D->GetXaxis()->GetBinCenter(i)+hist2D->GetXaxis()->GetBinWidth(i)/2.;
591  double Y = hist2D->GetYaxis()->GetBinCenter(j)+hist2D->GetYaxis()->GetBinWidth(j)/2.;
592  double Z = hist2D->GetBinContent(i,j);
593  if(i>imax) continue;
594  if(j<jmin) continue;
595 
596  // if((i<imin)&&(ri<4|| di<=4)) continue;
597  // if((j>jmax)&&(rj<4|| dj<=4)) continue;
598 
599  // if((i<imin)&&!(ri<4|| di<=4)) i=i+xdim;
600  // if((j>jmax)&&!(rj<4|| dj<=4)) j=j-ydim;
601 
602  if (i<imin && ripi==false) continue;
603  if (j<jmin && ripj==false) continue;
604 
605  if (i<imin && ripi==true) {i=i+xdim;ripi2=true;}
606  if (j<jmin && ripj==true) {j=j-ydim;ripj2=true;}
607 
608 //cout<<"i, j "<<i<<" "<<j<<" i min "<<imin<<" xdim "<<xdim<<" jmin "<<jmin<<" ydim "<<ydim<<endl;
609  int I = (i-imin)/xdim;
610  int J = (j-jmin)/ydim;
611 //cout << I << " " << J << endl;
612  FRAME[I][J]+=Z;
613 //cout << "sono ui" << endl;
614 
615  //if((i-xdim<imin)&&!(ri<4|| di<=4)) i=i-xdim;
616  //if((j+ydim>jmax)&&!(rj<4|| dj<=4)) j=j+ydim;
617 
618  if (ripi2==true) i=i-xdim;
619  if (ripj2==true) j=j+ydim;
620  //if(Z) cout << i << " " << j << " " << X << " " << Y << " " << Z << endl;
621  //if(Z) cout << I << " " << J << " " << FRAME[I][J] << endl;
622  }
623  }
624 
625  // fill nn frame
626  //for(int i=0;i<NDIM;i++) for(int j=0;j<NDIM;j++) nn_frame->push_back(FRAME[i][j]);
627  for(int i=0;i<NDIM;i++) for(int j=0;j<NDIM;j++) (*nn_frame)[i*NDIM+j] = FRAME[i][j];
628  // normalization
629  double norm=0;
630  int nn_size = nn_frame->size();
631  for(int i=0;i<nn_size;i++) norm+=(*nn_frame)[i];
632  for(int i=0;i<nn_size;i++) (*nn_frame)[i]/=norm;
633  if (SNRfi<=REF) ninf=ninf+1;
634  delete h;
635  return SNRfi0;
636 
637 //cout<<" fine funzione" <<endl;
638 }
639 
640 void
641 PlotFrame(std::vector<double>* nn_frame, int cid) {
642 
643  // get dimension of the frame
644  int nframe = nn_frame->size();
645  if(nframe != pow(int(sqrt(nframe)),2)) {
646  cout << "PlotFrame - Error : size is not a square number " << endl;
647  gSystem->Exit(1);
648  }
649  nframe = sqrt(nframe);
650 
651  if(hist2D) { delete hist2D; hist2D=NULL; }
652  hist2D=new TH2F("frame", "frame", nframe, 0, nframe, nframe, 0, nframe);
653 
654  hist2D->SetStats(kFALSE);
655  hist2D->SetTitleFont(12);
656  hist2D->SetFillColor(kWhite);
657 
658  hist2D->GetXaxis()->SetNdivisions(506);
659  hist2D->GetXaxis()->SetLabelFont(42);
660  hist2D->GetXaxis()->SetLabelOffset(0.014);
661  hist2D->GetXaxis()->SetTitleOffset(1.4);
662  hist2D->GetYaxis()->SetTitleOffset(1.2);
663  hist2D->GetYaxis()->SetNdivisions(506);
664  hist2D->GetYaxis()->SetLabelFont(42);
665  hist2D->GetYaxis()->SetLabelOffset(0.01);
666  hist2D->GetZaxis()->SetLabelFont(42);
667  hist2D->GetZaxis()->SetNoExponent(false);
668  hist2D->GetZaxis()->SetNdivisions(506);
669 
670  hist2D->GetXaxis()->SetTitleFont(42);
671  hist2D->GetXaxis()->SetTitle("Time");
672  hist2D->GetXaxis()->CenterTitle(true);
673  hist2D->GetYaxis()->SetTitleFont(42);
674  hist2D->GetYaxis()->SetTitle("Frequency");
675  hist2D->GetYaxis()->CenterTitle(true);
676 
677  hist2D->GetZaxis()->SetTitleOffset(0.6);
678  hist2D->GetZaxis()->SetTitleFont(42);
679  hist2D->GetZaxis()->CenterTitle(true);
680 
681  hist2D->GetXaxis()->SetLabelSize(0.03);
682  hist2D->GetYaxis()->SetLabelSize(0.03);
683  hist2D->GetZaxis()->SetLabelSize(0.03);
684 
685  for(int i=0;i<nframe;i++) {
686  for(int j=0;j<nframe;j++) {
687  //cout << i << " " << j << " " << (*nn_frame)[i*nframe+j] << endl;
688  hist2D->SetBinContent(i+1,j+1,(*nn_frame)[i*nframe+j]);
689  }
690  }
691 
692 
693  canvas= new TCanvas("mra", "mra", 200, 20, 600, 600);
694  canvas->Clear();
695  canvas->ToggleEventStatus();
696  canvas->SetGridx();
697  canvas->SetGridy();
698  canvas->SetFillColor(kWhite);
699  canvas->SetRightMargin(0.10);
700  canvas->SetLeftMargin(0.10);
701  canvas->SetBottomMargin(0.13);
702  canvas->SetBorderMode(0);
703 
704  // remove the red box around canvas
705  gStyle->SetFrameBorderMode(0);
706  gROOT->ForceStyle();
707 
708  gStyle->SetTitleH(0.050);
709  gStyle->SetTitleW(0.95);
710  gStyle->SetTitleY(0.98);
711  gStyle->SetTitleFont(12,"D");
712  gStyle->SetTitleColor(kBlue,"D");
713  gStyle->SetTextFont(12);
714  gStyle->SetTitleFillColor(kWhite);
715  gStyle->SetLineColor(kWhite);
716  gStyle->SetNumberContours(256);
717  gStyle->SetCanvasColor(kWhite);
718  gStyle->SetStatBorderSize(1);
719 
720  hist2D->Draw("COLZ");
721 
722  // change palette's width
723  canvas->Update();
724  TPaletteAxis *palette = (TPaletteAxis*)hist2D->GetListOfFunctions()->FindObject("palette");
725  palette->SetX1NDC(0.91);
726  palette->SetX2NDC(0.933);
727  palette->SetTitleOffset(0.92);
728  palette->GetAxis()->SetTickSize(0.01);
729  canvas->Modified();
730 
731  char fname[1024];
732  sprintf(fname, "nn_frame_display_%d.png",cid);
733  cout << "write " << fname << endl;
734  canvas->Print(fname);
735 }
736 
char ofile[1024]
std::vector< vector_int > cRate
Definition: netcluster.hh:398
TString ofName
double M
Definition: DrawEBHH.C:13
int irate
size_t frequency
Definition: netpixel.hh:111
int n
Definition: cwb_net.C:28
TString("c")
int palette
Definition: DrawGnetwork2.C:17
int ID
Definition: TestMDC.C:70
netpixel pix(nifo)
netcluster * pwc
Definition: cwb_job_obj.C:38
void ClusterToFrameNew01_dtmin_approx_cent(TString ifName, TString ofName)
size_t layers
Definition: netpixel.hh:112
int m
Definition: cwb_net.C:28
std::vector< vector_int > cList
Definition: netcluster.hh:397
int j
Definition: cwb_net.C:28
i drho i
double rate
Definition: netcluster.hh:378
int ncore
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:150
bool core
Definition: netpixel.hh:120
void clear()
Definition: watplot.hh:58
#define nIFO
TCanvas * canvas
Definition: watplot.hh:192
wavearray< double > freq
Definition: Regression_H1.C:79
wavearray< double > h
Definition: Regression_H1.C:25
watplot * WTS
Definition: ChirpMass.C:119
i() int(T_cor *100))
char fname[1024]
int k
void PlotFrame(std::vector< double > *nn_frame, int cid)
size_t time
Definition: netpixel.hh:110
TFile * ifile
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:413
int npix
double dt
#define RATE
int nifo
wavearray< int > index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double * itime
double dFreq
Definition: DrawPSD.C:17
double dT
Definition: testWDM_5.C:12
float rate
Definition: netpixel.hh:113
double getdata(char type='R', size_t n=0)
Definition: netpixel.hh:74
float wSNR[PE_MAX_EVENTS][NIFO_MAX]
Definition: cwb_pereport.C:336
double ConvertToFrame(netcluster *pwc, int cid, int nifo, char type, int irate, double &dT, double &dF, double &Fc, int &ninf, std::vector< double > *nn_frame, std::vector< double > *dt)