8 #include "Math/Vector3Dfwd.h"
9 #include "Math/Rotation3D.h"
12 #define COORDINATES "Geographic"
24 void ReadInjSet(vector<TString>& mdcType, vector<TString>& mdcSet);
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"));
73 gCANVAS =
new TCanvas(
"canvas",
"canvas",16,30,825,546);
74 gCANVAS->Range(-19.4801,-9.25,-17.4775,83.25);
80 gStyle->SetOptFit(kTRUE);
93 sprintf(cmd,
"network* net = new network;");
94 gROOT->ProcessLine(cmd);
97 sprintf(cmd,
"CWB::config* cfg = (CWB::config*)%p;",cfg);
98 gROOT->ProcessLine(cmd);
100 gROOT->ProcessLine(cmd);
110 gTRWAVE = (TTree*) gROOT->FindObject(
"waveburst");
113 gTRMDC = (TTree*) gROOT->FindObject(
"mdc");
125 int nSTART = multi ? 1 : 0;
126 int nSTOP = multi ? gMDC->
wfList.size() : 0;
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;
136 for(
int n=nSTART;
n<=nSTOP;
n++) {
150 if(bexit) gSystem->Exit(0);
else return;
170 sprintf(cut,
"type[0]==%d",mtype);
176 cout <<
"nmdc : " << nmdc << endl;
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()))));
190 sprintf(cut,
"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
193 sprintf(cut,
"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
202 cout <<
"nwave : " << nwave << endl;
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",
208 htemp->SetTitle(title);
219 cout << ofname << endl;
221 gStyle->SetOptStat(0);
244 sprintf(cut,
"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
247 sprintf(cut,
"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
255 cout <<
"nwave : " << nwave << endl;
257 double* erA0 =
gTRWAVE->GetV1();
262 TH1F*
hist =
new TH1F(
"hist",
"hist",100000,0,180*180);
265 double sarea = erA0[index[
i]]*erA0[index[
i]];
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));
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;
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(1
e-2,1);
288 hist->SetLineColor(kRed);
289 hist->SetLineWidth(2);
291 hist->SetStats(kFALSE);
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");
312 cout << ofname << endl;
331 if (!gROOT->GetClass(
"Polar3DVector")) gSystem->Load(
"libMathCore");
341 sprintf(sel,
"theta[0]:phi[0]:theta[1]:phi[1]");
343 sprintf(cut,
"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
346 sprintf(cut,
"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
354 cout <<
"nwave : " << nwave << endl;
356 double* th0 =
gTRWAVE->GetV1();
357 double* ph0 =
gTRWAVE->GetV2();
358 double* th1 =
gTRWAVE->GetV3();
359 double* ph1 =
gTRWAVE->GetV4();
361 TH1F*
hist =
new TH1F(
"hist",
"hist",24,-1,1);
367 Polar3DVector v0(1, th0[
i]*deg2rad, ph0[
i]*deg2rad);
368 Polar3DVector v1(1, th1[
i]*deg2rad, ph1[
i]*deg2rad);
370 double Dot = v0.Dot(v1);
371 double dOmega = 180.*TMath::ACos(Dot)/
TMath::Pi();
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);
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);
391 hist->SetStats(kFALSE);
404 cout << ofname << endl;
425 for (
int j=1;
j<10;
j++) {
427 cout <<
j*10 <<
" " << coverage[
j] << endl;
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);
446 double perc[11]={0,10,20,30,40,50,60,70,80,90,100};
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);
458 gr[1] =
new TGraph(2,x,y);
459 gr[1]->SetLineColor(1);
460 gr[1]->SetLineWidth(2);
461 gr[1]->SetLineStyle(2);
463 TMultiGraph*
mg =
new TMultiGraph();
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);
493 cout << ofname << endl;
507 sprintf(sel,
"erA[0]-erA[%d]>>hist_cumulative_erA1(2000,-100,100)",iderA);
509 sprintf(cut,
"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
512 sprintf(cut,
"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
520 cout <<
"nwave : " << nwave << endl;
521 if(nwave==0)
return 0;
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);
530 double perc = 100.*hist_cumulative_erA1->GetBinContent(1001);
531 delete hist_cumulative_erA1;
546 if (!gROOT->GetClass(
"Polar3DVector")) gSystem->Load(
"libMathCore");
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);
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]);
593 XYZVector D12 = D1-D2;
594 TVector3 vD12(D12.X(),D12.Y(),D12.Z());
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);
607 sprintf(sel,
"theta[0]:phi[0]:theta[1]:phi[1]");
609 sprintf(cut,
"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
612 sprintf(cut,
"type[1]==%d && abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
620 cout <<
"nwave : " << nwave << endl;
622 double* th0 =
gTRWAVE->GetV1();
623 double* ph0 =
gTRWAVE->GetV2();
624 double* th1 =
gTRWAVE->GetV3();
625 double* ph1 =
gTRWAVE->GetV4();
629 float markerSize = nwave>10000 ? 0.25 : 0.4;
630 for (
int i=0;i<
nwave;i++) {
632 if(binj) {ph=ph1[
i]; th=th1[
i];}
633 else {ph=ph0[
i]; th=th0[
i];}
635 #if (COORDINATES=="cWB")
636 gSM->
DrawMarker(ph, th, 20, markerSize, kBlack);
638 #if (COORDINATES=="Geographic")
641 gSM->
DrawMarker(phi,theta, 20, markerSize, kBlack);
652 cout << ofname << endl;
659 ReadInjSet(vector<TString>& mdcType, vector<TString>& mdcSet) {
671 imdc_name,imdc_fcentral,imdc_fbandwidth);
672 if(ninj==0) {cout <<
"Error - no injection - terminated" << endl;
exit(1);}
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];
684 for(
int j=0;
j<
ninj;
j++)
if(imdc_set[
j]==imdc_set_name[
i]) imdc_iset[
j]=
i;
688 mdcType.resize(ninj);
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];
700 char index_file[1024];
703 vector<TString> mdcType;
704 vector<TString> mdcSet;
707 int nset = multi ? mdcSet.size() : 1;
708 int ntype = multi ? mdcType.size() : 1;
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;
727 out <<
"<html>" << endl;
728 out <<
"<br>" << endl;
729 out <<
"<div class=\"tabber\">" << endl;
731 for(
int iset=0;iset<
nset;iset++) {
732 TString setName = multi ? mdcSet[iset] :
"ALL";
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;
741 out <<
"</div>" << endl;
742 out <<
"</html>" << endl;
748 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.css %s",
gODIR.Data());
750 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.js %s",
gODIR.Data());
757 for(
int iset=0;iset<
nset;iset++) {
759 TString setName = multi ? mdcSet[iset] :
"ALL";
761 sprintf(index_file,
"%s/%s.html",
gODIR.Data(),setName.Data());
766 out <<
"<html>" << endl;
770 if(multi && mdcType[
i]!=mdcSet[iset])
continue;
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;
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;
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;
804 out <<
"</div>" << endl;
805 out <<
"<br></td></tr></table>" << endl;
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;
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;
827 out <<
"</div>" << endl;
828 out <<
"<br></td></tr></table>" << endl;
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;
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;
850 out <<
"</div>" << endl;
851 out <<
"<br></td></tr></table>" << endl;
855 out <<
"</html>" << endl;
866 sprintf(cmd,
"root -n -l -b ${CWB_ROOTLOGON_FILE} '${HOME_CWB}/info/macros/MakeHTML.C\(\"%s\",false\)'",odir.Data());
868 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/ROOT.css %s",odir.Data());
870 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/ROOT.js %s",odir.Data());
872 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.css %s",odir.Data());
874 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.js %s",odir.Data());
size_t add(detector *)
param: detector structure return number of detectors in the network
void DrawAntennaPattern(int polarization=-1, int dpaletteId=0, bool btitle=true, int order=6)
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)
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
void PlotCoverageVsPercentage(int mtype)
void Import(TString umacro="")
size_t imdc_iset[NMDC_MAX]
void MakeHtmlIndex(TString *ptype, int psize, bool multi, bool header)
void AddHtmlHeader(TString odir)
size_t imdc_index[NMDC_MAX]
void PlotCosOmega(int mtype)
void ReadInjSet(vector< TString > &mdcType, vector< TString > &mdcSet)
char imdc_name[NMDC_MAX][128]
vector< waveform > wfList
void PlotSearchArea(int mtype)
TString sel("slag[1]:rho[1]")
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double GetPercentage(int mtype, int iderA)
void cwb_mkprc(TString odir="prc", bool multi=true, bool header=false, bool bexit=true)
size_t imdc_type[NMDC_MAX]
detectorParams detParms[4]
void Print(TString pname)
double imdc_fcentral[NMDC_MAX]
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
char imdc_set[NMDC_MAX][128]
double imdc_fbandwidth[NMDC_MAX]
double SpeedOfLightInVacuo()