19 TString pdir,
double P,
bool freq=
false,
bool showerr=
true);
34 #define IN_CWB_ASCII_FILE_L1 "L1_pe_wave.dat"
35 #define IN_CWB_ASCII_FILE_H1 "H1_pe_wave.dat"
37 #define OUT_CWB_ASCII_FILE_L1 "OUT_L1_pe_wave.dat"
38 #define OUT_CWB_ASCII_FILE_H1 "OUT_H1_pe_wave.dat"
76 bool exit =
label.Contains(
"EXIT") ?
true :
false;
77 label.ReplaceAll(
"EXIT",
"");
85 if(
label==
"GW150914") {tshift=-0.0069;tchirp=8.42;mchirp=30.0;}
86 if(
label==
"GW151226") {tshift=-0.0011;tchirp=7.65;mchirp=8.9;}
87 if(
label==
"LVT151012") {tshift= 0.0010;tchirp=7.43;mchirp=24.0;}
107 std::vector<TString> ifoname;
112 &frec[ifo], &fmed[ifo], &fl50[ifo], &fu50[ifo], &fl90[ifo], &fu90[ifo]);
114 &fu90[ifo], &wrec[ifo],
".",
P,
true);
119 ifoname =
ReadDataFromROOT(ipath,
n, &winj[
n], &wrec[n], &wwht[n], &wmed[n], &wl50[n], &wu50[n], &wl90[n], &wu90[n]);
120 nifo = ifoname.size();
123 if(
ifo>nifo-1) {cout <<
"DrawWaveformPE Error : input ifo=" <<
ifo <<
" parameter must be < nifo=" << nifo << endl;
exit(1);}
131 double tmed = wrec[0].
start()+(wrec[0].
size()/wrec[0].
rate())/2.;
133 double tstop = tmed+1;
139 PlotWaveformAsymmErrors(
"",
"", &wrec[
ifo], &wwht[ifo], &wl50[ifo], &wu50[ifo], &wl90[ifo], &wu90[ifo], &wrec[ifo],
"",
P,
false,
false);
144 PlotWaveformAsymmErrors(
"",
"", &wrec[
ifo], &wmed[ifo], &wl50[ifo], &wu50[ifo], &wl90[ifo], &wu90[ifo], &wrec[ifo],
"",
P);
149 PlotWaveformAsymmErrors(
"",
"", &winj[
ifo], &wmed[ifo], &wl50[ifo], &wu50[ifo], &wl90[ifo], &wu90[ifo], &wrec[ifo],
"",
P);
156 PlotWaveformAsymmErrors(
"",
"", &wdif[ifo], &wdif[ifo], &wl50[ifo], &wu50[ifo], &wl90[ifo], &wu90[ifo], &wrec[ifo],
"",
P);
162 TString ofile = TString::Format(
"%s_whitened_data_%s_Spectrogram.png",
label.Data(),ifoname[
ifo].Data());
164 if(exit)
exit(0);
else return;
168 TString title = TString::Format(
"%s whitened data %s + %s(inv/tshift)",
label.Data(),ifoname[0].Data(),ifoname[1].Data());
169 TString ofile = TString::Format(
"%s_whitened_data_%s_plus_%s_Spectrogram.png",
label.Data(),ifoname[0].Data(),ifoname[1].Data());
177 if(exit)
exit(0);
else return;
181 TString title = TString::Format(
"%s whitened data %s - %s(inv/tshift)",
label.Data(),ifoname[0].Data(),ifoname[1].Data());
182 TString ofile = TString::Format(
"%s_whitened_data_%s_minus_%s_Spectrogram.png",
label.Data(),ifoname[0].Data(),ifoname[1].Data());
190 if(exit)
exit(0);
else return;
194 DumpWaveform(
ifo, &wrec[
ifo], &wmed[ifo], &wl50[ifo], &wu50[ifo], &wl90[ifo], &wu90[ifo]);
229 TFile*
froot =
new TFile(ifile,
"READ");
230 if(froot==NULL) {cout <<
"ReadWaveformPE Error : opening input root file" << endl;
exit(1);}
231 TTree*
gTREE = froot->Get(
"waveburst");
232 if(gTREE==NULL) {cout <<
"ReadWaveformPE Error : no waveburst present in the file" << endl;
exit(1);}
235 TList*
list = gTREE->GetUserInfo();
236 int nIFO=list->GetSize();
237 std::vector<TString> ifoname;
238 if (nIFO==0) {cout <<
"ReadWaveformPE Error : no ifo present in the tree" << endl;
exit(1);}
239 for (
int n=0;
n<list->GetSize();
n++) {
242 ifoname.push_back(pDetector->
Name);
260 gTREE->SetBranchAddress(
"pe_trials",&pe_trials);
261 gTREE->SetBranchAddress(
"pe_erR",pe_erR);
262 gTREE->SetBranchAddress(
"pe_erF",pe_erF);
263 gTREE->SetBranchAddress(
"pe_pnul",pe_pnul);
264 gTREE->SetBranchAddress(
"pe_snet",pe_snet);
265 gTREE->SetBranchAddress(
"pe_ff",pe_ff);
266 gTREE->SetBranchAddress(
"pe_of",pe_of);
267 gTREE->SetBranchAddress(
"pe_mch",pe_mch);
270 gTREE->SetBranchAddress(TString::Format(
"pe_wINJ_%d",
n).Data(),&pe_wINJ[
n]);
271 gTREE->SetBranchAddress(TString::Format(
"pe_wREC_%d",n).Data(),&pe_wREC[n]);
272 gTREE->SetBranchAddress(TString::Format(
"pe_wWHT_%d",n).Data(),&pe_wWHT[n]);
273 gTREE->SetBranchAddress(TString::Format(
"pe_wMED_%d",n).Data(),&pe_wMED[n]);
274 gTREE->SetBranchAddress(TString::Format(
"pe_wL50_%d",n).Data(),&pe_wL50[n]);
275 gTREE->SetBranchAddress(TString::Format(
"pe_wU50_%d",n).Data(),&pe_wU50[n]);
276 gTREE->SetBranchAddress(TString::Format(
"pe_wL90_%d",n).Data(),&pe_wL90[n]);
277 gTREE->SetBranchAddress(TString::Format(
"pe_wU90_%d",n).Data(),&pe_wU90[n]);
278 gTREE->SetBranchAddress(TString::Format(
"pe_wAVR_%d",n).Data(),&pe_wAVR[n]);
279 gTREE->SetBranchAddress(TString::Format(
"pe_wRMS_%d",n).Data(),&pe_wRMS[n]);
285 cout <<
"---------------------> size " << pe_wMED[0]->
size() << endl;
286 cout <<
"---------------------> pe_trials " << pe_trials << endl;
299 for(
int n=0;n<
nIFO;n++) {
325 for(
int i=0;
i<
N;
i++) {avr+=
i*x[
i]*x[
i]; E+=x[
i]*x[
i];}
332 double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
333 for(
int j=1;
j<
N;
j++) {
334 a = ((M-
j>=0)&&(M-j<N)) ? x[M-j] : 0.;
335 b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
351 TString pdir,
double P,
bool freq,
bool showerr) {
367 for(
int i=0;
i<wrec->
size();
i++) {
368 if(time[
i]>bT && time[
i]<eT)
continue;
373 TGraphAsymmErrors* egr90 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl90->
data,wu90->
data);
374 egr90->SetLineColor(17);
375 egr90->SetFillStyle(1001);
376 egr90->SetFillColor(17);
377 egr90->GetXaxis()->SetTitle(
"time (s)");
378 if(freq) egr90->GetYaxis()->SetTitle(
"frequency (hz)");
379 else egr90->GetYaxis()->SetTitle(
"magnitude");
380 egr90->SetTitle(title);
382 TGraphAsymmErrors* egr50 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl50->
data,wu50->
data);
383 egr50->SetLineColor(15);
384 egr50->SetFillStyle(1001);
385 egr50->SetFillColor(15);
387 TGraph* agr =
new TGraph(size,time.
data,wmed->
data);
388 agr->SetLineWidth(1);
389 agr->SetLineColor(kWhite);
390 agr->SetLineStyle(1);
392 TGraph*
gr =
new TGraph(size,time.
data,wrec->
data);
396 TCanvas*
canvas =
new TCanvas(
"wrec",
"wrec",200,20,800,500);
401 egr90->GetXaxis()->SetRangeUser(bT, eT);
404 egr50->Draw(
"4same");
408 agr->GetXaxis()->SetTitle(
"time (s)");
409 if(freq) agr->GetYaxis()->SetTitle(
"frequency (hz)");
410 else agr->GetYaxis()->SetTitle(
"magnitude");
411 agr->SetTitle(title);
412 agr->GetXaxis()->SetRangeUser(bT, eT);
413 agr->SetLineColor(16);
421 canvas->Print(ofname);
422 cout <<
"write : " << ofname << endl;
440 if(wf1==NULL)
return wf;
441 if(wf1->
size()==0)
return wf;
443 double R=wf1->
rate();
445 double b_wf1 = wf1->
start();
446 double e_wf1 = wf1->
start()+wf1->
size()/R;
447 double b_wf2 = wf2->
start();
448 double e_wf2 = wf2->
start()+wf2->
size()/R;
450 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
451 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
453 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
454 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
455 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
457 for(
int i=0;
i<sizeXCOR;
i++) wf[
i+o_wf2] = wf1->
data[
i+o_wf1];
464 double R=wf1->
rate();
466 double b_wf1 = wf1->
start();
467 double e_wf1 = wf1->
start()+wf1->
size()/R;
468 double b_wf2 = wf2->
start();
469 double e_wf2 = wf2->
start()+wf2->
size()/R;
471 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
472 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
474 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
475 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
476 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
481 wfdif.
start(b_wf1+
double(o_wf1)/R);
483 for(
int i=0;
i<sizeXCOR;
i++) wfdif[
i] = wf1->
data[
i+o_wf1] - wf2->
data[
i+o_wf2];
506 for(
int i=0;
i<wrec->
size();
i++) {
507 if(time[
i]>bT && time[
i]<eT)
continue;
512 TGraphAsymmErrors* egr90 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl90->
data,wu90->
data);
513 egr90->SetLineColor(17);
514 egr90->SetFillStyle(1001);
515 egr90->SetFillColor(17);
516 egr90->GetXaxis()->SetTitle(
"time (s)");
517 if(freq) egr90->GetYaxis()->SetTitle(
"frequency (hz)");
518 else egr90->GetYaxis()->SetTitle(
"magnitude");
519 egr90->SetTitle(title);
521 TGraphAsymmErrors* egr50 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl50->
data,wu50->
data);
522 egr50->SetLineColor(15);
523 egr50->SetFillStyle(1001);
524 egr50->SetFillColor(15);
526 TGraph* agr =
new TGraph(size,time.
data,wmed->
data);
527 agr->SetLineWidth(1);
528 agr->SetLineColor(kWhite);
529 agr->SetLineStyle(1);
531 TGraph*
gr =
new TGraph(size,time.
data,wrec->
data);
535 TCanvas*
canvas =
new TCanvas(
"wrec",
"wrec",200,20,800,500);
540 egr90->GetXaxis()->SetRangeUser(bT, eT);
542 egr50->Draw(
"4same");
548 canvas->Print(ofname);
549 cout <<
"write : " << ofname << endl;
571 if (!out.good()) {cout <<
"Error Opening Output File : " << ofname << endl;
exit(1);}
572 cout <<
"Create Output File : " << ofname << endl;
576 out <<
"#whitened data : time, amp_point, amp_median, amp_lower_50_perc, amp_lower_90_perc, amp_upper_50_perc, amp_upper_90_perc" << endl;
580 double dt=1./wrec->
rate();
584 double vrec = wrec->
data[
i];
585 double vmed = wmed->
data[
i];
591 double l50 = vmed-vl50;
592 double u50 = vmed+vu50;
593 double l90 = vmed-vl90;
594 double u90 = vmed+vu90;
597 out << time <<
" " << vrec <<
" " << vmed <<
" " << l50 <<
" " << l90 <<
" " << u50 <<
" " << u90 << endl;
624 if(!in.good()) {cout <<
"Error Opening File : " << ifname << endl;
exit(1);}
629 in.getline(str,1024);
631 in.getline(str,1024);
632 if (!in.good())
break;
635 cout <<
"size " << size << endl;
636 in.clear(ios::goodbit);
637 in.seekg(0, ios::beg);
638 if (size==0) {cout <<
"Error : File " <<
ifile<<
" is empty" << endl;gSystem->Exit(1);}
641 double vrec, vavr, vrms, vmed, vl50, vl90, vu50, vu90;
642 double grec, gavr, grms, gmed, gl50, gl90, gu50, gu90;
673 in.getline(str,1024);
679 >> vrec >> vavr >> vrms >> vmed >> vl50 >> vl90 >> vu50 >> vu90
680 >> grec >> gavr >> grms >> gmed >> gl50 >> gl90 >> gu50 >> gu90;
682 if (!in.good())
break;
684 wrec->
data[
i] = vrec;
685 wmed->
data[
i] = vmed;
688 wl50->
data[
i] = vmed-vl50;
689 wu50->
data[
i] = vu50-vmed;
690 wl90->
data[
i] = vmed-vl90;
691 wu90->
data[
i] = vu90-vmed;
703 fmed->
data[
i] = gmed;
706 fl50->
data[
i] = gmed-gl50;
707 fu50->
data[
i] = gu50-gmed;
708 fl90->
data[
i] = gmed-gl90;
709 fu90->
data[
i] = gu90-gmed;
743 double dF=(x->
rate()/(double)x->
size())/2.;
745 for(
int j=0;
j<x->
size()/2;
j+=2) {
747 if(F<bF || F>eF) {x->
data[
j]=0;x->
data[
j+1]=0;}
754 if(tstart==0) tstart=x->
start();
755 if(tstop==0) tstop=x->
stop();
767 double p0 = 256.*Pi/5*pow(G*Mc*SM*Pi/C/C/C, 5./3);
774 fchirp =
new TF1(
"fchirp",
"pow([0]*([1]-x),-3./8.)", tmin, tmax);
775 fchirp->SetLineColor(kWhite);
778 fchirp->SetParameter(0, p0);
779 fchirp->SetParameter(1, p1);
788 double fparm=nfact*6;
793 int ysize=ystop-ystart;
800 stft =
new CWB::STFT(y,nfft,noverlap,
"energy",
"gauss",fparm);
803 double tstart = nfft/x->
rate()+ystart/x->
rate();
806 tstart+=0.9;tstop-=0.9;
807 stft->
Draw(tstart,tstop,16,1024,0,0,1);
812 canvas->SetLogy(
true);
814 if(ofname!=
"") stft->
Print(ofname);
detectorParams getDetectorParams()
virtual size_t size() const
virtual void rate(double r)
void Print(TString pname)
wavearray< double > a(hp.size())
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
watplot p1(const_cast< char * >("TFMap1"))
virtual void start(double s)
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
double GravitationalConstant()
static void TimeShift(wavearray< double > &x, double tShift=0.)
virtual void stop(double s)
double fabs(const Complex &x)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
virtual void resize(unsigned int)
gwavearray< double > * grec
double SpeedOfLightInVacuo()