27 #pragma GCC system_header 34 #include "TObjArray.h" 35 #include "TObjString.h" 36 #include "TPaletteAxis.h" 44 #include "TTreeFormula.h" 53 #define maxLayers 1025 55 void PlotFrame(std::vector<double>* nn_frame,
int cid);
57 double ConvertToFrame(
netcluster*
pwc,
int cid,
int nifo,
char type,
int irate,
double &
dT,
double &dF,
double &Fc,
int &ninf,std::vector<double>* nn_frame, std::vector<double>*
dt);
69 TFile*
ifile =
new TFile(ifName);
71 cout <<
"ClusterToFrame - Error : error opening file : " << ifName.Data() << endl;
76 TTree* nn_itree = (TTree *) ifile->Get(
"waveburst");
78 cout <<
"ClusterToFrame - Error : tree waveburst not found !!!" << endl;
83 nn_itree->SetBranchAddress(
"netcluster",&pwc);
84 int Ntot=nn_itree->GetEntries();
87 int nIFO = nn_itree->GetUserInfo()->GetSize();
89 cout<<
"ClusterToFrame - Error : wrong user detector list in header tree"<<endl;
94 TFile
ofile(ofName,
"RECREATE");
96 TTree* nn_otree = (TTree*)nn_itree->CloneTree(0);
101 std::vector<double>* nn_frame =
new vector<double>;
102 std::vector<double>*
dt =
new vector<double>;
112 for (
int iu=0; iu<
NPIX; iu++) deltat[iu]=0.;
116 nn_otree->Branch(
"nnframe", &nn_frame, NDIM*NDIM, 0);
117 nn_otree->Branch(
"dt_corepixels",&dt,NPIX,0);
118 nn_otree->Branch(
"index",&ind,
"ind/I");
119 nn_otree->Branch(
"t_duration",&dT,
"duration/D");
120 nn_otree->Branch(
"f_duration",&dF,
"f_duration/D");
121 nn_otree->Branch(
"central_f",&Fc,
"central_f/D");
124 int nn_isize = nn_itree->GetEntries();
125 cout <<
"nn_itree size : " << nn_isize << endl;
128 for(
int m=0;
m<nn_isize;
m++) {
129 if(
m%100==0) cout <<
"ClusterToFrame -> " <<
m <<
"/" << nn_isize << endl;
131 nn_itree->GetEntry(
m);
133 #ifdef WATPLOT // monster event display 136 sprintf(fname,
"monster_event_display_%d.png",
m);
137 cout <<
"write " << fname << endl;
138 WTS.
plot(pwc, 1, nIFO,
'L', 0, const_cast<char*>(
"COLZ"));
149 for (
int kk=0; kk<NDIM*
NDIM; kk++) (*nn_frame)[kk]=0.;
150 for (
int kk=0; kk<
NPIX; kk++) (*dt)[kk]=0.;
151 SNRfi=
ConvertToFrame(pwc,
m, nIFO,
'L', 0, dT, dF, Fc,ninf, nn_frame,dt);
154 #ifdef PLOT_FRAME // frame display 163 cout<<
"nn_size: "<<nn_isize<<endl;
169 cout<<
"number events with SNRfi<"<<
REF<<
": "<<ninf<<endl;
170 cout<<
"media importanza ultimo bin: "<<SNRf/nn_isize<<endl;
175 ConvertToFrame(
netcluster*
pwc,
int n_ev,
int nifo,
char type,
int irate,
double &
dT,
double &dF,
double &Fc,
int &ninf, std::vector<double>* nn_frame, std::vector<double>*
dt) {
178 if(type !=
'L' && type !=
'N')
return 0.;
182 std::vector<int>* vint = &(pwc->
cList[cid-1]);
184 int V = vint->size();
190 std::vector<int>* pv = pwc->
cRate.size() ? &(pwc->
cRate[cid-1]) : NULL;
191 irate = pv!=NULL ? (*pv)[0] : 0;
207 for(
int j=0;
j<V;
j++) {
209 if(!pix->
core)
continue;
211 if((irate)&&(irate !=
int(pix->
rate+0.5)))
continue;
223 double time=((double)pix->
time/rate)/pix->
layers;
224 if(time<minTime) minTime=time;
225 if(time>maxTime) {maxTime=time;
238 if(freq<minFreq) minFreq=
freq;
239 if(freq>maxFreq) maxFreq=
freq;
241 int const ncorepix=num_core;
242 double deltat[ncorepix];
243 for (
int u=0;u<ncorepix;u++) deltat[u]=0.;
250 int minRate=RATE/(pow(2,
MAXL));
251 int maxRate=RATE/(pow(2,
MINL));
256 Fc=(maxFreq+minFreq)/2;
258 double mindt = 1./maxRate;
259 double maxdt = 1./minRate;
260 double mindf = minRate/2.;
261 double maxdf = maxRate/2.;
269 int iminTime =
int(minTime);
270 int imaxTime =
int(maxTime+1);
271 int nTime = (imaxTime-iminTime)*maxRate;
274 hist2D=
new TH2F(
"WTF",
"WTF", nTime, iminTime, imaxTime, 2*
maxLayers-2, 0, RATE/2);
275 hist2D->SetXTitle(
"time, sec");
276 hist2D->SetYTitle(
"frequency, Hz");
278 double dFreq = (maxFreq-minFreq)/10.>2*maxdf ? (maxFreq-minFreq)/10. : 2*maxdf ;
279 double mFreq = minFreq-dFreq<0 ? 0 : minFreq-
dFreq;
280 double MFreq = maxFreq+dFreq>RATE/2 ? RATE/2 : maxFreq+
dFreq;
281 hist2D->GetYaxis()->SetRangeUser(mFreq, MFreq);
283 double dTime = (maxTime-minTime)/10.>2*maxdt ? (maxTime-minTime)/10. : 2*maxdt ;
284 double mTime = minTime-dTime<iminTime ? iminTime : minTime-dTime;
285 double MTime = maxTime+dTime>imaxTime ? imaxTime : maxTime+dTime;
286 hist2D->GetXaxis()->SetRangeUser(mTime,MTime);
290 double Likelihood2=0;
295 TTree* tS=
new TTree(
"Time_SNR",
"Time_SNR");
296 tS->Branch(
"amplitude",&,
"amplitude/D");
297 tS->Branch(
"central_time",&ttime,
"central_time/D");
298 for(
int n=0;
n<V;
n++) {
300 if(!pix->
core)
continue;
304 if((irate)&&(irate !=
int(pix->
rate+0.5)))
continue;
314 double null = wSNR-sSNR;
317 int iRATE = RATE/(pix->
layers-1);
324 int i=(itime-iminTime)*maxRate;
326 if(iRATE!=irate && irate!=0)
continue;
329 int L=0;
int R=1;
while (R < iRATE) {R*=2;L++;}
333 double deltaBIN=(maxTime-minTime)/(
NBIN-1);
335 TH1D*
h=
new TH1D(
"SNR_distribution",
"SNR_distribution",
NBIN,1,
NBIN+1);
337 double timec2[ncorepix];
340 double time_bin[
NBIN];
341 for (
int nbin=0;nbin<
NBIN;nbin++){
342 double mintime_bin=maxTime-(nbin)*deltaBIN;
343 if (nbin==0||nbin==NBIN-1) mintime_bin=mintime_bin-0.001;
344 time_bin[nbin]=mintime_bin;
346 sprintf(cuts,
"central_time>%5.5f",mintime_bin);
347 TTreeFormula*
treeformula=
new TTreeFormula(
"cuts",cuts,tS);
348 int err = treeformula->Compile(cuts);
350 cout <<
"DisplayROC.C - wrong input cuts " << cuts << endl;
354 tS->SetBranchAddress(
"central_time",&ttime);
355 tS->SetBranchAddress(
"amplitude",&);
356 for(
int n=0;
n<tS->GetEntries();
n++) {
361 if(treeformula!=NULL) {
362 if(treeformula->EvalInstance()==0) {
370 h->Fill(nbin+1,ybin);
373 std::sort(timec2,timec2+ncorepix);
376 TMath::Sort(ncorepix,timec2,index);
377 for (
int u=0;u<
NPIX;u++){
380 (*dt)[u] = deltat[index[u]];
389 TCanvas* h_canv=
new TCanvas(
"SNR",
"SNR",0,0,500,500);
392 sprintf(h_char,
"SNRisto/event%i.png",n_ev);
393 h_canv->SaveAs(h_char);
398 while (h->GetBinContent(js)==0&&js<
NBIN) js=js+1;
399 SNRfi0=(h->GetBinContent(js)/h->GetBinContent(NBIN));
403 double amp2[ncorepix];
404 for (
int u=0;u<ncorepix;u++) amp2[u]=0.;
406 for(
int iu=0;iu<ncorepix;iu++){
407 tS->SetBranchAddress(
"central_time",&ttime);
408 tS->SetBranchAddress(
"amplitude",&);
410 sprintf(cuts2,
"central_time>%5.5f",timec2[ncorepix-1-iu]-0.001);
411 TTreeFormula* treeformula2=
new TTreeFormula(
"cuts2",cuts2,tS);
412 int err2 = treeformula2->Compile(cuts2);
414 cout <<
"DisplayROC.C - wrong input cuts " << cuts2 << endl;
419 for(
int n=0;
n<tS->GetEntries();
n++) {
422 if(treeformula2!=NULL) {
423 if(treeformula2->EvalInstance()==0) {
433 if((ampTOT/h->GetBinContent(NBIN))>
REF) {
434 if(iu>0) ampFin=amp2[iu-1];
436 t_limit= timec2[ncorepix-1-iu];
445 SNRfi=ampFin/(h->GetBinContent(NBIN));
448 for(
int n=0;
n<V;
n++) {
450 if(!pix->
core)
continue;
451 if((irate)&&(irate !=
int(pix->
rate+0.5)))
continue;
457 if (SNRfi>0.&&SNRfi<=
REF&&(((
double)pix->
time/rate)/pix->
layers)>t_limit)
continue;
459 sSNR2 += pow(pix->
getdata(
'S',
m),2);
460 sSNR2 += pow(pix->
getdata(
'P',
m),2);
461 wSNR2 += pow(pix->
getdata(
'W',
m),2);
462 wSNR2 += pow(pix->
getdata(
'U',
m),2);
464 double null2 = wSNR2-sSNR2;
465 int iRATE = RATE/(pix->
layers-1);
470 int i=(itime-iminTime)*maxRate;
472 if(iRATE!=irate && irate!=0)
continue;
476 int L=0;
int R=1;
while (R < iRATE) {R*=2;L++;}
477 for(
int m=0;
m<
M;
m++) {
478 for(
int k=0;
k<
K;
k++) {
479 double SNRsum=
hist2D->GetBinContent(
i+1+
m,j+1+
k-K/2);
480 if(type==
'L')
hist2D->SetBinContent(
i+1+
m,j+1+
k-K/2,sSNR2+SNRsum);
481 else hist2D->SetBinContent(
i+1+
m,j+1+
k-K/2,null2+SNRsum);
486 int imin=
hist2D->GetNbinsX();
488 int jmin=
hist2D->GetNbinsY();
491 for (
int i=0;
i<=
hist2D->GetNbinsX();
i++) {
492 for (
int j=0;
j<=
hist2D->GetNbinsY();
j++) {
493 double X =
hist2D->GetXaxis()->GetBinCenter(
i)+
hist2D->GetXaxis()->GetBinWidth(
i)/2.;
494 double Y =
hist2D->GetYaxis()->GetBinCenter(
j)+
hist2D->GetYaxis()->GetBinWidth(
j)/2.;
495 double Z =
hist2D->GetBinContent(
i,
j);
510 int di = imax-imin+1;
511 int dj = jmax-jmin+1;
517 if (di%
NDIM==0) ri=0;
518 if (dj%
NDIM==0) rj=0;
576 int xdim = (imax-imin+1)/
NDIM;
577 int ydim = (jmax-jmin+1)/
NDIM;
583 for (
int i=0;
i<=
hist2D->GetNbinsX();
i++) {
585 for (
int j=0;
j<=
hist2D->GetNbinsY();
j++) {
590 double X =
hist2D->GetXaxis()->GetBinCenter(
i)+
hist2D->GetXaxis()->GetBinWidth(
i)/2.;
591 double Y =
hist2D->GetYaxis()->GetBinCenter(
j)+
hist2D->GetYaxis()->GetBinWidth(
j)/2.;
592 double Z =
hist2D->GetBinContent(
i,
j);
602 if (
i<imin && ripi==
false)
continue;
603 if (
j<jmin && ripj==
false)
continue;
605 if (
i<imin && ripi==
true) {
i=
i+xdim;ripi2=
true;}
606 if (
j<jmin && ripj==
true) {
j=
j-ydim;ripj2=
true;}
609 int I = (
i-imin)/xdim;
610 int J = (
j-jmin)/ydim;
618 if (ripi2==
true)
i=
i-xdim;
619 if (ripj2==
true)
j=
j+ydim;
630 int nn_size = nn_frame->size();
631 for(
int i=0;
i<nn_size;
i++) norm+=(*nn_frame)[
i];
632 for(
int i=0;
i<nn_size;
i++) (*nn_frame)[
i]/=norm;
633 if (SNRfi<=
REF) ninf=ninf+1;
644 int nframe = nn_frame->size();
645 if(nframe != pow(
int(sqrt(nframe)),2)) {
646 cout <<
"PlotFrame - Error : size is not a square number " << endl;
649 nframe = sqrt(nframe);
652 hist2D=
new TH2F(
"frame",
"frame", nframe, 0, nframe, nframe, 0, nframe);
656 hist2D->SetFillColor(kWhite);
658 hist2D->GetXaxis()->SetNdivisions(506);
659 hist2D->GetXaxis()->SetLabelFont(42);
660 hist2D->GetXaxis()->SetLabelOffset(0.014);
661 hist2D->GetXaxis()->SetTitleOffset(1.4);
662 hist2D->GetYaxis()->SetTitleOffset(1.2);
663 hist2D->GetYaxis()->SetNdivisions(506);
664 hist2D->GetYaxis()->SetLabelFont(42);
665 hist2D->GetYaxis()->SetLabelOffset(0.01);
666 hist2D->GetZaxis()->SetLabelFont(42);
667 hist2D->GetZaxis()->SetNoExponent(
false);
668 hist2D->GetZaxis()->SetNdivisions(506);
670 hist2D->GetXaxis()->SetTitleFont(42);
671 hist2D->GetXaxis()->SetTitle(
"Time");
672 hist2D->GetXaxis()->CenterTitle(
true);
673 hist2D->GetYaxis()->SetTitleFont(42);
674 hist2D->GetYaxis()->SetTitle(
"Frequency");
675 hist2D->GetYaxis()->CenterTitle(
true);
677 hist2D->GetZaxis()->SetTitleOffset(0.6);
678 hist2D->GetZaxis()->SetTitleFont(42);
679 hist2D->GetZaxis()->CenterTitle(
true);
681 hist2D->GetXaxis()->SetLabelSize(0.03);
682 hist2D->GetYaxis()->SetLabelSize(0.03);
683 hist2D->GetZaxis()->SetLabelSize(0.03);
685 for(
int i=0;
i<nframe;
i++) {
686 for(
int j=0;
j<nframe;
j++) {
688 hist2D->SetBinContent(
i+1,
j+1,(*nn_frame)[
i*nframe+
j]);
693 canvas=
new TCanvas(
"mra",
"mra", 200, 20, 600, 600);
695 canvas->ToggleEventStatus();
698 canvas->SetFillColor(kWhite);
699 canvas->SetRightMargin(0.10);
700 canvas->SetLeftMargin(0.10);
701 canvas->SetBottomMargin(0.13);
705 gStyle->SetFrameBorderMode(0);
708 gStyle->SetTitleH(0.050);
709 gStyle->SetTitleW(0.95);
710 gStyle->SetTitleY(0.98);
711 gStyle->SetTitleFont(12,
"D");
712 gStyle->SetTitleColor(kBlue,
"D");
713 gStyle->SetTextFont(12);
714 gStyle->SetTitleFillColor(kWhite);
715 gStyle->SetLineColor(kWhite);
716 gStyle->SetNumberContours(256);
717 gStyle->SetCanvasColor(kWhite);
718 gStyle->SetStatBorderSize(1);
724 TPaletteAxis *
palette = (TPaletteAxis*)
hist2D->GetListOfFunctions()->FindObject(
"palette");
725 palette->SetX1NDC(0.91);
726 palette->SetX2NDC(0.933);
727 palette->SetTitleOffset(0.92);
728 palette->GetAxis()->SetTickSize(0.01);
732 sprintf(fname,
"nn_frame_display_%d.png",cid);
733 cout <<
"write " << fname << endl;
std::vector< vector_int > cRate
void ClusterToFrameNew01_dtmin_approx_cent(TString ifName, TString ofName)
std::vector< vector_int > cList
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
void PlotFrame(std::vector< double > *nn_frame, int cid)
netpixel * getPixel(size_t n, size_t i)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double getdata(char type='R', size_t n=0)
float wSNR[PE_MAX_EVENTS][NIFO_MAX]
double ConvertToFrame(netcluster *pwc, int cid, int nifo, char type, int irate, double &dT, double &dF, double &Fc, int &ninf, std::vector< double > *nn_frame, std::vector< double > *dt)