Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ChirpMass.C
Go to the documentation of this file.
1 // this macro shows how to read the eBBH results from output ROOT file produced with CWB_Plugin_eBBH.C
2 // Author : G.Vedovato
3 
4 #define NCHIRP_MAX 4
5 
6 // ------------------------------------------------------------
7 // ebbh structure
8 // ------------------------------------------------------------
9 
10 struct EBBH { // structure for output for estimated parameters
11 
12  TString name;
13 
14  std::vector<TString> ifo;
15 
16  float isnr;
17  float osnr;
18 
19  double seggps; // segment start time
20 
21  float parms[NIFO_MAX]; // null pixel statistic, for each detector
22 
23  float mchirp;
24  float mass[2];
25  float spin[6]; // spin
26  float e0;
27  float rp0;
28  float dist;
29  float redshift;
30 
31  float ch_mass[NCHIRP_MAX]; // chirp mass
32  float ch_tmerger[NCHIRP_MAX]; // merger time
33  float ch_merr[NCHIRP_MAX]; // mchirp error
34  float ch_mchi2[NCHIRP_MAX]; // mchirp chi2
35  float ch_energy[NCHIRP_MAX]; // chirp energy
36 
37  float ei; // eccentricity index
38 
39  netcluster nc; // netcluster
40 
41  wavearray<double> wdat[NIFO_MAX]; // whitened data (in the vREC time range)
42  wavearray<double> winj[NIFO_MAX]; // whiten injected waveforms
43  wavearray<double> sinj[NIFO_MAX]; // strain injected waveforms
44  wavearray<double> wrec[NIFO_MAX]; // whiten reconstructed waveforms
45  wavearray<double> srec[NIFO_MAX]; // strain injected waveforms
46 
47  wavearray<double> wsim[2]; // strain simulated waveforms hp,hc
48 };
49 
50 // ------------------------------------------------------------
51 // functions
52 // ------------------------------------------------------------
53 
54 int ReadDataFromROOT(TString ifile, EBBH* ebbh);,
55 
56 void PlotMultiChirp(EBBH* ebbh);
57 void PlotMultiChirp2(EBBH* ebbh, int ifo=0, TString ptype="", TString wtype="");
58 void PlotMChirp(EBBH* ebbh);
59 void PlotChirp(EBBH* ebbh, int mch_order, double& tmerger, double& xmin, double& xmax, EColor color);
60 
61 double GetTimeBoundaries(wavearray<double> x, double P, double& bT, double& eT);
62 wavearray<double> GetAlignedWaveform(wavearray<double>* wf, wavearray<double>* wref);
63 wavearray<double> GetDifWaveform(wavearray<double>* wf1, wavearray<double>* wf2);
64 
65 double GetFitParameter(double mchirp);
66 void PlotCluster(watplot* WTS, netcluster* pwc, double fLow, double fHigh);
67 
68 void FindChirp(EBBH* ebbh, int mch_order, double& tmerger, double& xmin, double& xmax, int ifo, TString ptype, TString wtype);
69 
70 void FillHistogram(EBBH* ebbh, TH2F&* hN, TH2F&* hL, TH2F&* heT, TH2F&* heF, double zmax_thr=0.);
71 void PlotHistogram(EBBH* ebbh, TH2F&* hH);
72 double ChirpIntegralHistogram(EBBH* ebbh, TH2F&* h2d, double dMc);
73 
74 void DrawMChirpInstantaneousFrequency(EBBH* ebbh, TString wtype, int ifo);
75 
76 void GetCBC(EBBH* ebbh);
77 void GetEBBH(EBBH* ebbh);
78 
79 void Init();
80 void Clear();
81 
82 // -----------------------------------------------------------------------------------------------------------
83 // defines
84 // -----------------------------------------------------------------------------------------------------------
85 
86 #define FLOW 16
87 #define FHIGH 1024
88 
89 // default threshold (used in network::likelihoodWP)
90 /*
91 #define CHI2_THR 2.0
92 #define TMERGER_CUT 1.e20
93 #define ZMAX_THR 0.0
94 */
95 
96 
97 //#define CHI2_THR -2.5 // if negative the chirp pixels are tagged with pix->likelihood=0 after the selection
98 //#define TMERGER_CUT 0.1
99 #define CHI2_THR -4.5 // if negative the chirp pixels are tagged with pix->likelihood=0 after the selection
100 #define TMERGER_CUT 0.01
101 #define ZMAX_THR 0.2
102 
103 
104 
105 // ------------------------------------------------------------
106 // global variables
107 // ------------------------------------------------------------
108 
109 int gTYPE;
110 int gENTRY;
111 bool gPLOT;
112 TString gOFILE;
113 
114 CWB::mdc* MDC;
115 watplot* WTS;
116 watplot* PCH;
117 CWB::STFT* stft;
118 TF1* fchirp;
119 
120 TH2F* hN;
121 TH2F* hL;
122 TH2F* heT;
123 TH2F* heF;
124 
125 TF1* xfit[NCHIRP_MAX];
126 TGraphErrors* xgr[NCHIRP_MAX];
127 
128 TCanvas* canvas;
129 
130 int ChirpMass(TString ifroot, int gtype=0, int entry=0, int ifo=0, EBBH* xebbh=NULL) {
131  //
132  // ifroot : input root file
133  //
134  // gtype : plot type
135  // if gtype<0 the plot is saved with name with the following format : gTYPE_gtype_file_name.png
136  //
137  // gtype = 0 -> no plots, only computation of ei
138  // gtype = 1 -> plot chirp fit function with rec spectrogram
139  // gtype = 2 -> plot chirp fit function with likelihood
140  // gtype = 3 -> plot chirp fit function with likelihood without the main chirp
141  // show residuals after chirp1 removal ( only with PlotMultiChirp2 & "like" )
142  // gtype = 4 -> plot with inj spectrogram
143  // gtype = 5 -> plot likelihood TF after with pixels selected with zmax_thr
144  // gtype = 6 -> plot chirp mass fit in f^-3/8 vs time plane
145  // gtype = 7 -> plot chirp mass fit in f^-3/8 vs time plane for all chirp
146  // gtype = 8 -> plot ifo waveform in time
147  // gtype = 9 -> plot ifo waveform in frequency
148  // gtype = 10 -> plot whitened data spectrogram
149  // gtype = 11 -> plot inj waveform spectrogram
150  // gtype = 12 -> plot rec waveform spectrogram
151  // gtype = 13 -> plot inj waveform instantaneous frequency mchirp
152  // gtype = 14 -> plot rec waveform instantaneous frequency mchirp
153  // gtype = 15 -> plot sim waveform (GetCBC) instantaneous frequency mchirp
154  // gtype = 16 -> plot sim waveform (GetEBBH) instantaneous frequency mchirp
155  //
156  // entry : entry id of the event in the root file
157  //
158  // ifo : detector id : 0,1,... (see definition in user_parameters.C)
159  //
160  // xebbh : pointer to the EBBH structure : used by a macro to return the ebbh parameters
161  //
162  // Ex : root 'ChirpMass.C("root_file.root",1,0,0)'
163  //
164 
165  Init();
166 
167  if(gtype<-16 || gtype>16) {cout << "Error, gtype parameter value not allowed : [-1:11]" << endl;exit(1);}
168  if(ifroot=="") {cout << "Error, input root file name not defined" << endl;exit(1);}
169 
170  gPLOT = gtype<0 ? true : false;
171  gTYPE = abs(gtype);
172 
173  if(gPLOT) gROOT->SetBatch(true);
174 
175  gENTRY = entry;
176 
177  TString fDir = gSystem->DirName(ifroot);
178  TString fName = gSystem->BaseName(ifroot);
179 
180  if(fName.Contains("*")) {
181  TString fTag = fName;
182  fTag.ReplaceAll("*","");
183  vector<TString> fList = CWB::Toolbox::getFileListFromDir(fDir, ".root", "wave_", fTag, true);
184  if(fList.size()==0) {cout << "Error, no file selected" << endl;exit(1);}
185  if(fList.size()>1) {cout << "Error, multiple files selected" << endl;exit(1);}
186 
187  ifroot = fList[0];
188  fDir = gSystem->DirName(ifroot);
189  fName = gSystem->BaseName(ifroot);
190  }
191 
192  CWB::Toolbox::checkFile(ifroot);
193 
194  gOFILE = fName;
195  gOFILE.ReplaceAll(".root",".png"); // output file name for plot
196  gOFILE = TString::Format("gTYPE_%d_%s",gTYPE,gOFILE.Data());
197  if(gPLOT) cout << endl << "Output Plot File Name : " << gOFILE << endl << endl;
198 
199  CWB::Toolbox::checkFile(ifroot);
200 
201  EBBH ebbh;
202  MDC = new CWB::mdc();
203 
204  // read ebbh from input root file
205  int err = ReadDataFromROOT(ifroot, &ebbh);
206  if(err) return -100000;
207 
208  if(gTYPE==15) GetCBC(&ebbh); // get simulated CBC waveforms
209  if(gTYPE==16) GetEBBH(&ebbh); // get simulated EBBH waveforms
210 
211  // set waveforms start time at start segment time
212  ebbh.wdat[ifo].start(ebbh.seggps ? ebbh.wdat[ifo].start()-ebbh.seggps : 0);
213  ebbh.winj[ifo].start(ebbh.seggps ? ebbh.winj[ifo].start()-ebbh.seggps : 0);
214  ebbh.wrec[ifo].start(ebbh.seggps ? ebbh.wrec[ifo].start()-ebbh.seggps : 0);
215  ebbh.srec[ifo].start(ebbh.seggps ? ebbh.srec[ifo].start()-ebbh.seggps : 0);
216 
217  if(gTYPE==0) {PlotMultiChirp2(&ebbh); if(xebbh) {*xebbh=ebbh; Clear(); return 0;} else exit(0);}
218  if(gTYPE==1) {PlotMultiChirp2(&ebbh, ifo, "stft", "wrec");return 0;}
219  if(gTYPE==2 || gTYPE==3) {PlotMultiChirp2(&ebbh, ifo, "like", "wrec");return 0;}
220  if(gTYPE==4) {PlotMultiChirp2(&ebbh, ifo, "stft", "winj");return 0;}
221  if(gTYPE==5) {
222  FillHistogram(&ebbh, hN, hL, heT, heF, ZMAX_THR);
223  PlotHistogram(&ebbh, hL); return 0;
224 /*
225  FillHistogram(&ebbh, hN, hL, heT, heF, 0);
226  double cenergy = ChirpIntegralHistogram(&ebbh, hL, 0);
227  PlotHistogram(&ebbh, hL); return 0;
228  TH1F* hC = new TH1F("hC","hC",41,-20,20);
229  for(int i=0;i<=40;i++) {
230  double dMc = i-20;
231  double cenergy = ChirpIntegralHistogram(&ebbh, hL, dMc);
232  cout << "CHIRP ENERGY : " << dMc << " -> " << cenergy << endl;
233  hC->SetBinContent(i,cenergy);
234  }
235  hC->Draw();
236  return 0;
237 */
238  //PlotHistogram(&ebbh, hL); return 0;
239  //PlotHistogram(&ebbh, heF); return 0;
240  //PlotHistogram(&ebbh, heT); return 0;
241  //PlotHistogram(&ebbh, hN); return 0;
242  }
243  if(gTYPE==6) {PlotMChirp(&ebbh); return 0;}
244  if(gTYPE==7) {PlotMultiChirp(&ebbh);return 0;}
245  if(gTYPE==8) {
246  MDC->SetZoom(0.999);
247  watplot* plot=NULL;
248  if(ebbh.winj[ifo].size()!=0) {
249  plot=MDC->Draw(ebbh.winj[ifo],MDC_TIME,"ALP ZOOM",kRed);
250  MDC->Draw(ebbh.wrec[ifo],MDC_TIME,"SAME",kBlack);
251  } else {
252  plot=MDC->Draw(ebbh.wrec[ifo],MDC_TIME,"ALP ZOOM",kRed);
253  }
254  if(plot) plot->graph[0]->SetTitle(TString::Format("%s : Rec(black), Inj(red))",ebbh.name.Data()));
255  if(gPLOT) {stft->Print(gOFILE);exit(0);} else return 0;
256  }
257  if(gTYPE==9) {
258  MDC->SetZoom(0.999);
259  watplot* plot=NULL;
260  if(ebbh.winj[ifo].size()!=0) {
261  plot=MDC->Draw(ebbh.winj[ifo],MDC_FFT,"ALP ZOOM",kRed);
262  MDC->Draw(ebbh.wrec[ifo],MDC_FFT,"SAME",kBlack);
263  } else {
264  plot=MDC->Draw(ebbh.wrec[ifo],MDC_FFT,"ALP ZOOM",kRed);
265  }
266  if(plot) plot->graph[0]->SetTitle(TString::Format("%s : Rec(black), Inj(red))",ebbh.name.Data()));
267  if(gPLOT) {MDC->GetWatPlot()->canvas->Print(gOFILE);exit(0);} else return 0;
268  }
269  if(gTYPE==10) {MDC->Draw(ebbh.wdat[ifo],MDC_TF);if(gPLOT) {MDC->GetSTFT()->Print(gOFILE);exit(0);} else return 0;}
270  if(gTYPE==11) {MDC->Draw(ebbh.winj[ifo],MDC_TF);if(gPLOT) {MDC->GetSTFT()->Print(gOFILE);exit(0);} else return 0;}
271  if(gTYPE==12) {MDC->Draw(ebbh.wrec[ifo],MDC_TF);if(gPLOT) {MDC->GetSTFT()->Print(gOFILE);exit(0);} else return 0;}
272 
273  if(gTYPE==13) {DrawMChirpInstantaneousFrequency(&ebbh, "winj", ifo); return;}
274  if(gTYPE==14) {DrawMChirpInstantaneousFrequency(&ebbh, "wrec", ifo); return;}
275  if(gTYPE==15) {DrawMChirpInstantaneousFrequency(&ebbh, "wsim", ifo); return;}
276 
277  return 0;
278 }
279 
280 void Init() {
281 
282  MDC = NULL;
283  WTS = NULL;
284  PCH = NULL;
285  stft = NULL;
286 
287  hN = NULL;
288  hL = NULL;
289  heT = NULL;
290  heF = NULL;
291 
292  for(int i=0;i<NCHIRP_MAX;i++) {
293  xfit[i] = NULL;
294  xgr[i] = NULL;
295  }
296 
297  fchirp = NULL;
298 
299  canvas = NULL;
300 }
301 
302 void Clear() {
303 
304  if(MDC != NULL) delete MDC;
305  if(WTS != NULL) delete WTS;
306  if(PCH != NULL) delete PCH;
307  if(stft != NULL) delete stft;
308 
309  if(hN != NULL) delete hN;
310  if(hL != NULL) delete hL;
311  if(heT != NULL) delete heT;
312  if(heF != NULL) delete heF;
313 
314  for(int i=0;i<NCHIRP_MAX;i++) {
315  if(xfit[i] != NULL) delete xfit[i];
316  if(xgr[i] != NULL) delete xgr[i];
317  }
318 
319  if(fchirp != NULL) delete fchirp;
320 
321  if(canvas != NULL) delete canvas;
322 }
323 
324 int ReadDataFromROOT(TString ifile, EBBH* ebbh) {
325 //
326 // ----------------------------------------------------
327 // ROOT Output PE Parameters
328 // ----------------------------------------------------
329 
330  TFile* froot = new TFile(ifile,"READ");
331  TTree* itree = froot->Get("waveburst");
332  if(itree==NULL) {cout << "ReadDataFromROOT - Error : no waveburst present in the file" << endl;return 1;}
333 
334  // get detector list
335  TList* list = itree->GetUserInfo();
336  int nIFO=list->GetSize();
337  if (nIFO==0) {cout << "ReadDataFromROOT - Error : no ifo present in the tree" << endl;return 1;}
338  for (int n=0;n<list->GetSize();n++) {
339  detector* pDetector;
340  pDetector = (detector*)list->At(n);
341  ebbh->ifo.push_back(pDetector->Name);
342  detectorParams dParams = pDetector->getDetectorParams();
343  //pDetector->print();
344  }
345 
346  itree->SetBranchAddress("ndim",&nIFO);
347  itree->GetEntry(gENTRY);
348 
349  double seggps=0; // segment start time
350  float crate; // netcluster::rate
351  float parms[NIFO_MAX]; // ebbh parameters
352  float eBBH[4]; // ebbh parameters
353  float range[2]; // distance (Mpc)
354  float iSNR[NIFO_MAX];
355  float oSNR[NIFO_MAX];
356  std::vector<netpixel>* cluster;
357  wavearray<double>* wDAT[NIFO_MAX];
358  wavearray<double>* wINJ[NIFO_MAX];
359  wavearray<double>* wREC[NIFO_MAX];
360  wavearray<double>* sREC[NIFO_MAX];
361 
362  for(int n=0;n<nIFO;n++) {
363  wDAT[n] = new wavearray<double>;
364  wINJ[n] = new wavearray<double>;
365  wREC[n] = new wavearray<double>;
366  sREC[n] = new wavearray<double>;
367  }
368 
369  cluster = new std::vector<netpixel>;
370 
371  itree->SetBranchAddress("mass",ebbh->mass);
372  itree->SetBranchAddress("eBBH",eBBH);
373  itree->SetBranchAddress("range",range);
374  itree->SetBranchAddress("spin",ebbh->spin);
375  itree->SetBranchAddress("iSNR",iSNR);
376  itree->SetBranchAddress("oSNR",oSNR);
377 
378  itree->SetBranchAddress("ebbh_parms",parms);
379  for(int n=0;n<nIFO;n++) {
380  itree->SetBranchAddress(TString::Format("ebbh_wDAT_%d",n).Data(),&wDAT[n]);
381  itree->SetBranchAddress(TString::Format("ebbh_wINJ_%d",n).Data(),&wINJ[n]);
382  itree->SetBranchAddress(TString::Format("ebbh_wREC_%d",n).Data(),&wREC[n]);
383  itree->SetBranchAddress(TString::Format("ebbh_sREC_%d",n).Data(),&sREC[n]);
384  }
385 
386  itree->SetBranchAddress("ebbh_ch_mass", ebbh->ch_mass);
387  itree->SetBranchAddress("ebbh_ch_tmerger",ebbh->ch_tmerger);
388  itree->SetBranchAddress("ebbh_ch_merr", ebbh->ch_merr);
389  itree->SetBranchAddress("ebbh_ch_mchi2", ebbh->ch_mchi2);
390  itree->SetBranchAddress("ebbh_ch_energy", ebbh->ch_energy);
391 
392  itree->SetBranchAddress("ebbh_seggps",&seggps);
393  itree->SetBranchAddress("ebbh_crate",&crate);
394  itree->SetBranchAddress("ebbh_cluster",&cluster);
395 
396  itree->GetEntry(gENTRY);
397 
398  ebbh->e0 = eBBH[1];
399  ebbh->rp0 = eBBH[2];
400  ebbh->redshift = eBBH[3];
401  ebbh->dist = range[1];
402  ebbh->isnr = 0;
403  for(int n=0;n<nIFO;n++) ebbh->isnr += iSNR[n];
404  ebbh->isnr = sqrt(ebbh->isnr);
405  ebbh->osnr = 0;
406  for(int n=0;n<nIFO;n++) ebbh->osnr += oSNR[n];
407  ebbh->osnr = sqrt(ebbh->osnr);
408 
409  for(int n=0;n<nIFO;n++) {
410  ebbh->wdat[n] = GetAlignedWaveform(wDAT[n],wREC[n]);
411  ebbh->winj[n] = wINJ[n]->size() ? GetAlignedWaveform(wINJ[n],wREC[n]) : 0;
412  ebbh->wrec[n] = GetAlignedWaveform(wREC[n],wREC[n]);
413  ebbh->srec[n] = GetAlignedWaveform(sREC[n],wREC[n]);
414  }
415 
416  cout << "segment start time " << seggps << endl;
417  cout << "cluster rate " << crate << endl;
418  cout << "cluster size " << cluster->size() << endl;
419 
420  if(cluster->size()==0) {cout << "ReadDataFromROOT - Error : no clusters in the root file !!!" << endl;return 1;}
421 
422  ebbh->nc.rate=crate;
423  std::vector<int> pList;
424  for(int i=0;i<cluster->size();i++) {
425  ebbh->nc.pList.push_back((*cluster)[i]);
426  pList.push_back(i);
427  }
428  ebbh->nc.cList.push_back(pList);
429  clusterdata cd; // dummy cluster data
430  ebbh->nc.cData.push_back(cd);
431 
432  ebbh->seggps = seggps;
433 
434  ebbh->mchirp = pow(ebbh->mass[0]*ebbh->mass[1],3./5.)/pow(ebbh->mass[0]+ebbh->mass[1],1./5.);
435 
436  char name[1024];
437  sprintf(name,"source : m1=%0.1f m2=%0.1f mchirp=%0.1f e0=%0.1f rp0=%0.1f dist=%0.1f (Mpc) redshift=%0.1f isnr=%0.1f osnr=%0.1f",
438  ebbh->mass[0], ebbh->mass[1], ebbh->mchirp, ebbh->e0, ebbh->rp0, ebbh->dist, ebbh->redshift, ebbh->isnr, ebbh->osnr);
439  ebbh->name = name;
440 
441  for(int n=0;n<nIFO;n++) {
442  if(wDAT[n]) delete wDAT[n];
443  if(wINJ[n]) delete wINJ[n];
444  if(wREC[n]) delete wREC[n];
445  if(sREC[n]) delete sREC[n];
446  }
447  if(cluster) delete cluster;
448  if(itree) delete itree;
449  if(froot) delete froot;
450 
451  return 0;
452 }
453 
454 
455 double GetTimeBoundaries(wavearray<double> x, double P, double& bT, double& eT) {
456 
457  if(P<0) P=0;
458  if(P>1) P=1;
459 
460  int N = x.size();
461 
462  double E = 0; // signal energy
463  double avr = 0; // average
464  for(int i=0;i<N;i++) {avr+=i*x[i]*x[i]; E+=x[i]*x[i];}
465  int M=int(avr/E); // central index
466 
467  // search range which contains percentage P of the total energy E
468  int jB=0;
469  int jE=N-1;
470  double a,b;
471  double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
472  for(int j=1; j<N; j++) {
473  a = ((M-j>=0)&&(M-j<N)) ? x[M-j] : 0.;
474  b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
475  if(a) jB=M-j;
476  if(b) jE=M+j;
477  sum += a*a+b*b;
478  if(sum/E > P) break;
479  }
480 
481  bT = x.start()+jB/x.rate();
482  eT = x.start()+jE/x.rate();
483 
484  return eT-bT;
485 }
486 
487 wavearray<double> GetAlignedWaveform(wavearray<double>* wf1, wavearray<double>* wf2) {
488 
489  wavearray<double> wf = *wf2;
490  wf=0;
491 
492  if(wf1==NULL) return wf;
493  if(wf1->size()==0) return wf;
494 
495  double R=wf1->rate();
496 
497  double b_wf1 = wf1->start();
498  double e_wf1 = wf1->start()+wf1->size()/R;
499  double b_wf2 = wf2->start();
500  double e_wf2 = wf2->start()+wf2->size()/R;
501 
502  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
503  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
504 
505  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
506  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
507  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
508 
509  for(int i=0;i<sizeXCOR;i++) wf[i+o_wf2] = wf1->data[i+o_wf1];
510 
511  return wf;
512 }
513 
514 wavearray<double> GetDifWaveform(wavearray<double>* wf1, wavearray<double>* wf2) {
515 
516  double R=wf1->rate();
517 
518  double b_wf1 = wf1->start();
519  double e_wf1 = wf1->start()+wf1->size()/R;
520  double b_wf2 = wf2->start();
521  double e_wf2 = wf2->start()+wf2->size()/R;
522 
523  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
524  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
525 
526  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
527  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
528  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
529 
530  wavearray<double> wfdif(sizeXCOR);
531  wfdif=0.;
532  wfdif.rate(R);
533  wfdif.start(b_wf1+double(o_wf1)/R);
534 
535  for(int i=0;i<sizeXCOR;i++) wfdif[i] = wf1->data[i+o_wf1] - wf2->data[i+o_wf2];
536 
537  return wfdif;
538 }
539 
540 double GetFitParameter(double mchirp) {
541 
542  const double G = watconstants::GravitationalConstant();
543  const double SM = watconstants::SolarMass();
544  const double C = watconstants::SpeedOfLightInVacuo();
545  const double Pi = TMath::Pi();
546 
547 
548  double A = (96./5.) * pow(Pi,8./3.) * pow(G*SM*mchirp/C/C/C, 5./3);
549  double B = 3./8.;
550 
551  double p0 = -A/B;
552 
553  return p0;
554 }
555 
556 void PlotCluster(watplot* WTS, netcluster* pwc, double fLow, double fHigh) {
557 
558  double RATE = pwc->rate; // original rate
559  std::vector<int>* vint = &(pwc->cList[0]); // pixel list
560  int V = vint->size(); // cluster size
561  for(int j=0; j<V; j++) { // loop over the pixels
562  netpixel* pix = pwc->getPixel(1,j);
563  if(!pix->core) continue;
564  if(pix->likelihood==0) {
565  for(int m=0; m<2; m++) {
566  pix->setdata(0,'S',m); // snr whitened reconstructed signal 00
567  pix->setdata(0,'P',m); // snr whitened reconstructed signal 90
568  pix->setdata(0,'W',m); // snr whitened at the detector 00
569  pix->setdata(0,'U',m); // snr whitened at the detector 90
570  }
571  }
572  }
573  WTS->plot(pwc, 1, 2, 'L', 0, const_cast<char*>("COLZ"));
574  WTS->hist2D->GetYaxis()->SetRangeUser(fLow, fHigh);
575  WTS->canvas->SetLogy();
576  return;
577 }
578 
579 void PlotMultiChirp(EBBH* ebbh) {
580 
581  EColor color[NCHIRP_MAX] = {kBlack, kBlue, kGreen, kMagenta};
582 
583  double mchirp[NCHIRP_MAX];
584  double energy[NCHIRP_MAX];
585 
586  canvas = new TCanvas;
587  canvas->SetGridx();
588  canvas->SetGridy();
589 
590  // plot multi chirp likelihood tf
591 
592  double tmerger=0;
593  double xmin=0;
594  double xmax=0;
595 
596  for(int i=0; i<NCHIRP_MAX; i++) {
597  PlotChirp(ebbh, i, tmerger, xmin, xmax, color[i]);
598  mchirp[i] = ebbh->ch_mass[i];
599  energy[i] = ebbh->ch_energy[i];
600  }
601 
602  char title[256];
603  sprintf(title,"chirp mass : rec = %3.3f / %3.3f / %3.3f / %3.3f - chirp_energy = %3.2f / %3.2f / %3.2f / %3.2f",
604  mchirp[0],mchirp[1],mchirp[2],mchirp[3],energy[0],energy[1],energy[2],energy[3]);
605  xgr[0]->SetTitle(title);
606 
607  if(gPLOT) {canvas->Print(gOFILE);exit(0);}
608 
609  return;
610 }
611 
612 void PlotChirp(EBBH* ebbh, int mch_order, double& tmerger, double& xmin, double& xmax, EColor color) {
613 
614  double chi2_thr = CHI2_THR;
615  double tmerger_cut = TMERGER_CUT;
616  double zmax_thr = ZMAX_THR;
617 
618  double echirp = ebbh->nc.mchirp(1,chi2_thr,tmerger_cut,zmax_thr);
619 
620  clusterdata* pcd = &ebbh->nc.cData[0];
621 
622  ebbh->ch_mass[mch_order] = pcd->mchirp;
623  ebbh->ch_merr[mch_order] = pcd->mchirperr;
624  ebbh->ch_mchi2[mch_order] = pcd->chi2chirp;
625  ebbh->ch_energy[mch_order] = echirp;
626 
627  TGraphErrors* gr = &(pcd->chirp);
628 
629  xgr[mch_order] = new TGraphErrors(gr->GetN(),gr->GetX(),gr->GetY(),gr->GetEX(),gr->GetEY());
630 
631  char title[256];
632  sprintf(title,"chirp mass : rec = %3.3f [%3.2f] , chi2 = %3.2f",
633  pcd->mchirp,pcd->mchirperr,pcd->chi2chirp);
634 
635  xgr[mch_order]->SetTitle(title);
636  xgr[mch_order]->GetHistogram()->SetStats(kFALSE);
637  xgr[mch_order]->GetHistogram()->SetTitleFont(12);
638  xgr[mch_order]->SetFillColor(kWhite);
639  xgr[mch_order]->SetLineColor(kBlack);
640  xgr[mch_order]->GetXaxis()->SetNdivisions(506);
641  xgr[mch_order]->GetXaxis()->SetLabelFont(42);
642  xgr[mch_order]->GetXaxis()->SetLabelOffset(0.014);
643  xgr[mch_order]->GetXaxis()->SetTitleOffset(1.4);
644  xgr[mch_order]->GetYaxis()->SetTitleOffset(1.2);
645  xgr[mch_order]->GetYaxis()->SetNdivisions(506);
646  xgr[mch_order]->GetYaxis()->SetLabelFont(42);
647  xgr[mch_order]->GetYaxis()->SetLabelOffset(0.01);
648  xgr[mch_order]->GetXaxis()->SetTitleFont(42);
649  xgr[mch_order]->GetXaxis()->SetTitle("Time (sec)");
650  xgr[mch_order]->GetXaxis()->CenterTitle(true);
651  xgr[mch_order]->GetYaxis()->SetTitleFont(42);
652  xgr[mch_order]->GetYaxis()->SetTitle("Frequency^{-8/3}");
653  xgr[mch_order]->GetYaxis()->CenterTitle(true);
654  xgr[mch_order]->GetXaxis()->SetLabelSize(0.03);
655  xgr[mch_order]->GetYaxis()->SetLabelSize(0.03);
656 
657  xgr[mch_order]->SetLineColor(kGreen);
658  xgr[mch_order]->SetMarkerStyle(20);
659  xgr[mch_order]->SetMarkerColor(color);
660  xgr[mch_order]->SetMarkerSize(1);
661  mch_order==0 ? xgr[mch_order]->Draw("XAP") : xgr[mch_order]->Draw("XPsame");
662 
663  TF1* fit = &(pcd->fit);
664  double mchirp = pcd->mchirp;
665  if(tmerger==0) tmerger = -fit->GetParameter(1)/fit->GetParameter(0);
666  //if(xmin==0) xmin = fit->GetXmin();
667  //if(xmax==0) xmax = tmerger-0.001;
668  xmin = fit->GetXmin();
669  xmax = tmerger-0.001;
670  double p0 = GetFitParameter(mchirp);
671  double p1 = tmerger;
672 
673  xfit[mch_order] = new TF1(fit->GetName(), fit->GetTitle(), xmin, xmax);
674  xfit[mch_order]->SetLineColor(kRed);
675  xfit[mch_order]->SetParameter(0, p0);
676  xfit[mch_order]->SetParameter(1, p1);
677  xgr[mch_order]->Fit(xfit[mch_order],"Q");
678  xfit[mch_order]->Draw("same");
679 
680  ebbh->ch_tmerger[mch_order] = -xfit[mch_order]->GetParameter(1)/xfit[mch_order]->GetParameter(0);;
681 
682  return;
683 }
684 
685 void PlotMChirp(EBBH* ebbh) {
686 
687  double chi2_thr = CHI2_THR;
688  double tmerger_cut = TMERGER_CUT;
689  double zmax_thr = ZMAX_THR;
690 
691  ebbh->nc.cData[0].mchirp = 0;
692  ebbh->nc.cData[0].mchirperr = 0;
693  ebbh->nc.cData[0].tmrgr = 0;
694  ebbh->nc.cData[0].tmrgrerr = 0;
695  ebbh->nc.cData[0].chi2chirp = 0;
696 
697  PCH = new watplot(const_cast<char*>("chirp"));
698 
699  for(int i=0; i<NCHIRP_MAX; i++) {
700 
701  ebbh->nc.mchirp(1,chi2_thr,tmerger_cut,zmax_thr);
702 
703  clusterdata* pcd = &ebbh->nc.cData[0];
704  printf("mchirp : %.2e %.3f %.3f %.3f %.3f \n\n",
705  pcd->mchirp, pcd->mchirperr, pcd->tmrgr,
706  pcd->tmrgrerr, pcd->chi2chirp);
707 
708  // draw chirp : f^(-8/3) vs time
709  if(gPLOT) {
710  TString ofile = gOFILE;
711  ofile.ReplaceAll("gTYPE_6",TString::Format("gTYPE_6_%d",i+1).Data());
712 
713  PCH->canvas->cd();
714  PCH->plot(pcd, 0);
715  PCH->canvas->Print(ofile);
716  } else {
717  cout << "WARNING : Not Implemented : use gtype = -6" << endl;exit(0);
718  char title[256];
719  sprintf(title,"chirp mass : rec = %3.3f [%3.2f] , chi2 = %3.2f",
720  pcd->mchirp,pcd->mchirperr,pcd->chi2chirp);
721  TGraphErrors* gr = &(pcd->chirp);
722  xgr[i] = new TGraphErrors(gr->GetN(),gr->GetX(),gr->GetY(),gr->GetEX(),gr->GetEY());
723  xgr[i]->SetTitle(title);
724  xgr[i]->Draw("ALP");
725  return;
726  TF1* fit = &(pcd->fit);
727  double tmerger = -fit->GetParameter(1)/fit->GetParameter(0);
728  double xmin = fit->GetXmin()-0.1;
729  double xmax = tmerger-0.001;
730  xfit[i] = new TF1(fit->GetName(), "pow([0]*(x-[1]),-3./8.)", xmin, xmax);
731  double p0 = GetFitParameter(pcd->mchirp);
732  double p1 = tmerger;
733  xfit[i]->SetLineColor(kRed);
734  xfit[i]->SetParameter(0, p0);
735  xfit[i]->SetParameter(1, p1);
736  xfit[i]->Draw("same");
737  return;
738  }
739  }
740 
741  exit(0);
742 }
743 
744 void PlotMultiChirp2(EBBH* ebbh, int ifo, TString ptype, TString wtype) {
745 
746  double tmerger=0;
747  double xmin=0;
748  double xmax=0;
749 
750  // plot chirp
751 
752  cout << endl;
753  cout << "-------------------------------------------------------------------------------------------------------------------------------------------" << endl;
754  cout << ebbh->name << endl;
755  cout << endl << " ";
756  cout << "spin1x=" << ebbh->spin[0] << "\tspin1y=" << ebbh->spin[1] << "\tspin1z=" << ebbh->spin[2]
757  << "\tspin2x=" << ebbh->spin[3] << "\tspin2y=" << ebbh->spin[4] << "\tspin2z=" << ebbh->spin[2] << endl;
758  cout << "-------------------------------------------------------------------------------------------------------------------------------------------" << endl;
759 
760  for(int i=0;i<NCHIRP_MAX;i++) FindChirp(ebbh, i, tmerger, xmin, xmax, ifo, ptype, wtype);
761 
762  cout << "-------------------------------------------------------------------------------------------------------------------------------------------" << endl;
763 
764  // print final
765 
766  int imax=0;
767  float emax=0;
768  // find chirp with max energy
769  for(int i=0;i<NCHIRP_MAX;i++) if(ebbh->ch_energy[i]>emax) {emax=ebbh->ch_energy[i];imax=i;}
770  int imax2=imax;
771  // find chirp with highest chirp mass and energy > (max chirp energy)/4
772  for(int i=0;i<NCHIRP_MAX;i++) {
773  if(i==imax) continue;
774  if(ebbh->ch_mass[i]>0) {
775  if(ebbh->ch_mass[i]>ebbh->ch_mass[imax]) if(ebbh->ch_energy[i]>ebbh->ch_energy[imax]/4.) imax2=i;
776  }
777  }
778  imax=imax2;
779  float ebbh->ei=0;
780  for(int i=0;i<NCHIRP_MAX;i++) {
781  if(ebbh->ch_mass[i]>0) {
782  if(ebbh->ch_mass[i]<ebbh->ch_mass[imax]) ebbh->ei+=ebbh->ch_energy[i];
783  if(ebbh->ch_mass[i]>ebbh->ch_mass[imax]) ebbh->ei-=ebbh->ch_energy[i];
784  }
785  }
786  ebbh->ei/=ebbh->ch_energy[imax];
787  ebbh->ei*=100;
788 
789  cout << "CHIRP_M -> " << " MCH_MASS : " << " rec :" << ebbh->ch_mass[imax] << " [" << ebbh->ch_merr[imax] << "]"
790  << " / "<< " inj : " << ebbh->mchirp << "\t\t\t\tMCH_EI : " << ebbh->ei << "\t\tMCH_TIME : " << ebbh->ch_tmerger[imax] << " (sec)" << endl;
791  cout << "-------------------------------------------------------------------------------------------------------------------------------------------" << endl;
792 
793  return;
794 }
795 
796 void FindChirp(EBBH* ebbh, int mch_order, double& tmerger, double& xmin, double& xmax, int ifo, TString ptype, TString wtype) {
797 
798  if((ptype!="like")&&(ptype!="stft")&&(ptype!="")) {
799  cout << "PlotMultiChirp2 Error : wrong ptype, available value are : like/stft" << endl;
800  exit(1);
801  }
802 
803  double chi2_thr = CHI2_THR;
804  double tmerger_cut = TMERGER_CUT;
805  double zmax_thr = ZMAX_THR;
806 
807  EColor lcolor = (ptype=="like") ? kBlack : kWhite;
808  int lwidth=3;
809  int lstyle=2;
810 
811  if(mch_order==0 && ptype=="like") {
812  WTS = new watplot(const_cast<char*>("wts"));
813  WTS->plot(&ebbh->nc, 1, 2, 'L', 0, const_cast<char*>("COLZ"));
814  WTS->hist2D->GetYaxis()->SetRangeUser(FLOW, FHIGH);
815  WTS->canvas->SetLogy();
816  }
817 
818  double echirp = ebbh->nc.mchirp(1,chi2_thr,tmerger_cut,zmax_thr);
819 
820  if(mch_order==0 && ptype=="like") {
821  if(gTYPE==3) PlotCluster(WTS, &ebbh->nc, FLOW, FHIGH);
822  }
823 
824  clusterdata* pcd = &ebbh->nc.cData[0];
825  TF1* fit = &(pcd->fit);
826 
827  double mchirp = pcd->mchirp;
828 
829  if(tmerger==0) tmerger = -fit->GetParameter(1)/fit->GetParameter(0);
830  double Xmin = fit->GetXmin()-0.1;
831  double Xmax = tmerger-0.001;
832  if(mch_order==0) {xmin=Xmin;xmax=Xmax;}
833  if(Xmin<xmin) xmin = Xmin;
834  if(Xmax>xmax) xmax = Xmax;
835 
836  double p0 = GetFitParameter(mchirp);
837  double p1 = tmerger;
838 
839  ebbh->ch_mass[mch_order] = pcd->mchirp;
840  ebbh->ch_tmerger[mch_order] = -fit->GetParameter(1)/fit->GetParameter(0);
841  ebbh->ch_merr[mch_order] = pcd->mchirperr;
842  ebbh->ch_mchi2[mch_order] = pcd->chi2chirp;
843  ebbh->ch_energy[mch_order] = echirp;
844 
845  // fit function -> freq = p0 * (time-p1);
846  xfit[mch_order] = new TF1(fit->GetName(), "pow([0]*(x-[1]),-3./8.)", xmin, xmax);
847  xfit[mch_order]->SetLineColor(lcolor);
848  xfit[mch_order]->SetLineWidth(lwidth);
849  xfit[mch_order]->SetLineStyle(lstyle);
850  xfit[mch_order]->SetParameter(0, p0);
851  xfit[mch_order]->SetParameter(1, p1);
852  //mch_order==0 ? xfit[mch_order]->SetParameter(1, p1) : xfit[mch_order]->FixParameter(1, p1);
853 
854  if(mch_order==0 && ptype=="stft") {
855  if(wtype=="wdat") MDC->Draw(ebbh->wdat[ifo],MDC_TF);
856  if(wtype=="winj") MDC->Draw(ebbh->winj[ifo],MDC_TF);
857  if(wtype=="wrec") MDC->Draw(ebbh->wrec[ifo],MDC_TF);
858  stft = MDC->GetSTFT();
859  stft->GetHistogram()->SetTitle(ebbh->name);
860  stft->GetHistogram()->GetYaxis()->SetRangeUser(FLOW,FHIGH);
861  }
862  if(stft) stft->GetHistogram()->GetXaxis()->SetRangeUser(xmin,xmax+0.3);
863 
864  if(ptype!="") xfit[mch_order]->Draw("same");
865 
866  if(gPLOT && mch_order==NCHIRP_MAX-1 && ptype=="stft") {stft->Print(gOFILE);exit(0);}
867  if(gPLOT && mch_order==NCHIRP_MAX-1 && ptype=="like") {WTS->canvas->Print(gOFILE);exit(0);}
868 
869  // format standard output
870  cout.precision(1);
871  cout.setf(ios::fixed, ios::floatfield);
872 
873  cout << "CHIRP_" << mch_order+1 << " -> "
874  << " MCH_MASS : " << ebbh->ch_mass[mch_order] << "\t\tMCH_ERR : " << ebbh->ch_merr[mch_order]
875  << "\t\tMCH_CHI2 : " << ebbh->ch_mchi2[mch_order] << "\t\tMCH_EN : " << ebbh->ch_energy[mch_order]
876  << "\t\tMCH_TIME : " << setprecision(2) << ebbh->ch_tmerger[mch_order] << " (sec)" << endl;
877 
878  return;
879 }
880 
881 void FillHistogram(EBBH* ebbh, TH2F&* hN, TH2F&* hL, TH2F&* heT, TH2F&* heF, double zmax_thr) {
882 
883  netcluster* pwc = &(ebbh->nc);
884 
885  int cid = 1;
886  int nifo = ebbh->ifo.size();
887 
888  double RATE = pwc->rate; // original rate
889 
890  std::vector<int>* vint = &(pwc->cList[cid-1]); // pixel list
891 
892  int V = vint->size(); // cluster size
893  if(!V) return;
894 
895  double zmax = -1e100;
896  double etot = 0;
897  for(int j=0; j<V; j++) {
898  netpixel* pix = pwc->getPixel(cid, j);
899  if(pix->likelihood<1. || pix->frequency==0 ) continue;
900  if(pix->likelihood>zmax) zmax = pix->likelihood;
901  etot+=pix->likelihood;
902  }
903  double zthr = zmax_thr*zmax;
904 
905  int minLayers=1000;
906  int maxLayers=0;
907  double minTime=1e20;
908  double maxTime=0.;
909  double minFreq=1e20;
910  double maxFreq=0.;
911  for(int j=0; j<V; j++) { // loop over the pixels
912  netpixel* pix = pwc->getPixel(cid,j);
913  if(!pix->core) continue;
914 
915  if(pix->layers<minLayers) minLayers=pix->layers;
916  if(pix->layers>maxLayers) maxLayers=pix->layers;
917 
918  double dt = 1./pix->rate;
919  double time = int(pix->time/pix->layers)/double(pix->rate); // central bin time
920  time -= dt/2.; // begin bin time
921  if(time<minTime) minTime=time;
922  if(time+dt>maxTime) maxTime=time+dt;
923 
924  double freq = pix->frequency*pix->rate/2.;
925  if(freq<minFreq) minFreq=freq;
926  if(freq>maxFreq) maxFreq=freq;
927  }
928 
929  int minRate=RATE/(maxLayers-1);
930  int maxRate=RATE/(minLayers-1);
931 
932  double mindt = 1./maxRate;
933  double maxdt = 1./minRate;
934  double mindf = minRate/2.;
935  double maxdf = maxRate/2.;
936 
937  cout.precision(1);
938  cout.setf(ios::fixed, ios::floatfield);
939 
940 /*
941  cout << "minRate : " << minRate << "\t\t\t maxRate : " << maxRate << endl;
942  cout << "minTime : " << minTime << "\t\t\t maxTime : " << maxTime << endl;
943  cout << "minFreq : " << minFreq << "\t\t\t maxFreq : " << maxFreq << endl;
944  cout << "mindt : " << mindt << "\t\t\t maxdt : " << maxdt << endl;
945  cout << "mindf : " << mindf << "\t\t\t maxdf : " << maxdf << endl;
946 */
947 
948  double iminTime = minTime-maxdt;
949  double imaxTime = maxTime+maxdt;
950  int nTime = (imaxTime-iminTime)*maxRate;
951 
952  if(hN) { delete hN; hN=NULL; }
953  hN = new TH2F("hN", "hN", nTime, iminTime, imaxTime, 2*(maxLayers-1), 0, RATE/2);
954  if(hL) { delete hL; hL=NULL; }
955  hL = new TH2F("hL", "hL", nTime, iminTime, imaxTime, 2*(maxLayers-1), 0, RATE/2);
956  if(heT) { delete heT; heT=NULL; }
957  heT = new TH2F("heT", "heT", nTime, iminTime, imaxTime, 2*(maxLayers-1), 0, RATE/2);
958  if(heF) { delete heF; heF=NULL; }
959  heF = new TH2F("heF", "heF", nTime, iminTime, imaxTime, 2*(maxLayers-1), 0, RATE/2);
960 
961  double dFreq = (maxFreq-minFreq)/10.>2*maxdf ? (maxFreq-minFreq)/10. : 2*maxdf ;
962  double mFreq = minFreq-dFreq<0 ? 0 : minFreq-dFreq;
963  double MFreq = maxFreq+dFreq>RATE/2 ? RATE/2 : maxFreq+dFreq;
964 
965  double dTime = (maxTime-minTime)/10.>2*maxdt ? (maxTime-minTime)/10. : 2*maxdt ;
966  double mTime = minTime-dTime<iminTime ? iminTime : minTime-dTime;
967  double MTime = maxTime+dTime>imaxTime ? imaxTime : maxTime+dTime;
968 
969  int npix=0;
970  double Likelihood=0;
971  for(int n=0; n<V; n++) {
972  netpixel* pix = pwc->getPixel(cid,n);
973  if(!pix->core) continue;
974 
975  if(pix->likelihood<zthr) continue;
976 
977  double sSNR=0;
978  for(int m=0; m<nifo; m++) {
979  sSNR += pow(pix->getdata('S',m),2); // snr whitened reconstructed signal 00
980  sSNR += pow(pix->getdata('P',m),2); // snr whitened reconstructed signal 90
981  }
982 
983  int iRATE = int(pix->rate+0.5);
984  int M=maxRate/iRATE;
985  int K=2*(maxLayers-1)/(pix->layers-1);
986  double dt = 1./pix->rate;
987  double itime = int(pix->time/pix->layers)/double(pix->rate); // central bin time
988  itime -= dt/2.; // begin bin time
989  int i=(itime-iminTime)*maxRate;
990  int j=pix->frequency*K;
991  Likelihood+=sSNR;
992 // sSNR /= M*K; // Normalization (watplot not use normalization !!!)
993 //sSNR = pix->likelihood; // Use original oriinal likelihood before WP
994  for(int m=0;m<M;m++) {
995  for(int k=0;k<K;k++) {
996  double L = hL->GetBinContent(i+1+m,j+1+k-K/2);
997  hL->SetBinContent(i+1+m,j+1+k-K/2,sSNR+L);
998 
999  double N = hN->GetBinContent(i+1+m,j+1+k-K/2);
1000  hN->SetBinContent(i+1+m,j+1+k-K/2,N+1);
1001 
1002  double eT = heT->GetBinContent(i+1+m,j+1+k-K/2);
1003  //heT->SetBinContent(i+1+m,j+1+k-K/2,eT+M*M*L);
1004  heT->SetBinContent(i+1+m,j+1+k-K/2,eT+M);
1005 
1006  double eF = heF->GetBinContent(i+1+m,j+1+k-K/2);
1007  //heF->SetBinContent(i+1+m,j+1+k-K/2,eF+K*K*L);
1008  heF->SetBinContent(i+1+m,j+1+k-K/2,eF+K);
1009  }
1010  }
1011  npix++;
1012  }
1013  cout << "Likelihood : " << Likelihood << " npix : " << npix << endl;
1014 
1015  for(int i=0;i<=hL->GetNbinsX();i++) {
1016  for(int j=0;j<=hL->GetNbinsY();j++) {
1017 
1018  double L = hL->GetBinContent(i,j);
1019  if(L<=0) continue;
1020 
1021  double N = hN->GetBinContent(i,j);
1022 
1023  double eT = heT->GetBinContent(i,j);
1024  //heT->SetBinContent(i,j,mindt*sqrt(eT/L));
1025  //heT->SetBinContent(i,j,mindt);
1026  heT->SetBinContent(i,j,mindt*eT/N);
1027 
1028  double eF = heF->GetBinContent(i,j);
1029  //heF->SetBinContent(i,j,(mindf/2.)*sqrt(eF/L));
1030  //heF->SetBinContent(i,j,(mindf/2.));
1031  heF->SetBinContent(i,j,(mindf/2.)*eF/N);
1032  }
1033  }
1034 
1035  return;
1036 }
1037 
1038 void PlotHistogram(EBBH* ebbh, TH2F&* h2d) {
1039 
1040  if(!h2d) return;
1041 
1042  if(!canvas) canvas = new TCanvas;
1043  canvas->SetGridx();
1044  canvas->SetGridy();
1045  canvas->SetLogy();
1046 
1047  h2d->SetTitle(ebbh->name);
1048 
1049  h2d->GetYaxis()->SetRangeUser(FLOW,FHIGH);
1050 
1051  h2d->SetStats(kFALSE);
1052  h2d->SetTitleFont(12);
1053  h2d->SetFillColor(kWhite);
1054 
1055  h2d->GetXaxis()->SetNdivisions(506);
1056  h2d->GetXaxis()->SetLabelFont(42);
1057  h2d->GetXaxis()->SetLabelOffset(0.014);
1058  h2d->GetXaxis()->SetTitleOffset(1.4);
1059  h2d->GetYaxis()->SetTitleOffset(1.2);
1060  h2d->GetYaxis()->SetNdivisions(506);
1061  h2d->GetYaxis()->SetLabelFont(42);
1062  h2d->GetYaxis()->SetLabelOffset(0.01);
1063  h2d->GetZaxis()->SetLabelFont(42);
1064  h2d->GetZaxis()->SetNoExponent(false);
1065  h2d->GetZaxis()->SetNdivisions(506);
1066 
1067  h2d->GetXaxis()->SetTitleFont(42);
1068  h2d->GetXaxis()->SetTitle("Time (sec)");
1069  h2d->GetXaxis()->CenterTitle(true);
1070  h2d->GetYaxis()->SetTitleFont(42);
1071  h2d->GetYaxis()->SetTitle("Frequency (Hz)");
1072  h2d->GetYaxis()->CenterTitle(true);
1073 
1074  h2d->GetZaxis()->SetTitleOffset(0.6);
1075  h2d->GetZaxis()->SetTitleFont(42);
1076  h2d->GetZaxis()->CenterTitle(true);
1077 
1078  h2d->GetXaxis()->SetLabelSize(0.03);
1079  h2d->GetYaxis()->SetLabelSize(0.03);
1080  h2d->GetZaxis()->SetLabelSize(0.03);
1081 
1082  h2d->Draw("colz");
1083  if(fchirp) fchirp->Draw("same");
1084 
1085  if(gPLOT) {canvas->Print(gOFILE);exit(0);}
1086 
1087  return;
1088 }
1089 
1090 double ChirpIntegralHistogram(EBBH* ebbh, TH2F&* h2d, double dMc) {
1091 
1092  if(!h2d) return;
1093 
1094  const double G = watconstants::GravitationalConstant();
1095  const double SM = watconstants::SolarMass();
1096  const double C = watconstants::SpeedOfLightInVacuo();
1097  const double Pi = TMath::Pi();
1098 
1099  double Mc = ebbh->ch_mass[0];
1100  //cout << "Estimated Mc : " << Mc << endl;
1101 
1102  Mc+=dMc;
1103 
1104  double p0 = 256.*Pi/5*pow(G*Mc*SM*Pi/C/C/C, 5./3);
1105  double p1 = ebbh->ch_tmerger[0];
1106 
1107  double xmin = h2d->GetXaxis()->GetXmin();
1108  double xmax = h2d->GetXaxis()->GetXmax();
1109  if(p1<xmax) xmax=p1;
1110  cout.precision(3);
1111  //cout << "xmin : " << xmin << " xmax : " << xmax << endl;
1112 
1113  // fit function -> freq = p0 * (time-p1);
1114  fchirp = new TF1("fchirp", "pow([0]*([1]-x),-3./8.)", xmin, xmax);
1115  fchirp->SetLineColor(kBlack);
1116  fchirp->SetLineWidth(3);
1117  fchirp->SetLineStyle(2);
1118  fchirp->SetParameter(0, p0);
1119  fchirp->SetParameter(1, p1);
1120 
1121  double dT = hL->GetXaxis()->GetBinWidth(0);
1122  double dF = hL->GetYaxis()->GetBinWidth(0);
1123  //cout << "dT = " << dT << " dF = " << dF << endl;
1124 
1125  double cenergy=0.;
1126  for(int i=1;i<=hL->GetNbinsX();i++) {
1127  double T = hL->GetXaxis()->GetBinCenter(i)+hL->GetXaxis()->GetBinWidth(i)/2.;
1128  double DF = fchirp->Eval(T)-fchirp->Eval(T-dT);
1129  for(int j=0;j<=hL->GetNbinsY();j++) {
1130 
1131  double L = hL->GetBinContent(i,j);
1132  if(L<=0) continue;
1133 
1134  double F = hL->GetYaxis()->GetBinCenter(j)+hL->GetYaxis()->GetBinWidth(j)/2.;
1135 
1136  //if(F-fchirp->Eval(T-dT)<DF && (F-fchirp->Eval(T-dT)>0)) hL->SetBinContent(i,j,0);
1137  if(F-fchirp->Eval(T-dT)<DF && (F-fchirp->Eval(T-dT)>0)) cenergy+=L;
1138  }
1139  }
1140 
1141  //cout << "CHIRP ENERGY : " << cenergy << endl;
1142 
1143  delete fchirp; fchirp=NULL;
1144 
1145  return cenergy;
1146 }
1147 
1148 void DrawMChirpInstantaneousFrequency(EBBH* ebbh, TString wtype, int ifo) {
1149 
1150  wavearray<double> wf;
1151 
1152  if(wtype=="winj") wf = ebbh->winj[ifo]; // whiten injected waveforms
1153  if(wtype=="sinj") wf = ebbh->sinj[ifo]; // strain injected waveforms
1154  if(wtype=="wrec") wf = ebbh->wrec[ifo]; // whiten reconstructed waveforms
1155  if(wtype=="srec") wf = ebbh->srec[ifo]; // strain injected waveforms
1156  if(wtype=="wsim") wf = ebbh->wsim[0]; // simulated waveforms -> hp
1157 
1158  double start = wf.start();
1159  double imchirp = ebbh->ch_mass[0];
1160  //double tmerger = ebbh->ch_tmerger[0]-start;
1161  double tmerger = ebbh->ch_tmerger[0];
1162  cout << "mchirp - mchirp : " << imchirp << endl;
1163  cout << "mchirp - tmerger : " << tmerger << endl;
1164 
1165  //wf.start(0);
1166 
1167  // Get hp envelope
1168  wavearray<double> ewf = CWB::Toolbox::getHilbertEnvelope(wf);
1169  MDC->SetZoom(0.999);
1170  MDC->Draw(wf,MDC_TIME,"ALP ZOOM",kGray);
1171  MDC->Draw(ewf,MDC_TIME,"same",kRed);
1172 
1173  watplot* pts = MDC->GetWatPlot();
1174  TGraph* gr = pts->getGraph(0);
1175  gr->SetTitle(ebbh->name);
1176 
1177  // Get wf instantaneous frequency
1178  wavearray<double> fwf = CWB::Toolbox::getHilbertIFrequency(wf);
1179  //MDC->Draw(fwf,MDC_TIME,"ALP ZOOM",kBlack);return;
1180 
1181  // F -> F^-3/8
1182  double sF = 128;
1183  for(int i=0;i<fwf.size();i++) {
1184  fwf[i] = fwf[i]>FLOW ? pow(fwf[i]/sF,-8./3.) : 0;
1185  }
1186  //MDC->Draw(fwf,MDC_TIME,"ALP ZOOM",kBlack);return;
1187 
1188  const double G = watconstants::GravitationalConstant();
1189  const double SM = watconstants::SolarMass();
1190  const double C = watconstants::SpeedOfLightInVacuo();
1191  const double Pi = TMath::Pi();
1192 
1193  double M1=ebbh->mass[0];
1194  double M2=ebbh->mass[1];
1195  double Mc = pow(M1*M2,3./5.)/pow(M1+M2,1./5.);
1196  double MC = Mc;
1197  double mc = Mc;
1198  cout << "Mc " << Mc << endl;
1199  Mc *= SM;
1200  double kk = 256.*Pi/5*pow(G*SM*Pi/C/C/C, 5./3);
1201  double K = 256.*Pi/5*pow(G*Mc*Pi/C/C/C, 5./3);
1202  // scaling of frequency: in units of 128Hz
1203  kk *= pow(sF, 8./3);
1204  K *= pow(sF, 8./3);
1205 
1206  wavearray<double> tt = fwf;
1207  wavearray<double> ff = fwf;
1208  tt=0;
1209  ff=0;
1210 
1211  double dt = 1./ff.rate();
1212  double Tc = tmerger;
1213  for(int i=0;i<ff.size();i++) {
1214  double T = ff.start()+i*dt;
1215  //double F = pow(Tc-K*T,-3./8.);
1216  double F = T<Tc ? pow(K*(Tc-T),-3./8.) : 0;
1217  //double F = pow(K*(T),-3./8.);
1218  if(F<0 || F>1000) F=0;
1219  F = F>0 ? pow(F,-8./3.) : 0;
1220  tt.data[i] = T;
1221  ff.data[i] = F;
1222  if(F>1) F=0;
1223  }
1224  //MDC->Draw(ff,MDC_TIME,"ALP",kBlack);return;
1225 
1226  // remove data where the waveform amplitude (ewf) is negligeable
1227  double emax=ewf.max();
1228  int SS = 0;
1229  for(int i=0;i<fwf.size();i++) {
1230  if(ewf[i]/emax<0.1) SS++; else break;
1231  }
1232  int EE=0;
1233  for(int i=fwf.size()-1;i>=0;i--) {
1234  if(ewf[i]/emax<0.2) EE++; else break;
1235  }
1236 
1237  // ff errors
1238  wavearray<double> eff = ff;
1239  //eff=1.e20;
1240  //for(int i=0;i<fwf.size();i++) eff[i] = ewf[i]/emax>0.01 ? 1./ewf[i] : 1.e20;
1241  eff=0;
1242 
1243  TGraphErrors *grr = new TGraphErrors(fwf.size()-SS-EE, tt.data+SS, fwf.data+SS, 0, eff.data+SS);
1244  grr->GetHistogram()->SetTitle(ebbh->name);
1245  grr->GetHistogram()->GetXaxis()->SetTitle("time (sec)");
1246  grr->GetHistogram()->GetYaxis()->SetTitle("freq^-3/8 ");
1247  grr->GetHistogram()->GetXaxis()->SetTitleOffset(1.1);
1248  grr->GetHistogram()->GetYaxis()->SetTitleOffset(1.1);
1249  grr->SetMarkerStyle(7);
1250  grr->SetMarkerColor(kGray);
1251  grr->SetLineColor(kGray);
1252 
1253  TF1 *func = new TF1("func", "[0]*(x-[1])", 0.05,Tc);
1254  //func->FixParameter(0,-kk*pow(MC,1./0.6));
1255  //func->FixParameter(0,-kk*pow(imchirp,1./0.6));
1256  //func->FixParameter(1,Tc);
1257  func->SetParameter(0,-kk*pow(imchirp,1./0.6));
1258  func->SetParameter(1,Tc);
1259  func->SetLineColor(kRed);
1260  TCanvas *canvas = new TCanvas("canvas", "Linear and robust linear fitting");
1261  canvas->SetGrid();
1262  grr->Draw("ALP");
1263  printf("Ordinary least squares:\n");
1264  grr->Fit(func);
1265  func->Draw("same");
1266 
1267  double mch = -func->GetParameter(0)/kk;
1268  double mchirp = mch>0 ? pow(mch, 0.6) : -pow(-mch, 0.6);
1269 
1270  cout << "mchirp : " << mchirp << endl;
1271  cout << "mchirp-MC : " << mchirp-MC << endl;
1272  cout << "m1 m2 ichirp ochirp ochirp-ichirp : " << M1 << " " << M2 << " " << MC << " " << " " << mchirp << " " << mchirp-MC << endl;
1273 
1274 
1275  double mm = (pow(80.,-8./3.)-pow(30.,-8./3.))/(0.13);
1276  double mch1 = -mm/kk;
1277  double mchirp1 = mch1>0 ? pow(mch1, 0.6) : -pow(-mch1, 0.6);
1278  cout << "mchirp1 : " << mchirp1 << endl;
1279 
1280  return;
1281 }
1282 
1283 void GetCBC(EBBH* ebbh) {
1284 
1285  double To = 931072130;
1286 
1287  double Tc=ebbh->ch_tmerger[0];
1288 
1289  double M1 = ebbh->mass[0]<ebbh->mass[1] ? ebbh->mass[0] : ebbh->mass[1];
1290  double M2 = ebbh->mass[0]>ebbh->mass[1] ? ebbh->mass[0] : ebbh->mass[1];
1291 
1292  double S1 = fabs(ebbh->spin[0])<fabs(ebbh->spin[1]) ? fabs(ebbh->spin[0]) : fabs(ebbh->spin[1]);
1293  double S2 = fabs(ebbh->spin[0])>fabs(ebbh->spin[1]) ? fabs(ebbh->spin[0]) : fabs(ebbh->spin[1]);
1294 
1295  double DIST=ebbh->dist;
1296 
1297  // ---------------------------------
1298  // set inspiral parms
1299  // ---------------------------------
1300  TString inspOptions="";
1301  inspOptions+= TString::Format("--gps-start-time %f --gps-end-time %f ",Tc+To,Tc+To);
1302  inspOptions+= "--taper-injection start --seed 962488 ";
1303  inspOptions+= "--dir ./ ";
1304  inspOptions+= "--f-lower 24.000000 ";
1305  inspOptions+= "--time-step 60 ";
1306 
1307  inspOptions+= "--waveform IMRPhenomBthreePointFivePN ";
1308  inspOptions+= "--l-distr random ";
1309  inspOptions+= "--i-distr uniform ";
1310 
1311  inspOptions+= TString::Format("--min-mass1 %f --max-mass1 %f ",M1,M1);
1312  inspOptions+= TString::Format("--min-mass2 %f --max-mass2 %f ",M2,M2);
1313  inspOptions+= TString::Format("--min-mtotal %f --max-mtotal %f ",M1+M1,M2+M2);
1314 
1315  inspOptions+= "--m-distr componentMass ";
1316  inspOptions+= "--enable-spin --aligned ";
1317 
1318  inspOptions+= TString::Format("--min-spin1 %f --max-spin1 %f ",S1,S1);
1319  inspOptions+= TString::Format("--min-spin2 %f --max-spin2 %f ",S2,S2);
1320 
1321  inspOptions+= "--amp-order 0 ";
1322  inspOptions+= "--d-distr volume ";
1323  inspOptions+= TString::Format("--min-distance %f --max-distance %f ",DIST,DIST);
1324  //inspOptions+= "--output inspirals.xml "; // set output xml file
1325 
1326  MDC->SetInspiral("IMRPhenomBthreePointFivePN",inspOptions.Data());
1327 
1328  // Get the first waveform hp,hx components starting from gps = 931072130
1329  ebbh->wsim[0] = MDC->GetInspiral("hp",Tc+To-200,Tc+To+200);
1330  ebbh->wsim[1] = MDC->GetInspiral("hx",Tc-200,Tc+200);
1331  ebbh->wsim[0].start(ebbh->wsim[0].start()-To);
1332  ebbh->wsim[0].resample(2*FHIGH);
1333  ebbh->wsim[1].start(ebbh->wsim[1].start()-To);
1334  ebbh->wsim[1].resample(2*FHIGH);
1335  cout << "size : " << ebbh->wsim[0].size() << " rate : " << ebbh->wsim[0].rate() << " start : " << (int)ebbh->wsim[0].start() << endl;
1336 
1337  return;
1338 }
1339 
1340 void GetEBBH(EBBH* ebbh) {
1341 
1342  double To = ebbh->wrec[0].start();
1343 
1344  double Tc=ebbh->ch_tmerger[0];
1345 
1346  double M1 = ebbh->mass[0]<ebbh->mass[1] ? ebbh->mass[0] : ebbh->mass[1];
1347  double M2 = ebbh->mass[0]>ebbh->mass[1] ? ebbh->mass[0] : ebbh->mass[1];
1348 
1349  double S1 = fabs(ebbh->spin[0])<fabs(ebbh->spin[1]) ? fabs(ebbh->spin[0]) : fabs(ebbh->spin[1]);
1350  double S2 = fabs(ebbh->spin[0])>fabs(ebbh->spin[1]) ? fabs(ebbh->spin[0]) : fabs(ebbh->spin[1]);
1351 
1352  double DIST=ebbh->dist;
1353 
1354  double E0 = ebbh->e0;
1355  double RP0 = ebbh->rp0>00 ? ebbh->rp0 : 14;
1356  double REDSHIFT = ebbh->redshift;
1357 
1358  // get ebbh parameters from input command line
1359  TString ebbh_parms = TString::Format("1 %f %f %f %f %f %f",M1,M2,RP0,E0,DIST,REDSHIFT);
1360  cout << "ebbh_parms : " << ebbh_parms << endl;
1361 
1362  int gps = To+Tc;
1363 
1364  vector<mdcpar> par;
1365  par.resize(1);
1366  par[0].name=ebbh_parms;
1367  MDC->AddWaveform(MDC_EBBH, par);
1368  MDC->Print(0);
1369 
1370  // Get the first waveform hp,hx components starting from gps = 931072130
1371 
1372  waveform wf = MDC->GetWaveform("eBBH");
1373  ebbh->wsim[0] = wf.hp;
1374  ebbh->wsim[1] = wf.hx;
1375  ebbh->wsim[0].start(ebbh->wrec[0].start()-ebbh->seggps);
1376  ebbh->wsim[0].resample(2*FHIGH);
1377  ebbh->wsim[1].start(ebbh->wrec[1].start()-ebbh->seggps);
1378  ebbh->wsim[1].resample(2*FHIGH);
1379  cout << "size : " << ebbh->wsim[0].size() << " rate : " << ebbh->wsim[0].rate() << " start : " << (int)ebbh->wsim[0].start() << endl;
1380 
1381  return;
1382 }
spin
float spin[6]
Definition: cbc_plots.C:438
detector::getDetectorParams
detectorParams getDetectorParams()
Definition: detector.cc:201
ZMAX_THR
#define ZMAX_THR
Definition: ChirpMass.C:101
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
MDC_TIME
Definition: mdc.hh:172
C
static const double C
Definition: GNGen.cc:10
ChirpMass
int ChirpMass(TString ifroot, int gtype=0, int entry=0, int ifo=0, EBBH *xebbh=NULL)
Definition: ChirpMass.C:130
gPLOT
bool gPLOT
Definition: ChirpMass.C:111
CWB::Toolbox::getFileListFromDir
static vector< TString > getFileListFromDir(TString dir_name, TString endString="", TString beginString="", TString containString="", bool fast=false)
Definition: Toolbox.cc:4333
Init
void Init()
Definition: ChirpMass.C:280
CWB::mdc::Draw
watplot * Draw(TString name, int id=0, TString polarization="hp", MDC_DRAW type=MDC_TIME, TString options="ALP", Color_t color=kBlack)
Definition: mdc.cc:2288
wREC
std::vector< wavearray< double > > wREC[MAX_TRIALS]
Definition: CWB_Plugin_PE.C:406
waveform
Definition: mdc.hh:189
printf
printf("total live time: non-zero lags = %10.1f \n", liveTot)
FillHistogram
void FillHistogram(EBBH *ebbh, TH2F &*hN, TH2F &*hL, TH2F &*heT, TH2F &*heF, double zmax_thr=0.)
Definition: ChirpMass.C:881
dist
double dist
Definition: ConvertGWGC.C:48
fHigh
double fHigh
Definition: cwb2G_parameters.C:20
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:123
CWB::STFT::Print
void Print(TString pname)
Definition: STFT.cc:297
cluster
std::vector< pixel > cluster
Definition: wavepath.hh:43
ptype
TString ptype[nPLOT]
Definition: cwb_mkfad.C:94
hL
TH2F * hL
Definition: ChirpMass.C:121
xfit
TF1 * xfit[NCHIRP_MAX]
Definition: ChirpMass.C:125
a
wavearray< double > a(hp.size())
netpixel::frequency
size_t frequency
Definition: netpixel.hh:93
netpixel::likelihood
float likelihood
Definition: netpixel.hh:96
name
par[0] name
Definition: CWB_Plugin_HEN_SIM_Config.C:87
B
#define B
n
int n
Definition: cwb_net.C:10
CWB::mdc::GetSTFT
CWB::STFT * GetSTFT()
Definition: mdc.hh:400
MDC_FFT
Definition: mdc.hh:173
NCHIRP_MAX
#define NCHIRP_MAX
Definition: ChirpMass.C:4
FHIGH
#define FHIGH
Definition: ChirpMass.C:87
const_cast< char * >
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
TString
TString("c")
Definition: cwb_report_skymap.C:111
clusterdata::mchirp
float mchirp
Definition: netcluster.hh:66
clusterdata::tmrgrerr
float tmrgrerr
Definition: netcluster.hh:65
x
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
Definition: BaseFunctionWDM.C:51
watconstants::SolarMass
double SolarMass()
Definition: constants.hh:184
PlotChirp
void PlotChirp(EBBH *ebbh, int mch_order, double &tmerger, double &xmin, double &xmax, EColor color)
Definition: ChirpMass.C:612
pix
netpixel pix(nifo)
pwc
netcluster * pwc
Definition: cwb_job_obj.C:20
mdc
TString mdc[4]
Definition: DrawCoverageVsPercentagePRC.C:13
watplot::graph
std::vector< TGraph * > graph
Definition: watplot.hh:176
PlotMultiChirp
void PlotMultiChirp(EBBH *ebbh)
Definition: ChirpMass.C:579
P
TRandom3 P
Definition: compare_bkg.C:296
waveform::hp
wavearray< double > hp
Definition: mdc.hh:196
FLOW
#define FLOW
Definition: ChirpMass.C:86
mchirp
float mchirp
Definition: cbc_plots.C:436
CWB::Toolbox::getHilbertIFrequency
static wavearray< double > getHilbertIFrequency(wavearray< double > x)
Definition: Toolbox.cc:6151
wf
waveform wf
Definition: CWB_Plugin_HEN_SIM_Config.C:82
clusterdata::tmrgr
float tmrgr
Definition: netcluster.hh:64
p1
watplot p1(const_cast< char * >("TFMap1"))
M
#define M
Definition: UniqSLagsList.C:3
netpixel::layers
size_t layers
Definition: netpixel.hh:94
m
int m
Definition: cwb_net.C:10
oSNR
std::vector< double > oSNR[NIFO_MAX]
Definition: CWB_Plugin_PE.C:409
netcluster::cList
std::vector< vector_int > cList
Definition: netcluster.hh:379
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:119
ofile
char ofile[512]
Definition: cwb_condor_rescue.C:40
j
int j
Definition: cwb_net.C:10
cd
canvas cd(1)
i
i drho i
Definition: cwb_epparameters.C:85
netcluster::rate
double rate
Definition: netcluster.hh:360
list
TList * list
Definition: cwb_dump_events.C:31
Clear
void Clear()
Definition: ChirpMass.C:302
CWB::Toolbox::checkFile
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:3956
watplot::plot
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:132
watplot::getGraph
TGraph * getGraph(int n)
Definition: watplot.hh:168
N
#define N
cwb_online.ff
tuple ff
Definition: cwb_online.py:394
CWB::mdc::AddWaveform
mdcid AddWaveform(MDC_TYPE mdc_type, vector< mdcpar > par, TString uname="")
Definition: mdc.cc:445
CWB::mdc::Print
void Print(int level=0)
Definition: mdc.cc:2707
netpixel::core
bool core
Definition: netpixel.hh:102
ifo
char ifo[NIFO_MAX][8]
Definition: cwb1G_parameters.C:17
watplot::hist2D
TH2F * hist2D
Definition: watplot.hh:175
PlotMultiChirp2
void PlotMultiChirp2(EBBH *ebbh, int ifo=0, TString ptype="", TString wtype="")
Definition: ChirpMass.C:744
MDC_TF
Definition: mdc.hh:174
wavearray::max
virtual DataType_t max() const
Definition: wavearray.cc:1316
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
canvas
TCanvas * canvas
Definition: ChirpMass.C:128
M1
float M1
Definition: cbc_plots_sim4_ifar.C:2342
watplot::canvas
TCanvas * canvas
Definition: watplot.hh:174
hN
TH2F * hN
Definition: ChirpMass.C:120
plot
x plot
Definition: DrawPhaseShift.C:40
gENTRY
int gENTRY
Definition: ChirpMass.C:110
S2
Meyer< double > * S2
Definition: TestTObjectMeyer.C:26
time
double time[6]
Definition: cbc_plots.C:435
G
double G
Definition: DrawEBHH.C:12
DrawMChirpInstantaneousFrequency
void DrawMChirpInstantaneousFrequency(EBBH *ebbh, TString wtype, int ifo)
Definition: ChirpMass.C:1148
CWB::mdc
Definition: mdc.hh:216
gr
TGraph * gr
Definition: cwb_draw_sensitivity.C:4
heT
TH2F * heT
Definition: ChirpMass.C:122
mass
float mass[2]
Definition: cbc_plots.C:437
xgr
TGraphErrors * xgr[NCHIRP_MAX]
Definition: ChirpMass.C:126
WTS
watplot * WTS
Definition: ChirpMass.C:115
SM
static const double SM
Definition: GNGen.cc:8
PCH
watplot * PCH
Definition: ChirpMass.C:116
int
i() int(T_cor *100))
Pi
double Pi
Definition: DrawPhaseShift.C:12
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:4
CWB::STFT::GetHistogram
TH2D * GetHistogram()
Definition: STFT.hh:53
ReadDataFromROOT
int ReadDataFromROOT(TString ifile, EBBH *ebbh)
Definition: ChirpMass.C:324
ChirpIntegralHistogram
double ChirpIntegralHistogram(EBBH *ebbh, TH2F &*h2d, double dMc)
Definition: ChirpMass.C:1090
gOFILE
TString gOFILE
Definition: ChirpMass.C:112
inspOptions
TString inspOptions
Definition: CWB_Plugin_InjectCBC_Config.C:12
iSNR
std::vector< double > iSNR[NIFO_MAX]
Definition: CWB_Plugin_PE.C:408
k
int k
Definition: cwb_report_prod_2.C:150
MDC
CWB::mdc * MDC
Definition: ChirpMass.C:114
maxLayers
#define maxLayers
Definition: ClusterToFrameNew01_dtmin_approx_cent.C:35
PlotHistogram
void PlotHistogram(EBBH *ebbh, TH2F &*hH)
Definition: ChirpMass.C:1038
A
static double A
Definition: geodesics.cc:8
GetCBC
void GetCBC(EBBH *ebbh)
Definition: ChirpMass.C:1283
clusterdata::chi2chirp
float chi2chirp
Definition: netcluster.hh:68
heF
TH2F * heF
Definition: ChirpMass.C:123
F
double F
Definition: DrawPhaseShift.C:14
netpixel::time
size_t time
Definition: netpixel.hh:92
par
vector< mdcpar > par
Definition: CWB_Plugin_HEN_SIM_Config.C:83
entry
double * entry
Definition: cwb_setcuts.C:206
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
FindChirp
void FindChirp(EBBH *ebbh, int mch_order, double &tmerger, double &xmin, double &xmax, int ifo, TString ptype, TString wtype)
Definition: ChirpMass.C:796
froot
TFile * froot
Definition: TestReadWriteDetectorObject.C:11
netcluster::getPixel
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:395
npix
int npix
Definition: ReadHistFromJobFile.C:66
dt
double dt
Definition: DrawPhaseShift.C:13
watconstants::GravitationalConstant
double GravitationalConstant()
Definition: constants.hh:113
RATE
#define RATE
Definition: BaseFunctionWDM.C:6
CHI2_THR
#define CHI2_THR
Definition: ChirpMass.C:99
clusterdata::mchirperr
float mchirperr
Definition: netcluster.hh:67
nifo
int nifo
Definition: cwb_dump_events.C:32
waveform::hx
wavearray< double > hx
Definition: mdc.hh:197
title
char title[256]
Definition: SSeriesExample.C:1
e0
double e0
Definition: RescaleEBBH.C:17
PlotMChirp
void PlotMChirp(EBBH *ebbh)
Definition: ChirpMass.C:685
gps
double gps
Definition: cwb_report_prod_2.C:194
PlotCluster
void PlotCluster(watplot *WTS, netcluster *pwc, double fLow, double fHigh)
Definition: ChirpMass.C:556
CWB::mdc::GetWatPlot
watplot * GetWatPlot()
Definition: mdc.hh:401
T
double T
Definition: testWDM_4.C:11
clusterdata::chirp
TGraphErrors chirp
chirp fit parameters (don't remove ! fix crash when exit from CINT)
Definition: netcluster.hh:75
fLow
double fLow
Definition: cwb2G_parameters.C:19
ebbh
EBBH * ebbh
Definition: LoopChirpMass.C:26
gTYPE
int gTYPE
Definition: ChirpMass.C:109
detector::Name
char Name[16]
Definition: detector.hh:309
fabs
double fabs(const Complex &x)
Definition: numpy.cc:37
fchirp
TF1 * fchirp
Definition: ChirpMass.C:118
CWB::mdc::GetWaveform
waveform GetWaveform(int ID, int id=0)
Definition: mdc.cc:2512
fList
vector< TString > fList
Definition: LoopChirpMass.C:30
GetDifWaveform
wavearray< double > GetDifWaveform(wavearray< double > *wf1, wavearray< double > *wf2)
Definition: ChirpMass.C:514
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
CWB::mdc::SetZoom
void SetZoom(double epzoom=EPZOOM)
Definition: mdc.hh:396
itree
TTree * itree
Definition: cwb_dump_events.C:27
wavearray::data
DataType_t * data
Definition: wavearray.hh:301
M2
float M2
Definition: cbc_plots_sim4_ifar.C:2342
K
int K
Definition: CWB_Plugin_Sim4_BRST_Config.C:87
itime
double * itime
Definition: cwb_condor_create_ced.C:53
dFreq
double dFreq
Definition: DrawPSD.C:17
CWB::Toolbox::getHilbertEnvelope
static wavearray< double > getHilbertEnvelope(wavearray< double > x)
Definition: Toolbox.cc:6131
dT
double dT
Definition: testWDM_5.C:12
wINJ
wavearray< double > wINJ[NIFO_MAX]
Definition: CWB_Plugin_xWRC.C:73
netpixel::rate
float rate
Definition: netpixel.hh:95
netpixel::getdata
double getdata(char type='R', size_t n=0)
Definition: netpixel.hh:56
TMERGER_CUT
#define TMERGER_CUT
Definition: ChirpMass.C:100
GetFitParameter
double GetFitParameter(double mchirp)
Definition: ChirpMass.C:540
fName
char fName[256]
Definition: CreateListEBBH.C:25
GetAlignedWaveform
wavearray< double > GetAlignedWaveform(wavearray< double > *wf, wavearray< double > *wref)
Definition: ChirpMass.C:487
rp0
double rp0
Definition: RescaleEBBH.C:16
GetEBBH
void GetEBBH(EBBH *ebbh)
Definition: ChirpMass.C:1340
GetTimeBoundaries
double GetTimeBoundaries(wavearray< double > x, double P, double &bT, double &eT)
Definition: ChirpMass.C:455
MDC_EBBH
Definition: mdc.hh:129
stft
CWB::STFT * stft
Definition: ChirpMass.C:117
exit
exit(0)
watconstants::SpeedOfLightInVacuo
double SpeedOfLightInVacuo()
Definition: constants.hh:96
netpixel::setdata
bool setdata(double a, char type='R', size_t n=0)
Definition: netpixel.hh:40
avr
double avr
Definition: TestNormalization.C:46