Logo Coherent WaveBurst  
Reference Guide
Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cwb_mkprc.C
Go to the documentation of this file.
1 // this macro compute the PRC FOM
2 // it is used in post production
3 // can be run in standalone or it is called from the sim report
4 
5 #include <vector>
6 #include "TVector3.h"
7 #include "TRotation.h"
8 #include "Math/Vector3Dfwd.h"
9 #include "Math/Rotation3D.h"
10 
11 #define RESOLUTION 2
12 #define COORDINATES "Geographic"
13 #define PROJECTION ""
14 #define NMDC_MAX 64
15 
16 void RECvsINJ(int mtype);
17 void PlotSearchArea(int mtype);
18 void PlotEventToAntPat(int mtype, bool binj=true, int polarization=3);
19 void PlotCosOmega(int mtype);
21 void MakeHtmlIndex(TString* ptype, int psize, bool multi, bool header);
23 double GetPercentage(int mtype, int iderA);
24 void ReadInjSet(vector<TString>& mdcType, vector<TString>& mdcSet);
25 
26 // --------------------------------------------------------
27 // Global variables
28 // --------------------------------------------------------
29 TString gODIR = pp_dir; // output dir
30 TString gPTYPE = ""; // plot type
31 TCanvas* gCANVAS = NULL; // canvas object
32 CWB::mdc* gMDC = NULL; // MDC object
33 TTree* gTRWAVE = NULL; // wave tree
34 TTree* gTRMDC = NULL; // mdc tree
35 
36 // -----------------------------------------
37 // plot type : ptype
38 // -----------------------------------------
39 
40 #define nPLOT 6
42  "RECvsINJ", // Injected vs Reconstructed SNRnet distribution
43  "RECvsAP", // Sky Source Distribution of injected MDC
44  "INJvsAP", // Sky Source Distribution of reconstructed MDC
45  "CPROBvsSAREA", // Cumulative Probablitily vs Searched Area
46  "PDENSvsANGLE", // Probability Density vs largest angle between
47  // the source's true location and any point within the searched area
48  "COVvsPERC" // PP-Plot : Coverage vs Percentage
49  };
50 
51 
52 void cwb_mkprc(TString odir="prc", bool multi=true, bool header=false, bool bexit=true) {
53 //
54 // odir : output fad directory (default="prc")
55 // multi : if true (def) the prc plots are produced for each MDC type
56 // otherwise the plots include all MDC types
57 // header : if true (def=false) then the html index includes the cWB header logo
58 // bexit : if true (def) then the macro exit at the end of the execution
59 //
60 
62 
63  TB.checkFile(gSystem->Getenv("CWB_ROOTLOGON_FILE"));
64  TB.checkFile(gSystem->Getenv("CWB_PARAMETERS_FILE"));
65  TB.checkFile(gSystem->Getenv("CWB_UPARAMETERS_FILE"));
66  TB.checkFile(gSystem->Getenv("CWB_PPARAMETERS_FILE"));
67  TB.checkFile(gSystem->Getenv("CWB_UPPARAMETERS_FILE"));
68  TB.checkFile(gSystem->Getenv("CWB_EPPARAMETERS_FILE"));
69 
70  // ------------------------------------------------
71  // create canvas
72  // ------------------------------------------------
73  gCANVAS = new TCanvas("canvas", "canvas",16,30,825,546);
74  gCANVAS->Range(-19.4801,-9.25,-17.4775,83.25);
75  gCANVAS->SetBorderSize(2);
76  gCANVAS->SetFrameFillColor(0);
77  gCANVAS->SetGridx();
78  gCANVAS->SetGridy();
79 
80  gStyle->SetOptFit(kTRUE);
81 
82  // ------------------------------------------------
83  // define output dir
84  // ------------------------------------------------
86 
87  // ------------------------------------------------
88  // export to CINT net,cfg
89  // exec config plugin
90  // ------------------------------------------------
91  char cmd[128];
92  //sprintf(cmd,"network* net = (network*)%p;",net);
93  sprintf(cmd,"network* net = new network;");
94  gROOT->ProcessLine(cmd);
96  cfg->Import();
97  sprintf(cmd,"CWB::config* cfg = (CWB::config*)%p;",cfg);
98  gROOT->ProcessLine(cmd);
99  sprintf(cmd,"int gIFACTOR=1;");
100  gROOT->ProcessLine(cmd);
101  configPlugin.Exec();
102  //MDC.Print();
103  gMDC = &MDC;
104 
105  // ------------------------------------------------
106  // open wave/mdc trees
107  // ------------------------------------------------
108  // open wave file
109  TFile* fwave = TFile::Open(sim_file_name);
110  gTRWAVE = (TTree*) gROOT->FindObject("waveburst");
111  // open mdc file
112  TFile *fmdc = TFile::Open(mdc_file_name);
113  gTRMDC = (TTree*) gROOT->FindObject("mdc");
114 
115  // ------------------------------------------------
116  // create output dir
117  // ------------------------------------------------
118  TB.mkDir(gODIR,!pp_batch);
119 
120  // ------------------------------------------------
121  // create plots
122  // ------------------------------------------------
123  gCANVAS->cd(); // must be done because MakeMultiPlotFAD set a different canvas
124 
125  int nSTART = multi ? 1 : 0;
126  int nSTOP = multi ? gMDC->wfList.size() : 0;
127 
128  if(nSTOP>NMDC_MAX) {
129  cout << endl;
130  cout << "cwb_mkprc.C - Error : max number of MDC types allowed in multi mode is " << NMDC_MAX << endl;
131  cout << " use 'single' option, Ex : cwb_mkprc M1 single header" << endl;
132  cout << endl;
133  exit(1);
134  }
135 
136  for(int n=nSTART;n<=nSTOP;n++) { // loop over mdc types
137  gPTYPE="RECvsINJ"; RECvsINJ(n);
138  gPTYPE="CPROBvsSAREA"; PlotSearchArea(n);
139  gPTYPE="PDENSvsANGLE"; PlotCosOmega(n);
140  gPTYPE="INJvsAP"; PlotEventToAntPat(n, true, 3);
141  gPTYPE="RECvsAP"; PlotEventToAntPat(n, false, 3);
142  gPTYPE="COVvsPERC"; PlotCoverageVsPercentage(n);
143  }
144 
145  // ------------------------------------------------
146  // make prc.html
147  // ------------------------------------------------
148  MakeHtmlIndex(ptype,nPLOT,multi,header);
149 
150  if(bexit) gSystem->Exit(0); else return;
151 }
152 
153 void RECvsINJ(int mtype) {
154 
155  // -------------------------------------------------------------
156  // Plot Injected vs Reconstructed SNRnet distribution
157  // -------------------------------------------------------------
158 
159  char sel[256]="";
160  char cut[256]="";
161  char title[256]="";
162 
163  // INJECTED
164  sprintf(sel,"snr[0]*snr[0]");
165  for(int n=1;n<nIFO;n++) sprintf(sel,"%s+snr[%d]*snr[%d]",sel,n,n);
166  sprintf(sel,"sqrt(%s)",sel);
167  if(mtype==0) {
168  sprintf(cut,"");
169  } else {
170  sprintf(cut,"type[0]==%d",mtype);
171  }
172  //gTRMDC->SetFillColor(kBlue);
173  //gTRMDC->SetFillStyle(3001);
174  gTRMDC->Draw(sel,cut,"");
175  int nmdc = gTRMDC->GetSelectedRows();
176  cout << "nmdc : " << nmdc << endl;
177 
178  TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
179  htemp->GetXaxis()->SetTitle("SNRnet");
180  htemp->GetXaxis()->SetTitleOffset(1.2);
181  htemp->GetYaxis()->SetTitle("counts");
182  htemp->GetYaxis()->SetRangeUser(0.1,pow(10.,TMath::Ceil(TMath::Log10(htemp->GetMaximum()))));
183 
184  // DETECTED
185  sprintf(sel,"range[1]");
186  sprintf(sel,"iSNR[0]");
187  for(int n=1;n<nIFO;n++) sprintf(sel,"%s+iSNR[%d]",sel,n);
188  sprintf(sel,"sqrt(%s)",sel);
189  if(mtype==0) {
190  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
192  } else {
193  sprintf(cut,"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
194  mtype,nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
195  }
196  if(T_vED>0) sprintf(cut,"%s && neted[0]/ecor<%f",cut,T_vED);
197  if(T_pen>0) sprintf(cut,"%s && penalty>%f",cut,T_pen);
198  gTRWAVE->SetFillColor(kRed);
199  //gTRWAVE->SetFillStyle(3001);
200  gTRWAVE->Draw(sel,cut,"same");
201  int nwave = gTRWAVE->GetSelectedRows();
202  cout << "nwave : " << nwave << endl;
203  sprintf(title,"%s",cut);
204 
205  TString mdcName = mtype ? gMDC->wfList[mtype-1].name : "ALL";
206  sprintf(title,"%s - Injected(blue) vs Reconstructed(red) SNRnet distribution : inj = %d : rec : %d",
207  mdcName.Data(),nmdc,nwave);
208  htemp->SetTitle(title);
209 
210  // print plot
211  char ofname[256];
212  if(mtype==0) {
213  sprintf(ofname,"%s/%s_%s.png",
214  gODIR.Data(),gPTYPE.Data(),"ALL");
215  } else {
216  sprintf(ofname,"%s/%s_%s.png",
217  gODIR.Data(),gPTYPE.Data(),gMDC->wfList[mtype-1].name.Data());
218  }
219  cout << ofname << endl;
220 
221  gStyle->SetOptStat(0);
222  gCANVAS->SetLogx(true);
223  gCANVAS->SetLogy(true);
224  gCANVAS->Print(ofname);
225 
226  return;
227 }
228 
229 void
231 
232  // ------------------------------------------------------------
233  // Plot Cumulative Probablitily vs Searched Area
234  // ------------------------------------------------------------
235 
236  using namespace ROOT::Math;
237 
238  char sel[256];
239  char cut[256];
240 
241  // DETECTED
242  sprintf(sel,"erA[0]");
243  if(mtype==0) {
244  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
246  } else {
247  sprintf(cut,"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
249  }
250  if(T_vED>0) sprintf(cut,"%s && neted[0]/ecor<%f",cut,T_vED);
251  if(T_pen>0) sprintf(cut,"%s && penalty>%f",cut,T_pen);
252  gTRWAVE->SetFillColor(kRed);
253  gTRWAVE->Draw(sel,cut,"goff");
254  int nwave = gTRWAVE->GetSelectedRows();
255  cout << "nwave : " << nwave << endl;
256 
257  double* erA0 = gTRWAVE->GetV1();
258  Int_t *index = new Int_t[nwave];
259 
260  TMath::Sort(nwave,gTRWAVE->GetV1(),index,false);
261 
262  TH1F* hist = new TH1F("hist","hist",100000,0,180*180);
263 
264  for(int i=0;i<nwave;i++) {
265  double sarea = erA0[index[i]]*erA0[index[i]];
266  hist->Fill(sarea);
267  }
268 
269  // Normalization
270  double integral = 0;
271  for (int i=0;i<=hist->GetNbinsX();i++) integral+=hist->GetBinContent(i);
272  for (int i=0;i<=hist->GetNbinsX();i++) hist->SetBinContent(i,hist->GetBinContent(i)/integral);
273  for (int i=2;i<=hist->GetNbinsX();i++) hist->SetBinContent(i,hist->GetBinContent(i)+hist->GetBinContent(i-1));
274  // find area at 50%
275  double search_area50=0;
276  for (int i=0;i<=hist->GetNbinsX();i++) if(hist->GetBinContent(i)>0.5) {search_area50=hist->GetBinCenter(i);break;}
277  cout << "search_area50 : " << search_area50 << endl;
278 
279  char title[256];
280  TString mdcName = mtype ? gMDC->wfList[mtype-1].name : "ALL";
281  sprintf(title,"%s - Searched Area : Ndet = %d",mdcName.Data(),nwave);
282  hist->SetTitle(title);
283  hist->GetXaxis()->SetTitle("searched area (deg^{2})");
284  hist->GetXaxis()->SetTitleOffset(1.2);
285  hist->GetYaxis()->SetTitle("cumulative probability");
286  hist->GetXaxis()->SetRangeUser(1,180*180);
287  hist->GetYaxis()->SetRangeUser(1e-2,1);
288  hist->SetLineColor(kRed);
289  hist->SetLineWidth(2);
290 
291  hist->SetStats(kFALSE);
292  hist->Draw();
293 
294  // draw the legend
295  TLegendEntry* entry;
296  TLegend *legend = new TLegend(0.4942529,0.1413502,0.8850575,0.2447257,NULL,"brNDC");
297  legend->SetTextFont(42);
298  legend->SetTextSize(0.03);
299  legend->SetLineColor(kBlack);
300  legend->SetFillColor(kWhite);
301  char label[256];sprintf(label,"searched area 50% : %d (deg^{2})",(int)search_area50);
302  legend->AddEntry(hist,label,"lp");
303  legend->Draw();
304 
305  // print plot
306  char ofname[256];
307  if(mtype==0) {
308  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),"ALL");
309  } else {
310  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),gMDC->wfList[mtype-1].name.Data());
311  }
312  cout << ofname << endl;
313  gCANVAS->SetGridx();
314  gCANVAS->SetGridy();
315  gCANVAS->SetLogx();
316  gCANVAS->SetLogy();
317  gCANVAS->Print(ofname);
318 
319  delete [] index;
320  delete hist;
321 }
322 
323 void
325 
326  // ------------------------------------------------------------------------
327  // Plot Probability Density vs largest angle between
328  // the source's true location and any point within the searched area
329  // ------------------------------------------------------------------------
330 
331  if (!gROOT->GetClass("Polar3DVector")) gSystem->Load("libMathCore");
332 
333  using namespace ROOT::Math;
334 
335  // open reconstructed event file
336 
337  char sel[256];
338  char cut[256];
339 
340  // DETECTED
341  sprintf(sel,"theta[0]:phi[0]:theta[1]:phi[1]");
342  if(mtype==0) {
343  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
345  } else {
346  sprintf(cut,"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
348  }
349  if(T_vED>0) sprintf(cut,"%s && neted[0]/ecor<%f",cut,T_vED);
350  if(T_pen>0) sprintf(cut,"%s && penalty>%f",cut,T_pen);
351  gTRWAVE->SetFillColor(kRed);
352  gTRWAVE->Draw(sel,cut,"goff");
353  int nwave = gTRWAVE->GetSelectedRows();
354  cout << "nwave : " << nwave << endl;
355 
356  double* th0 = gTRWAVE->GetV1();
357  double* ph0 = gTRWAVE->GetV2();
358  double* th1 = gTRWAVE->GetV3();
359  double* ph1 = gTRWAVE->GetV4();
360 
361  TH1F* hist = new TH1F("hist","hist",24,-1,1);
362 
363  double deg2rad = TMath::Pi()/180.;
364 
365  for(int i=0;i<nwave;i++) {
366 
367  Polar3DVector v0(1, th0[i]*deg2rad, ph0[i]*deg2rad);
368  Polar3DVector v1(1, th1[i]*deg2rad, ph1[i]*deg2rad);
369 
370  double Dot = v0.Dot(v1);
371  double dOmega = 180.*TMath::ACos(Dot)/TMath::Pi();
372 
373  hist->Fill(Dot);
374  }
375 
376  // Normalization
377  double integral = 0;
378  for (int i=0;i<=hist->GetNbinsX();i++) integral+=hist->GetBinContent(i);
379  for (int i=0;i<=hist->GetNbinsX();i++) hist->SetBinContent(i,hist->GetBinContent(i)/integral);
380 
381  char title[256];
382  TString mdcName = mtype ? gMDC->wfList[mtype-1].name : "ALL";
383  sprintf(title,"%s - angle offset : #events = %d",mdcName.Data(),nwave);
384  hist->SetTitle(title);
385  hist->GetXaxis()->SetTitle("cos(angle offset)");
386  hist->GetXaxis()->SetTitleOffset(1.2);
387  hist->GetYaxis()->SetTitle("probability density");
388  hist->GetXaxis()->SetRangeUser(-1,1);
389  hist->GetYaxis()->SetRangeUser(0.001,1);
390 
391  hist->SetStats(kFALSE);
392  gCANVAS->SetGridx();
393  gCANVAS->SetGridy();
394  gCANVAS->SetLogy();
395  hist->Draw();
396 
397  // print plot
398  char ofname[256];
399  if(mtype==0) {
400  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),"ALL");
401  } else {
402  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),gMDC->wfList[mtype-1].name.Data());
403  }
404  cout << ofname << endl;
405  gCANVAS->Print(ofname);
406 
407  delete hist;
408 }
409 
410 void
412 
413  // ----------------------------------------------------------
414  // PP-Plot : Coverage vs Percentage
415  // ----------------------------------------------------------
416 
417  TString TITLE = "pp-plot";
418 
419  // if coverage > percentage then erA is too large, the true erA is smaller (erA is conservative)
420  // if coverage < percentage then erA is too small, the true erA is greater (erA is optimistic)
421 
422  double coverage[11];
423  coverage[0]=0.;
424  coverage[10]=100.;
425  for (int j=1;j<10;j++) {
426  coverage[j]=GetPercentage(mtype, j);
427  cout << j*10 << " " << coverage[j] << endl;
428  }
429 
430  gCANVAS->Clear();
431  gCANVAS->ToggleEventStatus();
432  gCANVAS->SetFillColor(0);
433  gCANVAS->SetGridx();
434  gCANVAS->SetGridy();
435  gCANVAS->SetLogy(false);
436  gCANVAS->SetLogx(false);
437 
438  gStyle->SetTitleH(0.032);
439  gStyle->SetTitleW(0.98);
440  gStyle->SetTitleY(0.98);
441  gStyle->SetTitleFont(72);
442  gStyle->SetLineColor(kWhite);
443  gStyle->SetPalette(1,0);
444  gStyle->SetNumberContours(256);
445 
446  double perc[11]={0,10,20,30,40,50,60,70,80,90,100};
447  TGraph* gr[2];
448 
449  gr[0] = new TGraph(11,perc,coverage);
450  gr[0]->SetLineColor(kRed);
451  gr[0]->SetLineWidth(1);
452  gr[0]->SetMarkerSize(1);
453  gr[0]->SetMarkerColor(kRed);
454  gr[0]->SetMarkerStyle(20);
455 
456  double x[2]={0,100};
457  double y[2]={0,100};
458  gr[1] = new TGraph(2,x,y);
459  gr[1]->SetLineColor(1);
460  gr[1]->SetLineWidth(2);
461  gr[1]->SetLineStyle(2);
462 
463  TMultiGraph* mg = new TMultiGraph();
464  mg->Add(gr[0],"lp");
465  mg->Add(gr[1],"lp");
466  mg->Paint("ap");
467  char title[256];
468  TString mdcName = mtype ? gMDC->wfList[mtype-1].name : "ALL";
469  sprintf(title,"%s - %s",mdcName.Data(),TITLE.Data());
470  mg->GetHistogram()->SetTitle(title);
471  mg->GetHistogram()->GetXaxis()->SetTitle("Percentage");
472  mg->GetHistogram()->GetYaxis()->SetTitle("Coverage");
473  mg->GetHistogram()->GetXaxis()->SetRangeUser(0,100);
474  mg->GetHistogram()->GetYaxis()->SetRangeUser(0,100);
475  mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.3);
476  mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);
477  mg->GetHistogram()->GetXaxis()->CenterTitle(true);
478  mg->GetHistogram()->GetYaxis()->CenterTitle(true);
479  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
480  mg->GetHistogram()->GetXaxis()->SetTitleFont(42);
481  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
482  mg->GetHistogram()->GetYaxis()->SetTitleFont(42);
483  mg->GetHistogram()->GetXaxis()->SetNdivisions(511);
484  mg->Draw("a");
485 
486  // print plot
487  char ofname[256];
488  if(mtype==0) {
489  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),"ALL");
490  } else {
491  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),gMDC->wfList[mtype-1].name.Data());
492  }
493  cout << ofname << endl;
494  gCANVAS->Print(ofname);
495 
496  delete mg;
497 }
498 
499 double
500 GetPercentage(int mtype, int iderA) {
501 
502  // open reconstructed event file
503  char sel[256];
504  char cut[256];
505 
506  // DETECTED
507  sprintf(sel,"erA[0]-erA[%d]>>hist_cumulative_erA1(2000,-100,100)",iderA);
508  if(mtype==0) {
509  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
511  } else {
512  sprintf(cut,"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
514  }
515  if(T_vED>0) sprintf(cut,"%s && neted[0]/ecor<%f",cut,T_vED);
516  if(T_pen>0) sprintf(cut,"%s && penalty>%f",cut,T_pen);
517  gTRWAVE->SetFillColor(kRed);
518  gTRWAVE->Draw(sel,cut,"goff");
519  int nwave = gTRWAVE->GetSelectedRows();
520  cout << "nwave : " << nwave << endl;
521  if(nwave==0) return 0;
522 
523  double integral = hist_cumulative_erA1->ComputeIntegral();
524  if (integral==0) {cout << "Empty histogram !!!" << endl;exit(0);}
525  double* cumulative = hist_cumulative_erA1->GetIntegral();
526  for (int i=0;i<hist_cumulative_erA1->GetNbinsX();i++)
527  hist_cumulative_erA1->SetBinContent(i,cumulative[i]/integral);
528  //cout << "integral " << integral << endl;
529 
530  double perc = 100.*hist_cumulative_erA1->GetBinContent(1001);
531  delete hist_cumulative_erA1;
532  return perc;
533 }
534 
535 void
536 PlotEventToAntPat(int mtype, bool binj, int polarization) {
537 
538  // ------------------------------------------------------------
539  // Plot Sky Source Distribution of injected/reconstructed MDC
540  // ------------------------------------------------------------
541 
542  // binj=true -> select injected directions
543  // binj=false -> select reconstructed directions
544  // polarization is the antenna pattern type
545 
546  if (!gROOT->GetClass("Polar3DVector")) gSystem->Load("libMathCore");
547 
548  using namespace ROOT::Math;
549  using namespace CWB;
550 
551  double r2d = 180./TMath::Pi();
552  double d2r = TMath::Pi()/180.;
553 
555  TString dummy=""; // introduced to avoid error message in SetOptions, to be fixed !!!
556 
557  // create network
558  gnetwork* gNET = new gnetwork;
559 
560  // create detector objects
561  detector* pD[NIFO_MAX]; //! pointers to detectors
562  for(int i=0; i<nIFO; i++) {
563  if(strlen(ifo[i])>0) pD[i] = new detector(ifo[i]); // built in detector
564  else pD[i] = new detector(detParms[i]); // user define detector
565  }
566  // add detector object to network
567  for(int i=0; i<nIFO; i++) gNET->add(pD[i]);
568 
569  // setup skymap options
570  gskymap* gSM = gNET->GetGskymap();
572 
573 
574  TH2D* h2 = (TH2D*)gSM->GetHistogram();
575  h2->GetXaxis()->SetTitleSize(0.05);
576  h2->GetXaxis()->SetLabelSize(0.05);
577  h2->GetYaxis()->SetTitleSize(0.05);
578  h2->GetYaxis()->SetLabelSize(0.05);
579  h2->GetYaxis()->SetLabelFont(42);
580  h2->GetYaxis()->SetLabelFont(42);
581  h2->GetXaxis()->SetTitleFont(42);
582  h2->GetYaxis()->SetTitleFont(42);
583  h2->GetZaxis()->SetRangeUser(0,1.0);
584 
585  // draw antenna pattern
586  gNET->DrawAntennaPattern(polarization,0,true);
587 
588  // compute direction of D1 D2 axis -> vector v12
589  XYZVector D1(pD[0]->Rv[0],pD[0]->Rv[1],pD[0]->Rv[2]);
590  XYZVector D2(pD[1]->Rv[0],pD[1]->Rv[1],pD[1]->Rv[2]);
591  D1=D1/speedlight;
592  D2=D2/speedlight;
593  XYZVector D12 = D1-D2;
594  TVector3 vD12(D12.X(),D12.Y(),D12.Z());
595 
596  double th12 = r2d*vD12.Theta();
597  double ph12 = r2d*vD12.Phi();
598  cout << "coordinates D12 " << ph12 << " " << th12 << endl;
599  Polar3DVector v12(1, d2r*th12, d2r*ph12);
600 
601  // open reconstructed event file
602 
603  char sel[256];
604  char cut[256];
605 
606  // DETECTED
607  sprintf(sel,"theta[0]:phi[0]:theta[1]:phi[1]");
608  if(mtype==0) {
609  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
611  } else {
612  sprintf(cut,"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
613  mtype,nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
614  }
615  if(T_vED>0) sprintf(cut,"%s && neted[0]/ecor<%f",cut,T_vED);
616  if(T_pen>0) sprintf(cut,"%s && penalty>%f",cut,T_pen);
617  gTRWAVE->SetFillColor(kRed);
618  gTRWAVE->Draw(sel,cut,"goff");
619  int nwave = gTRWAVE->GetSelectedRows();
620  cout << "nwave : " << nwave << endl;
621 
622  double* th0 = gTRWAVE->GetV1();
623  double* ph0 = gTRWAVE->GetV2();
624  double* th1 = gTRWAVE->GetV3();
625  double* ph1 = gTRWAVE->GetV4();
626 
627  // draw events
628  double ph,th;
629  float markerSize = nwave>10000 ? 0.25 : 0.4;
630  for (int i=0;i<nwave;i++) {
631 
632  if(binj) {ph=ph1[i]; th=th1[i];}
633  else {ph=ph0[i]; th=th0[i];}
634 
635 #if (COORDINATES=="cWB")
636  gSM->DrawMarker(ph, th, 20, markerSize, kBlack); // cWB
637 #endif
638 #if (COORDINATES=="Geographic")
639  double phi,theta;
640  CwbToGeographic(ph,th,phi,theta);
641  gSM->DrawMarker(phi,theta, 20, markerSize, kBlack); // Geographic
642 #endif
643  }
644 
645  // print plot
646  char ofname[256];
647  if(mtype==0) {
648  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),"ALL");
649  } else {
650  sprintf(ofname,"%s/%s_%s.png", gODIR.Data(),gPTYPE.Data(),gMDC->wfList[mtype-1].name.Data());
651  }
652  cout << ofname << endl;
653  gSM->Print(ofname);
654 
655  delete gNET;
656 }
657 
658 void
659 ReadInjSet(vector<TString>& mdcType, vector<TString>& mdcSet) {
660 
661  // read injection file types
662  char imdc_set[NMDC_MAX][128]; // injection set
663  size_t imdc_type[NMDC_MAX]; // injection type
664  char imdc_name[NMDC_MAX][128]; // injection name
665  double imdc_fcentral[NMDC_MAX]; // injection central frequencies
666  double imdc_fbandwidth[NMDC_MAX]; // injection bandwidth frequencies
667  size_t imdc_index[NMDC_MAX]; // type reference array
668  size_t imdc_iset[NMDC_MAX]; // injection set index
669 
670  int ninj=ReadInjType(mdc_inj_file,NMDC_MAX,imdc_set,imdc_type,
671  imdc_name,imdc_fcentral,imdc_fbandwidth);
672  if(ninj==0) {cout << "Error - no injection - terminated" << endl;exit(1);}
673 
675  int nset=0;
676  for(int i=0;i<ninj;i++) {
677  bool bnew=true;
678  for(int j=0;j<nset;j++) if(imdc_set[i]==imdc_set_name[j]) bnew=false;
679  if(bnew) imdc_set_name[nset++]=imdc_set[i];
680  }
681  //cout << "nset : " << nset << endl;
682 
683  for(int i=0;i<nset;i++) {
684  for(int j=0;j<ninj;j++) if(imdc_set[j]==imdc_set_name[i]) imdc_iset[j]=i;
685  }
686  //for (int iset=0;iset<nset;iset++) cout << iset << " " << imdc_set_name[iset].Data() << endl;
687 
688  mdcType.resize(ninj);
689  mdcSet.resize(nset);
690 
691  for(int i=0;i<ninj;i++) mdcType[i]=imdc_set[i];
692  for(int iset=0;iset<nset;iset++) mdcSet[iset]=imdc_set_name[iset];
693 
694  return;
695 }
696 
697 void
698 MakeHtmlIndex(TString* ptype, int psize, bool multi, bool header) {
699 
700  char index_file[1024];
701 
702  // read mdc sets
703  vector<TString> mdcType;
704  vector<TString> mdcSet;
705  ReadInjSet(mdcType, mdcSet);
706 
707  int nset = multi ? mdcSet.size() : 1;
708  int ntype = multi ? mdcType.size() : 1;
709 
710  // make index html file
711 
712  sprintf(index_file,"%s/prc.html",gODIR.Data());
713 
714  ofstream out;
715  out.open(index_file,ios::out); // create index.html
716  char oline[1024];
717 
718  out << "<head>" << endl;
719  out << "<!-- Include the tabber code -->" << endl;
720  out << "<script type=\"text/javascript\" src=\"tabber.js\"></script>" << endl;
721  out << "<link rel=\"stylesheet\" href=\"tabber.css\" TYPE=\"text/css\" MEDIA=\"screen\">" << endl;
722  out << "<script type=\"text/javascript\">" << endl;
723  out << "document.write('<style type=\"text/css\">.tabber{display:none;}<\/style>');" << endl;
724  out << "</script>" << endl;
725  out << "</head>" << endl;
726 
727  out << "<html>" << endl;
728  out << "<br>" << endl;
729  out << "<div class=\"tabber\">" << endl;
730 
731  for(int iset=0;iset<nset;iset++) {
732  TString setName = multi ? mdcSet[iset] : "ALL";
733  int height=1400;
734  if(multi) for(int i=1;i<ntype;i++) if(mdcType[i]==setName) height+=1400;
735  out << "<div class=\"tabbertab\">" << endl;
736  out << "<h2>" << setName << "</h2>" << endl;
737  out << "<iframe src=\""<< setName<<".html\" width=\"100%\" height=\""
738  <<height<<"px\" frameborder=\"0\"></iframe>" << endl;
739  out << "</div>" << endl;
740  }
741  out << "</div>" << endl;
742  out << "</html>" << endl;
743  out.close();
744 
745  // copy tabber javascript
746 
747  char cmd[1024];
748  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.css %s",gODIR.Data());
749  gSystem->Exec(cmd);
750  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.js %s",gODIR.Data());
751  gSystem->Exec(cmd);
752 
753  if(header) AddHtmlHeader(gODIR);
754 
755  // make frame html files
756 
757  for(int iset=0;iset<nset;iset++) {
758 
759  TString setName = multi ? mdcSet[iset] : "ALL";
760 
761  sprintf(index_file,"%s/%s.html",gODIR.Data(),setName.Data());
762 
763  ofstream out;
764  out.open(index_file,ios::out); // create index.html
765  char oline[1024];
766  out << "<html>" << endl;
767 
768  for(int i=0;i<ntype;i++) {
769 
770  if(multi && mdcType[i]!=mdcSet[iset]) continue;
771 
772  char fname[1024];
773  TString mdcName = multi ? gMDC->wfList[i].name : "ALL";
774 
775  // title
776 
777  out << "<p><center>" << endl;
778  out << "<font color=\"red\" style=\"font-weight:bold;\"><h1>"
779  << "<a id="<< mdcName.Data() << ">"
780  << mdcName.Data()<<"</a></h1></font>" << endl;
781  out << "</center>" << endl;
782  out << "<hr><br>" << endl;
783 
784  // injected vs reconstructed source locations
785 
786  out << "<p><table summary=\"\"><tr align=\"left\"><td valign=\"top\" width=\"50%\">"
787  << "<div align=\"center\"><ul><br/>" << endl;
788  out << "<div align=\"center\"><strong>injected source locations</strong></div><br>" << endl;
789  out << "<a class=\"image\" title=\"injected_source_locations\">" << endl;
790  sprintf(fname,"%s_%s.png",ptype[2].Data(),mdcName.Data());
791  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
792  out << oline << endl;
793  out << "</br></ul>" << endl;
794 
795  out << "</div>" << endl;
796  out << "<p></td><td valign=\"top\" width=\"50%\"><div align=\"center\"><ul><br/>" << endl;
797  out << "<div align=\"center\"><strong>reconstructed source locations</strong></div><br>" << endl;
798  out << "<a class=\"image\" title=\"reconstructed_source_locations\">" << endl;
799  sprintf(fname,"%s_%s.png",ptype[1].Data(),mdcName.Data());
800  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
801  out << oline << endl;
802  out << "</br></ul>" << endl;
803 
804  out << "</div>" << endl;
805  out << "<br></td></tr></table>" << endl;
806 
807  // search area & cos(omega) distribution
808 
809  out << "<p><table summary=\"\"><tr align=\"left\"><td valign=\"top\" width=\"50%\">"
810  << "<div align=\"center\"><ul><br/>" << endl;
811  out << "<div align=\"center\"><strong>searched area</strong></div><br>" << endl;
812  out << "<a class=\"image\" title=\"searched_area\">" << endl;
813  sprintf(fname,"%s_%s.png",ptype[3].Data(),mdcName.Data());
814  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
815  out << oline << endl;
816  out << "</br></ul>" << endl;
817 
818  out << "</div>" << endl;
819  out << "<p></td><td valign=\"top\" width=\"50%\"><div align=\"center\"><ul><br/>" << endl;
820  out << "<div align=\"center\"><strong>angle offset</strong></div><br>" << endl;
821  out << "<a class=\"image\" title=\"cos_omega_distribution\">" << endl;
822  sprintf(fname,"%s_%s.png",ptype[4].Data(),mdcName.Data());
823  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
824  out << oline << endl;
825  out << "</br></ul>" << endl;
826 
827  out << "</div>" << endl;
828  out << "<br></td></tr></table>" << endl;
829 
830  // injected vs reconstructed distribution & pp-plot
831 
832  out << "<p><table summary=\"\"><tr align=\"left\"><td valign=\"top\" width=\"50%\">"
833  << "<div align=\"center\"><ul><br/>" << endl;
834  out << "<div align=\"center\"><strong>injected vs reconstructed snr distribution</strong></div>" << endl;
835  out << "<a class=\"image\" title=\"inj_vs_rec_distribution\"><br>" << endl;
836  sprintf(fname,"%s_%s.png",ptype[0].Data(),mdcName.Data());
837  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
838  out << oline << endl;
839  out << "</br></ul>" << endl;
840 
841  out << "</div>" << endl;
842  out << "<p></td><td valign=\"top\" width=\"50%\"><div align=\"center\"><ul><br/>" << endl;
843  out << "<div align=\"center\"><strong>pp-plot</strong></div><br>" << endl;
844  out << "<a class=\"image\" title=\"pp_plot\">" << endl;
845  sprintf(fname,"%s_%s.png",ptype[5].Data(),mdcName.Data());
846  sprintf(oline,"<img src=\"%s\" width=460></a>",fname);
847  out << oline << endl;
848  out << "</br></ul>" << endl;
849 
850  out << "</div>" << endl;
851  out << "<br></td></tr></table>" << endl;
852 
853  }
854 
855  out << "</html>" << endl;
856  out.close();
857  }
858 
859  return;
860 }
861 
863 
864  char cmd[1024];
865 
866  sprintf(cmd,"root -n -l -b ${CWB_ROOTLOGON_FILE} '${HOME_CWB}/info/macros/MakeHTML.C\(\"%s\",false\)'",odir.Data());
867  gSystem->Exec(cmd);
868  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/ROOT.css %s",odir.Data());
869  gSystem->Exec(cmd);
870  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/ROOT.js %s",odir.Data());
871  gSystem->Exec(cmd);
872  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.css %s",odir.Data());
873  gSystem->Exec(cmd);
874  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.js %s",odir.Data());
875  gSystem->Exec(cmd);
876 }
877 
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
TH2D * GetHistogram()
Definition: gskymap.hh:120
gskymap * gSM
char cut[512]
double T_pen
gnetwork * gNET
TString gPTYPE
Definition: cwb_mkprc.C:30
char mdc_inj_file[1024]
size_t add(detector *)
param: detector structure return number of detectors in the network
Definition: network.cc:2528
char cmd[1024]
int nwave
Definition: cbc_plots.C:820
Definition: ced.hh:24
void DrawAntennaPattern(int polarization=-1, int dpaletteId=0, bool btitle=true, int order=6)
Definition: gnetwork.cc:655
int n
Definition: cwb_net.C:10
TString("c")
#define TITLE
Definition: TestSiteJ1.C:44
int nset
Definition: cwb_mkeff.C:57
double T_cor
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
void PlotEventToAntPat(int mtype, bool binj=true, int polarization=3)
Definition: cwb_mkprc.C:536
char odir[1024]
float theta
size_t ntype
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
Definition: gskymap.cc:724
TString gODIR
Definition: cwb_mkprc.C:29
TH2F * ph
CWB::mdc * MDC
i pp_inetcc
void PlotCoverageVsPercentage(int mtype)
Definition: cwb_mkprc.C:411
int polarization
int j
Definition: cwb_net.C:10
i drho i
void Import(TString umacro="")
Definition: config.cc:334
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:3956
size_t imdc_iset[NMDC_MAX]
Definition: cwb_mkeff.C:50
void MakeHtmlIndex(TString *ptype, int psize, bool multi, bool header)
Definition: cwb_mkprc.C:698
CWB::Toolbox TB
Definition: ComputeSNR.C:5
char ifo[NIFO_MAX][8]
#define RESOLUTION
Definition: cwb_mkprc.C:11
#define nIFO
void AddHtmlHeader(TString odir)
Definition: cwb_mkprc.C:862
ofstream out
Definition: cwb_merge.C:196
vector< int > mtype
Definition: cwb_report_pe.C:48
int nmdc
Definition: cbc_plots.C:791
float phi
TString ptype[nPLOT]
Definition: cwb_mkprc.C:41
size_t imdc_index[NMDC_MAX]
Definition: cwb_mkeff.C:49
void RECvsINJ(int mtype)
Definition: cwb_mkprc.C:153
Definition: mdc.hh:216
#define COORDINATES
Definition: cwb_mkprc.C:12
TGraph * gr
TTree * gTRMDC
Definition: cwb_mkprc.C:34
detector D1
double deg2rad
void PlotCosOmega(int mtype)
Definition: cwb_mkprc.C:324
TString label
Definition: MergeTrees.C:21
double Pi
void ReadInjSet(vector< TString > &mdcType, vector< TString > &mdcSet)
Definition: cwb_mkprc.C:659
char oline[1024]
const int NIFO_MAX
Definition: wat.hh:4
TTree * gTRWAVE
Definition: cwb_mkprc.C:33
int ninj
Definition: cwb_mkeff.C:52
char fname[1024]
#define speedlight
Definition: watfun.hh:15
#define NMDC_MAX
Definition: cwb_mkprc.C:14
TCanvas * gCANVAS
Definition: cwb_mkprc.C:31
#define nPLOT
Definition: cwb_mkprc.C:40
char imdc_name[NMDC_MAX][128]
Definition: cwb_mkeff.C:46
double * entry
Definition: cwb_setcuts.C:206
double e
vector< waveform > wfList
Definition: mdc.hh:355
int ReadInjType(TString ifName, int ntype_max, char set[][128], size_t type[], char name[][128], double fcentral[], double fbandwidth[])
Definition: Toolfun.hh:740
char sim_file_name[1024]
gskymap * GetGskymap()
Definition: gnetwork.hh:26
void PlotSearchArea(int mtype)
Definition: cwb_mkprc.C:230
char title[256]
Definition: SSeriesExample.C:1
wavearray< int > index
double T_win
char pp_dir[512]
Definition: test_config1.C:155
TString sel("slag[1]:rho[1]")
static void mkDir(TString dir, bool question=false, bool remove=true)
Definition: Toolbox.cc:4000
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:396
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double T_cut
double GetPercentage(int mtype, int iderA)
Definition: cwb_mkprc.C:500
TMacro configPlugin
void cwb_mkprc(TString odir="prc", bool multi=true, bool header=false, bool bexit=true)
Definition: cwb_mkprc.C:52
size_t imdc_type[NMDC_MAX]
Definition: cwb_mkeff.C:45
double T_vED
char mdc_file_name[1024]
#define PROJECTION
Definition: cwb_mkprc.C:13
TString config
detectorParams detParms[4]
TString * imdc_set_name
Definition: cwb_mkeff.C:56
void Print(TString pname)
Definition: gskymap.cc:1104
wavearray< double > y
Definition: Test10.C:31
double imdc_fcentral[NMDC_MAX]
Definition: cwb_mkeff.C:47
TMultiGraph * mg
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:66
i drho pp_irho
detector ** pD
char imdc_set[NMDC_MAX][128]
Definition: cwb_mkeff.C:44
double imdc_fbandwidth[NMDC_MAX]
Definition: cwb_mkeff.C:48
exit(0)
double SpeedOfLightInVacuo()
Definition: constants.hh:96
CWB::mdc * gMDC
Definition: cwb_mkprc.C:32