Logo coherent WaveBurst  
Library Reference Guide
Logo
Test0_dANN_var.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Serena Vinciguerra
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #define XIFO 4
20 
21 #pragma GCC system_header
22 
23 #include "cwb.hh"
24 #include "config.hh"
25 #include "network.hh"
26 #include "wavearray.hh"
27 #include "TString.h"
28 #include "TObjArray.h"
29 #include "TObjString.h"
30 #include "TPaletteAxis.h"
31 #include "TMultiLayerPerceptron.h"
32 #include "TMLPAnalyzer.h"
33 #include "TRandom.h"
34 #include "TComplex.h"
35 #include "TMath.h"
36 #include "mdc.hh"
37 #include "watplot.hh"
38 #include <vector>
39 
40 
41 //#define nINP 64
42 #define nIFO 3
43 #define nRHO 3
44 #define nANN 10
45 #define deltaANN 0.1
46 #define nCC 4
47 #define NPIX 20
48 #define RHOth 5.5
49 #define nth 10
50 using namespace std;
51 void graph(TString ifile);
52 void Annth(TString ifile);
53 void Plots(TString ifile);
54 //void Test0(TString NN_FILE,TString TEST_FILE,TString ofile, int TS=0, int TB=0, int s=0, int b=0, int uf=0){
55 void Test0_dANN_var(TString NN_FILE,TString TEST_FILE, int TS=0, int TB=0, int s=0, int b=0, int uf=0){
56 
57  int p=0;
58  char NNi[1024];
59  char NNi2[1024];
60  sprintf(NNi2,"%s",NN_FILE.Data());
61  while (NNi2[p]){
62  //cout<<NNi2[p]<<endl;
63  if (NNi2[p]=='N'){
64 // cout<<p<<endl;
65  if((NNi2[p+1]=='N')&&(NNi2[p+2]=='\/')) {
66 // cout<<" dentro if"<<endl;
67  int hh=p+3;
68  while (NNi2[hh]!='\0'){NNi[hh-p-3]=NNi2[hh];hh=hh+1;}
69  break;
70  }
71  }
72  p=p+1;
73  }
74 
75  int q=0;
76  char Filei[1024];
77 // for (int i=0;i<1024;i++)Filei[i]='';
78  char Filei2[1024];
79  sprintf(Filei2,"%s",TEST_FILE.Data());
80  while (Filei2[q]){
81 // cout<<Filei2[q]<<endl;
82  if(Filei2[q]=='n'&&Filei2[q+1]=='n'&&Filei2[q+2]=='T'&&Filei2[q+3]=='R'&&Filei2[q+4]=='E'&&(Filei2[q+5]=='E')&&(Filei2[q+6]=='\/')) {
83  int hh=q+7;
84  //while (Filei2[hh]!='\0') {Filei[hh-p-8]=Filei2[hh];hh=hh+1;cout<<Filei2[hh]<<endl;}
85  while (Filei2[hh]!='\0') {Filei[hh-p-7]=Filei2[hh];hh=hh+1;
86  //cout<<Filei2[hh]<<endl;
87  }
88  for (int h0=hh-p-7;h0<1024;h0++) Filei[h0]='\0';
89  break;
90  }
91  q=q+1;
92  }
93 
94  cout<<NNi<<" original String: "<<NNi2<<endl;
95  cout<<Filei<<" original String: "<<Filei2<<endl;
96 
97  TString NNi0(NNi);
98  TString Filei0(Filei);
99 
100  NNi0.ReplaceAll(".root","");
101  Filei0.ReplaceAll(".root","");
102 
103 
104  TFile* fnet =TFile::Open(NN_FILE.Data());
105  TMultiLayerPerceptron *mlp = (TMultiLayerPerceptron*)fnet->Get("TMultiLayerPerceptron");
106  if(mlp==NULL) {cout << "Error getting mlp" << endl;exit(1);}
107  TTree* infot=(TTree*)fnet->Get("info");
108  int NNs;
109  int NNb;
110  int NNnTS;
111  int NNnTB;
112  infot->SetBranchAddress("Rand_start_Sig",&NNs);
113  infot->SetBranchAddress("Rand_start_Bg",&NNb);
114  infot->SetBranchAddress("#trainSig",&NNnTS);
115  infot->SetBranchAddress("#trainBg",&NNnTB);
116  infot->GetEntry(0);
117 
118  TFile* fTEST =TFile::Open(TEST_FILE.Data());
119  TTree* NNTree=(TTree*)fTEST->Get("nnTree");
120  int entries=NNTree->GetEntries();
121  cout<<"entries: "<<entries<<endl;
122 
123 //estraggo le info che servono anche al tree di def dell'mlp
124  int ndim;
125  int ninp;
126  int y;
127  int entriesTot;
128  int bg_entries;
129  int sig_entries;
130 
131  NNTree->SetBranchAddress("#Entries_type",&entriesTot);
132  NNTree->SetBranchAddress("Matrix_dim",&ndim);
133  NNTree->SetBranchAddress("#inputs",&ninp);
134  NNTree->SetBranchAddress("amplitude_mode",&y);
135 
136  NNTree->GetEntry(0);
137  sig_entries=entriesTot;
138  NNTree->GetEntry(entries-1);
139  bg_entries=entriesTot;
140 
141  int const NDIM=ndim;
142  int const nINP=ninp;
143 
144  cout<<"NDIM "<<NDIM<<endl;
145  cout<<"nINP"<<nINP<<endl;
146  cout<<"sig e "<<sig_entries<<endl;
147  cout<<"bg e "<<bg_entries<<endl;
148  int minevents=0;
149  if (sig_entries>bg_entries) minevents=bg_entries;
150  else minevents=sig_entries;
151 
152  if(b==0) b=sig_entries;
153  if(TB==0 && TS==0){
154  TS=minevents;
155  TB=minevents;
156  }
157  if (b<sig_entries){
158  cout<<"Error: Bg index<sig_entries"<<endl;
159  exit(0);
160  }
161  if((TS>sig_entries||TB>bg_entries)&&(TS==TB)) {TS=minevents-s;TB=minevents-b+sig_entries;}
162  if((TS>sig_entries||TB>bg_entries)&&(TS!=TB)) {TS=sig_entries-s;TB=bg_entries-b+sig_entries;}
163 
164 
165  char NOMEtot[1024];
166  sprintf(NOMEtot,"NN_%s_FILE_%s_TS%i_TB%i_st%i_bt%i_dANN%1.2f_uf%i",NNi0.Data(),Filei0.Data(),TS,TB,s,b,deltaANN,uf);
167  cout<<"nome: "<<NOMEtot<<endl;
168 
169  TString SIG_FILE;
170  TString BG_FILE;
171  char FILE_NAME[516];
172  NNTree->SetBranchAddress("Files_name",&FILE_NAME);
173  NNTree->GetEntry(0);
174  SIG_FILE=FILE_NAME;
175  NNTree->GetEntry(entries-1);
176  BG_FILE=FILE_NAME;
177  cout<<"fine ifdef RHO_CC"<<endl;
178 
179  TChain sigTree("waveburst");//cerca il Tree "waveburst nei file
180  sigTree.Add(SIG_FILE.Data());//determina i file
181  netevent signal(&sigTree,nIFO);
182  int sig_entries2 = signal.GetEntries();
183  cout << "sig entries2 : " << sig_entries2 << endl;
184 
185  TChain bgTree("waveburst");
186  bgTree.Add(BG_FILE.Data());
187  netevent background(&bgTree,nIFO);
188  int bg_entries2 = background.GetEntries();
189  cout << "bg entries2 : " << bg_entries2 << endl;
190 
191  cout<<"b: "<<b<<endl;
192  cout<<"s: "<<s<<endl;
193 
194  // add leaf
195  Float_t x[nINP];
196  for (int jj=0; jj<nINP;jj++) x[jj]=0.;
197  char ilabel[nINP][16];
198 
199  //costruzione una foglia per ogni input
200  for(int i=0;i<nINP;i++) {
201  sprintf(ilabel[i],"x%i",i+1);
202  NNTree->SetBranchAddress(ilabel[i], &x[i]);
203  }
204 
205 char ofile[1024];
206 sprintf(ofile,"outfile/%s.root",NOMEtot);
207 //TFile*f=new TFile(ofile.Data(),"RECREATE");
208 TFile*f=new TFile(ofile,"RECREATE");
209  TTree* NNTree2=new TTree("Parameters","Parameters");
210 NNTree2->SetDirectory(f);
211  double out=0.;
212  double NNcc=0.;
213  double NNrho=0.;
214  NNTree2->Branch("ANNout",&out,"ANNout/D");
215  NNTree2->Branch("cc",&NNcc,"cc/D");
216  NNTree2->Branch("rho",&NNrho,"rho/D");
217  char NNfilename[1024];
218  NNTree2->Branch("NNname",&NNfilename,"NNname/C");
219  sprintf(NNfilename,"%s",NN_FILE.Data());
220  char Testf[1024];
221  NNTree2->Branch("TestFile",&Testf,"TestFile/C");
222  sprintf(Testf,"%s",TEST_FILE.Data());
223  int nTestS;
224  NNTree2->Branch("#TestSig",&nTestS,"#TestSig/I");
225 // nTestS=TS;
226  cout<<"nTestS: "<<nTestS<<" TS: "<<TS<<endl;
227  cout<<"dopo def tree"<<endl;
228 /* TChain TreeS("waveburst");//cerca il Tree "waveburst nei file
229  TreeS.Add(TEST_FILE.Data());//determina i file
230  netevent entryS(&TreeS,nIFO);
231  int entriesTot=0;
232  entriesTot = entryS.GetEntries();
233  cout << "Sig entries : " << entriesTot << endl;
234  std::vector<double>* frameS = new vector<double>;
235  entryS.fChain->SetBranchAddress("nnframe", &frameS);//fa esattamente come ha fatto per le altre variabili nella macro netevent.cc
236 
237  //cout<<"cc"<<(double)entryS.netcc[1]<<endl;
238 */ //cout<<"rho"<<(double)entryS.rho[0]<<endl;
239 
240 int scount=0;
241  if(uf==0) nTestS=TS;
242  else {
243 for(int n=s;n<s+TS;n++) {
244  if (uf!=0&&n>=NNs&&n<=(NNs+NNnTS)) continue;
245  scount=scount+1;
246  }
247 }
248 
249  double params[nINP];
250  int sig_05=0;
251  for (int i=0; i<nINP; i++) params[i]=0.;
252  //for(int n=0;n<entryS.GetEntries();n++) {
253  for(int n=s;n<s+TS;n++) {
254  if (uf!=0&&n>=NNs&&n<=(NNs+NNnTS)) continue;
255  NNTree->GetEntry(n);
256  signal.GetEntry(n);
257  NNcc=(double)signal.netcc[1];
258  NNrho=(double)signal.rho[0];
259  for (int i=0; i<nINP; i++){
260  //x[i]=(*frameS)[i];
261  params[i]=x[i];
262  }
263  double output=mlp->Evaluate(0,params);
264  out=output;
265  if (out<0.6) sig_05=sig_05+1;
266  cout<<"rho: "<<signal.rho[0]<<" cc "<<signal.netcc[1]<<" out: "<<out<<endl;
267  NNTree2->Fill();
268  }
269 
270 cout<<"riempito sig"<<endl;
271 int bg_05=0;
272  //for(int n=sig_entries+b;n<sig_entries+b+TB;n++) {
273  for(int n=b;n<b+TB;n++) {
274  if (uf!=0&&n>=NNb&&n<=(NNb+NNnTB)) continue;
275  NNTree->GetEntry(n);
276  cout<<"n: "<<n<<"Bg index"<<(n-sig_entries)<<endl;
277  background.GetEntry(n-sig_entries);
278  NNcc=(double)background.netcc[1];
279  NNrho=(double)background.rho[0];
280  for (int i=0; i<nINP; i++){
281  //x[i]=(*frameS)[i];
282  params[i]=x[i];
283  }
284 
285  double output=mlp->Evaluate(0,params);
286  out=output;
287  if(out>0.6) bg_05=bg_05+1;
288  cout<<"rho: "<<background.rho[0]<<" cc "<<background.netcc[1]<<" out: "<<out<<endl;
289  NNTree2->Fill();
290  }
291 cout<<"riempito bg"<<endl;
292 int Realentries=0;
293 Realentries=NNTree2->GetEntries();
294 NNTree2->Write();
295 f->Close();
296 cout<<"chiuso file"<<endl;
297 
298 //graph(ofile.Data());
299 graph(ofile);
300 cout<<"dopo richiamo funzione"<<endl;
301 Annth(ofile);
302 Plots(ofile);
303 cout<<" Sig with out<06: "<<sig_05<<endl;
304 cout<<" Bg with out>06: "<<bg_05<<endl;
305 cout<<" Realentries "<<Realentries<<endl;
306 
307 }
308 
310  TString name(ifile);
311  name.ReplaceAll("outfile/","");
312  TFile* fTEST =TFile::Open(ifile.Data());
313  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
314  double out;
315  double cc;
316  double rho;
317  int nSi;
318  NNTree2->SetBranchAddress("ANNout",&out);
319  NNTree2->SetBranchAddress("cc",&cc);
320  NNTree2->SetBranchAddress("rho",&rho);
321  NNTree2->SetBranchAddress("#TestSig",&nSi);
322  NNTree2->GetEntry(0);
323  int const nSig=nSi;
324  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
325  int const ncurve=nANN*nCC;
326  cout<<"dentro funzione dopodef"<<endl;
327 //LOG(NBg)vs RHO------------------------------------------
328  double* rhoSig[ncurve];
329  for (int i=0;i<ncurve;i++) rhoSig[i]=new double[nSig];
330  cout<<"dopo def rhoSig"<<endl;
331  //double rhoSig[20][10];
332  int NSig[ncurve];
333 // int nBg=0;
334  int const nBg=NNTree2->GetEntries()-nSig;
335  for (int i=0;i<ncurve;i++) {
336  NSig[i]=0;
337  for (int j=0;j<nSig;j++) rhoSig[i][j]=0.;
338  }
339  cout<<"dopo def rhoSig"<<endl;
340  double* rhoBg[ncurve];
341  for (int i=0;i<ncurve;i++) rhoBg[i]=new double[nBg];
342  //double rhoBg[20][10];
343  int NBg[ncurve];
344  for (int i=0;i<ncurve;i++) {
345  NBg[i]=0;
346  for (int j=0;j<nBg;j++) rhoBg[i][j]=0.;
347  }
348  cout<<"dopo def rhoBg"<<endl;
349  double ccTh[nCC];
350  for (int i=0;i<nCC;i++) ccTh[i]=0.;
351  double NNTh[nANN];
352  for (int i=0;i<nANN;i++) NNTh[i]=0.;
353  double deltacc=0.;
354 // double deltaANN=0.;
355  deltacc=0.2/nCC;
356  //deltaANN=0.6/(nANN-1);
357 
358  cout<<NNTree2->GetEntries()<<endl;
359  for(int n=0;n<NNTree2->GetEntries();n++){
360  cout<<n<<endl;
361  NNTree2->GetEntry(n);
362  cout<<"rho "<<rho<<" cc "<<cc<<" out "<<out<<endl;
363  for(int i=0; i<nCC;i++){
364  ccTh[i]=0.5+i*deltacc;
365  if(cc<ccTh[i]) continue;
366  for(int j=0; j<nANN;j++){
367  if(j==0) NNTh[j]=-1000.;
368  //else NNTh[j]=0.5+(j-1)*deltaANN;
369  //else NNTh[j]=0.1+(j-1)*deltaANN;
370  else NNTh[j]=0.+(j)*deltaANN;
371  //else NNTh[j]=0.+(j-1)*deltaANN/1000.;
372  //else NNTh[j]=1.;
373  if(out<NNTh[j]) continue;
374  int ni=0;
375  if(n>nSig) {
376  NBg[i*nANN+j]= NBg[i*nANN+j]+1;
377  while(rhoBg[i*nANN+j][ni]!=0)ni=ni+1;
378  rhoBg[i*nANN+j][ni]=rho;
379  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
380 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
381  }
382  else {
383  NSig[i*nANN+j]= NSig[i*nANN+j]+1;
384  while(rhoSig[i*nANN+j][ni]!=0)ni=ni+1;
385  rhoSig[i*nANN+j][ni]=rho;
386  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
387 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
388  }
389  // }
390  }
391  }
392  }
393  cout<<"dopo riempimento variabili"<<endl;
394  int* indexS[ncurve];
395  for (int i=0;i<ncurve;i++) indexS[i]=new int[nSig];
396 
397  TGraph * gS[ncurve];
398  for (int y=0;y<ncurve;y++) {
399  int igS=0;
400  int igS_p=0;
401  /* int indexS[nSig];
402  double rhoS[nSig];
403  for(int i=0;i<nSig;i++) {
404  rhoS[i]=0.;
405  indexS[i]=0;
406  }*/
407  // ig=1;
408  // for(int i=0;i<nSig;i++) rhoS[i]=rhoSig[y][i];
409  gS[y]=new TGraph();
410 // gS[y]->SetMarkerStyle(7);
411  TMath::Sort(nSig,rhoSig[y],indexS[y],false);
412 // cout<<"dopo Sort "<<y<<endl;
413  for (int k=0;k<nSig;k++) {
414  int ii=indexS[y][k];
415  int yy=0;
416  if (k>0){
417 // cout<<"k "<<k<<endl;
418  int ij=indexS[y][k-1];
419  //if(rhoSig[y][ii]!=0&&rhoSig[y][ii]!=rhoSig[y][ij]) {
420  if(rhoSig[y][ii]!=0) {
421  yy=NSig[y]-igS;
422  //gS[y]->SetPoint(igS,rhoSig[y][ii],yy);
423  if(rhoSig[y][ii]!=rhoSig[y][ij]) gS[y]->SetPoint(igS_p++,rhoSig[y][ii],yy);
424  cout<<"igS"<<igS<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
425  igS=igS+1;
426  }
427  }
428  else {
429  if(rhoSig[y][ii]!=0){
430  yy=NSig[y]-igS;
431  if(rhoSig[y][ii]>100)continue;
432  gS[y]->SetPoint(0,rhoSig[y][ii],yy);
433  igS=igS+1;
434 // cout<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
435  }
436 
437  }
438  }
439  }
440  // cout<<"dopo inserimento puntiiS"<<endl;
441 // cout<<ncurve<<endl;
442  int* indexB[ncurve];
443  for (int i=0;i<ncurve;i++) indexB[i]=new int[nBg];
444 
445  TGraph * gB[ncurve];
446  for (int y=0;y<ncurve;y++) {
447  int igB=0;
448  int igB_p=0;
449  gB[y]=new TGraph();
450  TMath::Sort(nBg,rhoBg[y],indexB[y],false);
451 // cout<<"dopo Sort "<<y<<endl;
452  for (int k=0;k<nBg;k++) {
453  int ii=indexB[y][k];
454  int yy=0;
455  if (k>0){
456  int ij=indexB[y][k-1];
457  //if(rhoBg[y][ii]!=0&&rhoBg[y][ii]!=rhoBg[y][ij]) {
458  if(rhoBg[y][ii]!=0) {
459  yy=NBg[y]-igB;
460  //gB[y]->SetPoint(igB,rhoBg[y][ii],yy);
461  if(rhoBg[y][ii]!=rhoBg[y][ij]) gB[y]->SetPoint(igB_p++,rhoBg[y][ii],yy);
462  igB=igB+1;
463 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
464  }
465  }
466  else {
467  if(rhoBg[y][ii]!=0) {
468  yy=NBg[y]-igB;
469  gB[y]->SetPoint(0,rhoBg[y][ii],yy);
470  igB=igB+1;
471 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
472  }
473  }
474  }
475  }
476  cout<<"dopo inserimento puntiB"<<endl;
477  //gB[1]->SetMarkerStyle(7);
478  //gB[1]->SetPoint(1,1,2);
479  TCanvas* cS=new TCanvas("Efficiency_vs_rho","Efficiency_vs_rho",0,0,1200,700);
480  cS->Divide(2,2);
481  cS->cd(1)->SetLogy();
482  TMultiGraph* mg1=new TMultiGraph();
483 // gS[0]->SetMarkerStyle(7);
484  gS[0]->SetMarkerColor(2);
485  gS[0]->SetLineColor(2);
486  mg1->SetTitle("cc=0.5;rho;#Events");
487  if(gS[1]->GetN()!=0) mg1->Add(gS[1]);
488  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
489 // gS[0]->Draw("apl");
490  for (int h=2;h<nANN;h++){
491  gS[h]->SetMarkerColor(3);
492  gS[h]->SetLineColor(3);
493 // gS[h]->SetMarkerStyle(h+1);
494  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
495  // gS[h]->Draw("apl,same");
496  }
497  mg1->Draw("apl");
498  cS->cd(2)->SetLogy();
499  TMultiGraph* mg2=new TMultiGraph();
500 // gS[nANN]->SetMarkerStyle(7);
501  gS[nANN]->SetMarkerColor(2);
502  gS[nANN]->SetLineColor(2);
503  mg2->SetTitle("cc=0.55;rho;#Events");
504  if(gS[nANN]->GetN()!=0) mg2->Add(gS[nANN]);
505  for (int h=1;h<nANN;h++){
506  gS[nANN+h]->SetMarkerColor(3);
507  gS[nANN+h]->SetLineColor(3);
508 // gS[nANN+h]->SetMarkerStyle(h+1);
509  if(gS[nANN+h]->GetN()!=0) mg2->Add(gS[nANN+h]);
510  // gS[nANN+h]->Draw("apl,same");
511  }
512  mg2->Draw("apl");
513  cS->cd(3)->SetLogy();
514  TMultiGraph* mg3=new TMultiGraph();
515 // gS[nANN*2]->SetMarkerStyle(7);
516  gS[nANN*2]->SetMarkerColor(2);
517  gS[nANN*2]->SetLineColor(2);
518  mg3->SetTitle("cc=0.6;rho;#Events");
519  if(gS[nANN*2+1]->GetN()!=0) mg3->Add(gS[nANN*2+1]);
520  if(gS[nANN*2]->GetN()!=0) mg3->Add(gS[nANN*2]);
521  for (int h=2;h<nANN;h++){
522  gS[2*nANN+h]->SetMarkerColor(3);
523  gS[2*nANN+h]->SetLineColor(3);
524 // gS[2*nANN+h]->SetMarkerStyle(h+1);
525  if(gS[2*nANN+h]->GetN()!=0) mg3->Add(gS[2*nANN+h]);
526  // gS[2*nANN+h]->Draw("apl,same");
527  }
528  mg3->Draw("apl");
529  cS->cd(4)->SetLogy();
530  TMultiGraph* mg4=new TMultiGraph();
531 // gS[nANN*3]->SetMarkerStyle(7);
532  gS[nANN*3]->SetMarkerColor(2);
533  gS[nANN*3]->SetLineColor(2);
534  mg4->SetTitle("cc=0.65;rho;#Events");
535  if(gS[nANN*3]->GetN()!=0) mg4->Add(gS[nANN*3]);
536  // gS[nANN*3]->Draw("apl");
537  for (int h=1;h<nANN;h++){
538  gS[3*nANN+h]->SetMarkerColor(3);
539  gS[3*nANN+h]->SetLineColor(3);
540 // gS[3*nANN+h]->SetMarkerStyle(h+1);
541  if(gS[3*nANN+h]->GetN()!=0) mg4->Add(gS[3*nANN+h]);
542  // gS[3*nANN+h]->Draw("apl,same");
543  }
544  mg4->Draw("apl");
545  cout<<"nuovo canv"<<endl;
546  TCanvas* cB=new TCanvas("Number_vs_rho","Number_vs_rho",0,0,1200,700);
547  cB->Divide(2,2);
548  cB->cd(1)->SetLogy();
549  TMultiGraph* mg1B=new TMultiGraph();
550 // gB[0]->SetMarkerStyle(7);
551  gB[0]->SetMarkerColor(2);
552  gB[0]->SetLineColor(2);
553  mg1B->SetTitle("cc=0.5;rho;#Events");
554  if(gB[0]->GetN()!=0) mg1B->Add(gB[0]);
555  for (int h=1;h<nANN;h++){
556  gB[h]->SetMarkerColor(3);
557  gB[h]->SetLineColor(3);
558  // gB[h]>SetMarkerStyle(h+1);
559  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
560  //gB[h]->Draw("apl,same");
561  }
562  mg1B->Draw("apl");
563  cB->cd(2)->SetLogy();
564  TMultiGraph* mg2B=new TMultiGraph();
565  //gB[nANN]->SetMarkerStyle(7);
566  gB[nANN]->SetMarkerColor(2);
567  gB[nANN]->SetLineColor(2);
568  mg2B->SetTitle("cc=0.55;rho;#Events");
569  if(gB[nANN]->GetN()!=0) mg2B->Add(gB[nANN]);
570  for (int h=1;h<nANN;h++){
571  gB[nANN+h]->SetMarkerColor(3);
572  gB[nANN+h]->SetLineColor(3);
573  //gB[nANN+h]>SetMarkerStyle(h+1);
574  if(gB[nANN+h]->GetN()!=0) mg2B->Add(gB[nANN+h]);
575  //gB[h]->Draw("apl,same");
576  }
577  mg2B->Draw("apl");
578  cB->cd(3)->SetLogy();
579  TMultiGraph* mg3B=new TMultiGraph();
580 // gB[2*nANN]->SetMarkerStyle(7);
581  gB[2*nANN]->SetMarkerColor(2);
582  gB[2*nANN]->SetLineColor(2);
583  mg3B->SetTitle("cc=0.6;rho;#Events");
584  if(gB[2*nANN]->GetN()!=0) mg3B->Add(gB[2*nANN]);
585  for (int h=1;h<nANN;h++){
586  gB[2*nANN+h]->SetMarkerColor(3);
587  gB[2*nANN+h]->SetLineColor(3);
588 // gB[2*nANN+h]>SetMarkerStyle(h+1);
589  if(gB[2*nANN+h]->GetN()!=0) mg3B->Add(gB[2*nANN+h]);
590  //gB[h]->Draw("apl,same");
591  }
592  mg3B->Draw("apl");
593  cB->cd(4)->SetLogy();
594  TMultiGraph* mg4B=new TMultiGraph();
595 // gB[3*nANN]->SetMarkerStyle(7);
596  gB[3*nANN]->SetMarkerColor(2);
597  gB[3*nANN]->SetLineColor(2);
598  mg4B->SetTitle("cc=0.65;rho;#Events");
599  if(gB[3*nANN]->GetN()!=0) mg4B->Add(gB[3*nANN]);
600  for (int h=1;h<nANN;h++){
601  gB[3*nANN+h]->SetMarkerColor(3);
602  gB[3*nANN+h]->SetLineColor(3);
603 // gB[3*nANN+h]>SetMarkerStyle(h+1);
604  if(gB[3*nANN+h]->GetN()!=0) mg4B->Add(gB[3*nANN+h]);
605  //gB[h]->Draw("apl,same");
606  }
607  mg4B->Draw("apl");
608 
609 // cout<<"dopo Draw()"<<endl;
610  TString CnameS(name);
611  TString CnameB(name);
612  TString CnameSroot(name);
613  TString CnameBroot(name);
614  char CnameS2[1024];
615  char CnameB2[1024];
616  char CnameS2root[1024];
617  char CnameB2root[1024];
618  CnameS.ReplaceAll(".root",".png");
619  CnameB.ReplaceAll(".root",".png");
620  sprintf(CnameS2,"logN_rho/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameS.Data());
621  sprintf(CnameB2,"logN_rho/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameB.Data());
622  sprintf(CnameS2root,"logN_rho/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameSroot.Data());
623  sprintf(CnameB2root,"logN_rho/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameBroot.Data());
624  cS->SaveAs(CnameS2);
625  cB->SaveAs(CnameB2);
626  cS->Print(CnameS2root);
627  cB->Print(CnameB2root);
628 // cout<<"fine"<<endl;
629 
630 }
631 
632 
634  TString name(ifile);
635  name.ReplaceAll("outfile/","");
636  TFile* fTEST =TFile::Open(ifile.Data());
637  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
638  double out;
639  double cc;
640  double rho;
641  int nSi;
642  NNTree2->SetBranchAddress("ANNout",&out);
643  NNTree2->SetBranchAddress("cc",&cc);
644  NNTree2->SetBranchAddress("rho",&rho);
645  NNTree2->SetBranchAddress("#TestSig",&nSi);
646  NNTree2->GetEntry(0);
647  int const nSig=nSi;
648 // cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
649  int const ncurve2=nRHO*nCC;
650 
651 
652  double* ANNSig[ncurve2];
653  for (int i=0;i<ncurve2;i++) ANNSig[i]=new double[nSig];
654  //double rhoSig[20][10];
655  int NSig[ncurve2];
656 // int nBg=0;
657  int const nBg=NNTree2->GetEntries()-nSig;
658  for (int i=0;i<ncurve2;i++) {
659  NSig[i]=0;
660  for (int j=0;j<nSig;j++) ANNSig[i][j]=0.;
661  }
662  double* ANNBg[ncurve2];
663  for (int i=0;i<ncurve2;i++) ANNBg[i]=new double[nBg];
664 
665  //double rhoBg[20][10];
666  int NBg[ncurve2];
667  for (int i=0;i<ncurve2;i++) {
668  NBg[i]=0;
669  for (int j=0;j<nBg;j++) ANNBg[i][j]=0.;
670  }
671  double ccTh[nCC];
672  for (int i=0;i<nCC;i++) ccTh[i]=0.;
673  double rhoTh[nRHO];
674  for (int i=0;i<nRHO;i++) rhoTh[i]=0.;
675  double deltacc=0.;
676  double deltarho=0.;
677  deltacc=0.2/nCC;
678  deltarho=1./(nRHO);
679 
680 
681  for(int n=0;n<NNTree2->GetEntries();n++){
682  NNTree2->GetEntry(n);
683  cout<<"rho "<<rho<<" cc "<<cc<<" out "<<out<<endl;
684  for(int i=0; i<nCC;i++){
685  ccTh[i]=0.5+i*deltacc;
686  if(cc<ccTh[i]) continue;
687  for(int j=0; j<nRHO;j++){
688  rhoTh[j]=5+j*deltarho;
689  if(rho<rhoTh[j]) continue;
690  int ni=0;
691  if(n>nSig) {
692  NBg[i*nRHO+j]= NBg[i*nRHO+j]+1;
693  while(ANNBg[i*nRHO+j][ni]!=0)ni=ni+1;
694  ANNBg[i*nRHO+j][ni]=out;
695  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
696  }
697  else {
698  NSig[i*nRHO+j]= NSig[i*nRHO+j]+1;
699  while(ANNSig[i*nRHO+j][ni]!=0)ni=ni+1;
700  ANNSig[i*nRHO+j][ni]=out;
701  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
702  }
703  }
704  }
705  }
706 
707  int* indexS[ncurve2];
708  for (int i=0;i<ncurve2;i++) indexS[i]=new int[nSig];
709 
710  TGraph * gS[ncurve2];
711  for (int y=0;y<ncurve2;y++) {
712  int igS=0;
713  int igS_p=0;
714  gS[y]=new TGraph();
715  TMath::Sort(nSig,ANNSig[y],indexS[y],false);
716  for (int k=0;k<nSig;k++) {
717  int ii=indexS[y][k];
718  int yy=0;
719  if (k>0){
720  //cout<<"k "<<k<<endl;
721  int ij=indexS[y][k-1];
722  //if(ANNSig[y][ii]!=0&&ANNSig[y][ii]!=ANNSig[y][ij]) {
723  if(ANNSig[y][ii]!=0) {
724  yy=NSig[y]-igS;
725  //gS[y]->SetPoint(igS,ANNSig[y][ii],yy);
726  if(ANNSig[y][ii]!=ANNSig[y][ij]) gS[y]->SetPoint(igS_p++,ANNSig[y][ii],yy);
727  // cout<<"igS"<<igS<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
728  igS=igS+1;
729 
730  }
731  }
732  else {
733  if(ANNSig[y][ii]!=0){
734  yy=NSig[y]-igS;
735  gS[y]->SetPoint(0,ANNSig[y][ii],yy);
736  igS=igS+1;
737  // cout<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
738  }
739 
740  }
741  }
742  }
743 
744 
745 
746  int* indexB[ncurve2];
747  for (int i=0;i<ncurve2;i++) indexB[i]=new int[nBg];
748 
749  TGraph * gB[ncurve2];
750  for (int y=0;y<ncurve2;y++) {
751  int igB=0;
752  int igB_p=0;
753  gB[y]=new TGraph();
754  TMath::Sort(nBg,ANNBg[y],indexB[y],false);
755  // cout<<"dopo Sort "<<y<<endl;
756  for (int k=0;k<nBg;k++) {
757  int ii=indexB[y][k];
758  int yy=0;
759  if (k>0){
760  int ij=indexB[y][k-1];
761  //if(ANNBg[y][ii]!=0&&ANNBg[y][ii]!=ANNBg[y][ij]) {
762  if(ANNBg[y][ii]!=0) {
763  yy=NBg[y]-igB;
764  //gB[y]->SetPoint(igB,ANNBg[y][ii],yy);
765  if(ANNBg[y][ii]!=ANNBg[y][ij]) gB[y]->SetPoint(igB_p++,ANNBg[y][ii],yy);
766  igB=igB+1;
767  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
768  }
769  }
770  else {
771  if(ANNBg[y][ii]!=0) {
772  yy=NBg[y]-igB;
773  gB[y]->SetPoint(0,ANNBg[y][ii],yy);
774  igB=igB+1;
775  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
776  }
777  }
778  }
779  }
780 
781 
782  TCanvas* cS=new TCanvas("Efficiency_vs_ANN","Efficiency_vs_ANN",0,0,1200,700);
783  cS->Divide(2,2);
784  //cS->cd(1)->SetLogy();
785  cS->cd(1);
786  TMultiGraph* mg1=new TMultiGraph();
787  mg1->SetTitle("cc=0.5;ANN;#Events");
788  for (int h=0;h<nRHO;h++){
789  gS[h]->SetLineColor(4);
790  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
791  }
792  mg1->Draw("al");
793  cS->cd(2);
794  // cS->cd(2)->SetLogy();
795  TMultiGraph* mg2=new TMultiGraph();
796  mg2->SetTitle("cc=0.55;ANN;#Events");
797  for (int h=0;h<nRHO;h++){
798  gS[nRHO+h]->SetLineColor(4);
799  if(gS[nRHO+h]->GetN()!=0) mg2->Add(gS[nRHO+h]);
800  }
801  mg2->Draw("al");
802 
803  //cS->cd(3)->SetLogy();
804  cS->cd(3);
805  TMultiGraph* mg3=new TMultiGraph();
806  mg3->SetTitle("cc=0.6;ANN;#Events");
807  for (int h=0;h<nRHO;h++){
808  gS[2*nRHO+h]->SetLineColor(4);
809  if(gS[2*nRHO+h]->GetN()!=0) mg3->Add(gS[2*nRHO+h]);
810  }
811  mg3->Draw("al");
812  //cS->cd(4)->SetLogy();
813  cS->cd(4);
814  TMultiGraph* mg4=new TMultiGraph();
815  mg4->SetTitle("cc=0.65;ANN;#Events");
816  for (int h=0;h<nRHO;h++){
817  gS[3*nRHO+h]->SetLineColor(4);
818  if(gS[3*nRHO+h]->GetN()!=0) mg4->Add(gS[3*nRHO+h]);
819  }
820  mg4->Draw("al");
821 
822  TCanvas* cB=new TCanvas("Number_vs_ANN","Number_vs_ANN",0,0,1200,700);
823  cB->Divide(2,2);
824  cB->cd(1)->SetLogy();
825  TMultiGraph* mg1B=new TMultiGraph();
826  mg1B->SetTitle("cc=0.5;ANN;#Events");
827  for (int h=0;h<nRHO;h++){
828  gB[h]->SetLineColor(4);
829  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
830  }
831  mg1B->Draw("al");
832  cB->cd(2)->SetLogy();
833  TMultiGraph* mg2B=new TMultiGraph();
834  mg2B->SetTitle("cc=0.55;ANN;#Events");
835  for (int h=0;h<nRHO;h++){
836  gB[nRHO+h]->SetLineColor(4);
837  if(gB[nRHO+h]->GetN()!=0) mg2B->Add(gB[nRHO+h]);
838  }
839  mg2B->Draw("al");
840  cB->cd(3)->SetLogy();
841  TMultiGraph* mg3B=new TMultiGraph();
842  mg3B->SetTitle("cc=0.6;ANN;#Events");
843  for (int h=0;h<nRHO;h++){
844  gB[2*nRHO+h]->SetLineColor(4);
845  if(gB[2*nRHO+h]->GetN()!=0) mg3B->Add(gB[2*nRHO+h]);
846  }
847  mg3B->Draw("al");
848  cB->cd(4)->SetLogy();
849  TMultiGraph* mg4B=new TMultiGraph();
850  mg4B->SetTitle("cc=0.6;ANN;#Events");
851  for (int h=0;h<nRHO;h++){
852  gB[3*nRHO+h]->SetLineColor(4);
853  if(gB[3*nRHO+h]->GetN()!=0) mg4B->Add(gB[3*nRHO+h]);
854  }
855  mg4B->Draw("al");
856 
857 
858  //cout<<"dopo Draw()"<<endl;
859  TString CnameS(name);
860  TString CnameB(name);
861  TString CnameSroot(name);
862  TString CnameBroot(name);
863  char CnameS2[1024];
864  char CnameB2[1024];
865  char CnameS2root[1024];
866  char CnameB2root[1024];
867  CnameS.ReplaceAll(".root",".png");
868  CnameB.ReplaceAll(".root",".png");
869  sprintf(CnameS2,"ANNthres/N_ANN_S_%s",CnameS.Data());
870  sprintf(CnameB2,"ANNthres/N_ANN_B_%s",CnameB.Data());
871  sprintf(CnameS2root,"ANNthres/N_ANN_S_%s",CnameSroot.Data());
872  sprintf(CnameB2root,"ANNthres/N_ANN_B_%s",CnameBroot.Data());
873  cS->SaveAs(CnameS2);
874  cB->SaveAs(CnameB2);
875  cS->Print(CnameS2root);
876  cB->Print(CnameB2root);
877  //cout<<"fine"<<endl;
878 
879 //CARTELLA Annth
880 
881 
882 
883 }
884 
885 
887  TString name(ifile);
888  name.ReplaceAll("outfile/","");
889  TFile* fTEST =TFile::Open(ifile.Data());
890  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
891  double out;
892  double cc;
893  double rho;
894  int nSi;
895  NNTree2->SetBranchAddress("ANNout",&out);
896  NNTree2->SetBranchAddress("cc",&cc);
897  NNTree2->SetBranchAddress("rho",&rho);
898  NNTree2->SetBranchAddress("#TestSig",&nSi);
899  NNTree2->GetEntry(0);
900  int const nSig=nSi;
901  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
902  int const nBg=NNTree2->GetEntries()-nSig;
903 
904 /* TH2D* hglitch=new TH2D("Purezza","Purezza",NCC,minCC,maxCC,NRHO,minRHO,maxRHO);
905  TH2D* h2BgO_R=new TH2D("rho_cc_out","rho_cc_out",50,0.,50,50,-0.5,1.5);
906  TH2D* h2SigO_R=new TH2D("rho_cc_out","rho_cc_out",50,0.,50,50,-0.5,1.5);
907  TH2D* h2BgR_C=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,0.,50);
908  TH2D* h2SigR_C=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,0.,50);
909  TH2D* h2Bg=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,-0.5,1.5);
910  TH2D* h2Sig=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,-0.5,1.5);
911  h2BgO_R->SetDirectory(0);
912  h2SigO_R->SetDirectory(0);
913  h2BgR_C->SetDirectory(0);
914  h2SigR_C->SetDirectory(0);
915  h2Bg->SetDirectory(0);
916  h2Sig->SetDirectory(0);
917  hglitch->SetDirectory(0);
918 */
919 
920  TGraph * gS[3];
921  TGraph * gB[3];
922  gB[0]=new TGraph();
923  gS[0]=new TGraph();
924  gB[1]=new TGraph();
925  gS[1]=new TGraph();
926  gB[2]=new TGraph();
927  gS[2]=new TGraph();
928  double aANNS[nSig];
929  double aRHOS[nSig];
930  double aCCS[nSig];
931  for (int i=0;i<nSig;i++){
932  aCCS[i]=0.;
933  aRHOS[i]=0.;
934  aANNS[i]=0.;
935  }
936  double aANNB[nBg];
937  double aRHOB[nBg];
938  double aCCB[nBg];
939  for (int i=0;i<nBg;i++){
940  aCCB[i]=0.;
941  aRHOB[i]=0.;
942  aANNB[i]=0.;
943  }
944 // cout<<"nSig: "<<nSig<<endl;
945  for (int n = 0; n <NNTree2->GetEntries(); n++){
946  NNTree2->GetEntry(n);
947  if(n<nSig) {
948  aRHOS[n]=rho;
949  aANNS[n]=out;
950  aCCS[n]=cc;
951 /* gS[0]->SetPoint(n,cc,rho);
952  gS[1]->SetPoint(n,cc,out);
953  gS[2]->SetPoint(n,rho,out);
954  cout<<"Sig_graph1: x="<<cc<<" y: "<<rho<<endl;
955  cout<<"Sig_graph2: x="<<cc<<" y: "<<out<<endl;
956  cout<<"Sig_graph3: x="<<rho<<" y: "<<out<<endl;
957 */
958  }
959  else {
960  aRHOB[n-nSig]=rho;
961  aANNB[n-nSig]=out;
962  aCCB[n-nSig]=cc;
963 /* gB[0]->SetPoint(n-nSig,cc,rho);
964  gB[1]->SetPoint(n-nSig,cc,out);
965  gB[2]->SetPoint(n-nSig,rho,out);
966  cout<<"Bg_graph1: x="<<cc<<" y: "<<rho<<endl;
967  cout<<"Bg_graph2: x="<<cc<<" y: "<<out<<endl;
968  cout<<"Bg_graph3: x="<<rho<<" y: "<<out<<endl;
969 */
970  }
971  }
972 
973  int indexB[nBg];
974  for (int k=0;k<nBg;k++)indexB[k]=0;
975  TMath::Sort(nBg,aRHOB,indexB,false);
976  int Z[(nth+1)*nth];
977  for (int k=0;k<(nth+1)*nth;k++)Z[k]=0;
978  double cc1th[nth];
979  double ANN1th[nth+1];
980  ANN1th[0]=0.;
981  for (int k=0;k<(nth);k++){
982  ANN1th[k+1]=0.;
983  cc1th[k]=0.;
984  }
985  double dANNth=0.;
986  double dccth=0.;
987  dANNth=1./nth;
988  dccth=0.5/nth;
989  ANN1th[0]=-100.;
990  for (int k=0;k<(nth);k++){
991  ANN1th[k+1]=0.+dANNth*(k+1);
992  cc1th[k]=0.5+dccth*k;
993  }
994 
995  /*for (int a=0;a<nBg;a++){
996  if(aRHOB[a]>=RHOth){
997  for (int b=0;b<nth;b++){
998  if(aCCB[a]>=cc1th[b]){
999  for(int c=0;c<nth+1;c++){
1000  if(aANNB[a]>=ANN1th[c]) Z[b*nth+c]=Z[b*nth+c]+1;
1001  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1002  }
1003  }
1004  }
1005  }
1006  }*/
1007  for (int a=0;a<nBg;a++){
1008  cout<<"tutti: a "<<a<<" rho "<<aRHOB[a]<<" aNNB "<<aANNB[a]<<" cc "<<aCCB[a]<<endl;
1009  if(aRHOB[a]<RHOth)continue;
1010  cout<<"dopo rho: a "<<a<<" rho "<<aRHOB[a]<<" aNNB "<<aANNB[a]<<" cc "<<aCCB[a]<<endl;
1011  for (int b=0;b<nth;b++){
1012  if(aCCB[a]<cc1th[b])continue;
1013  cout<<"dopo cc: "<<cc1th[b]<<"a "<<a<<" rho "<<aRHOB[a]<<" aNNB "<<aANNB[a]<<" cc "<<aCCB[a]<<endl;
1014  for(int c=0;c<nth+1;c++){
1015  if(aANNB[a]<ANN1th[c]) continue;
1016  cout<<"dopo ANN:"<<ANN1th[c]<<" a "<<a<<" rho "<<aRHOB[a]<<" aNNB "<<aANNB[a]<<" cc "<<aCCB[a]<<endl;
1017  Z[b*(nth+1)+c]=Z[b*(nth+1)+c]+1;
1018  cout<<"Zindex: "<<b*(nth+1)+c<<" Z: "<<Z[b*(nth+1)+c]<<endl;
1019 
1020  //if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*(nth+1)+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1021  }
1022 
1023  }
1024 
1025  }
1026 
1027  //TH2D* hglitch=new TH2D("Survived_glitches(RHO_cut:5.5)","Survived_glitches(RHO_cut:5.5)",nth,0.5,1,nth+1,-1,1);
1028  double bins[nth+1];
1029  for (int d=0; d<nth+1;d++) bins[d]=0.;
1030  //bins[0]=0.;
1031  for (int d=0; d<nth+1;d++) bins[d+1]=0.+(d+1)*dANNth;
1032  TH2D* hglitch=new TH2D("Survived_glitches(RHO_cut:5.5)","Survived_glitches(RHO_cut:5.5)",nth,0.5,1,nth+1,bins);
1033 // TH2D* hglitch=new TH2D("Survived_glitches(RHO_cut:5.5)","Survived_glitches(RHO_cut:5.5)",nth,0.5,1,nth+1,0.,1.);
1034  hglitch->SetDirectory(0);
1035  //if(tot[NRHO*ii+ji]>0) hglitch->Fill(cc1[ji]+deltacc/2,rho1[ii]+deltarho/2,pur[NRHO*ii+ji]);
1036  for (int b=0;b<nth;b++) for(int c=0;c<nth+1;c++) {
1037  //hglitch->Fill(cc1th[b]+dccth/2.,ANN1th[c]+dANNth/2.,Z[(nth+1)*b+c]);
1038  hglitch->Fill(cc1th[b]+dccth/2.,bins[c]+dANNth/2.,Z[(nth+1)*b+c]);
1039  cout<<" Zindex: "<<(nth+1)*b+c<< " x: "<<cc1th[b]<<" y: "<<ANN1th[c]<<" z: "<<Z[(nth+1)*b+c]<<endl;
1040  }
1041 
1042  gB[0]=new TGraph(nBg,aCCB,aRHOB);
1043  gS[0]=new TGraph(nSig,aCCS,aRHOS);
1044  gB[1]=new TGraph(nBg,aCCB,aANNB);
1045  gS[1]=new TGraph(nSig,aCCS,aANNS);
1046  gB[2]=new TGraph(nBg,aRHOB,aANNB);
1047  gS[2]=new TGraph(nSig,aRHOS,aANNS);
1048 
1049  for(int i=0;i<3;i++){
1050  gS[i]->SetMarkerColor(2);
1051  gB[i]->SetMarkerColor(4);
1052  gS[i]->SetMarkerStyle(6);
1053  gB[i]->SetMarkerStyle(7);
1054  }
1055 
1056  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1057  c->Divide(2,2);
1058  c->cd(1);
1059  TMultiGraph* mg1=new TMultiGraph();
1060  mg1->SetTitle("cc_rho;cc;rho");
1061  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1062  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1063  mg1->Draw("ap");
1064  c->cd(2);
1065  TMultiGraph* mg2=new TMultiGraph();
1066  mg2->SetTitle("cc_ANNoutput;cc;ANNoutput");
1067  if(gB[1]->GetN()!=0) mg2->Add(gB[1]);
1068  if(gS[1]->GetN()!=0) mg2->Add(gS[1]);
1069  mg2->Draw("ap");
1070  c->cd(3);
1071  TMultiGraph* mg3=new TMultiGraph();
1072  mg3->SetTitle("rho_ANNoutput;rho;ANNoutput");
1073  if(gB[2]->GetN()!=0) mg3->Add(gB[2]);
1074  if(gS[2]->GetN()!=0) mg3->Add(gS[2]);
1075  mg3->Draw("ap");
1076  c->cd(4);
1077  TText* text=new TText(0.37,0.0,"no cuts on ANN");
1078  hglitch->SetStats(0);
1079  hglitch->GetXaxis()->SetTitle("cc");
1080  hglitch->GetYaxis()->SetTitle("ANNoutput");
1081  hglitch->GetZaxis()->SetTitle("count");
1082  hglitch->Draw("colz");
1083  text->Draw();
1084  gPad->SetLogz();
1085 
1086  TString Cname(name);
1087  TString Cnameroot(name);
1088  char Cname2[1024];
1089  char Cname2root[1024];
1090  Cname.ReplaceAll(".root",".png");
1091  sprintf(Cname2,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_%s",Cname.Data());
1092  sprintf(Cname2root,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_%s",Cnameroot.Data());
1093  c->SaveAs(Cname2);
1094  c->Print(Cname2root);
1095 
1096  TCanvas* c2=new TCanvas("Plots_Bkg_on_Sig","Plots_Bkg_on_Sig",0,0,1200,700);
1097  c2->Divide(2,2);
1098  c2->cd(1);
1099  gS[0]->Draw("ap");
1100  gB[0]->Draw("p,same");
1101 
1102  c2->cd(2);
1103  gS[1]->Draw("ap");
1104  gB[1]->Draw("p,same");
1105  c2->cd(3);
1106  gS[2]->Draw("ap");
1107  gB[2]->Draw("p,same");
1108  c2->cd(4);
1109  TText* text2=new TText(0.37,0.0,"no cuts on ANN");
1110  hglitch->SetStats(0);
1111  hglitch->GetXaxis()->SetTitle("cc");
1112  hglitch->GetYaxis()->SetTitle("ANNoutput");
1113  hglitch->GetZaxis()->SetTitle("count");
1114  hglitch->Draw("colz");
1115  text2->Draw();
1116  gPad->SetLogz();
1117 
1118  TString Cname_2(name);
1119  TString Cname_2root(name);
1120  char Cname2_2[1024];
1121  char Cname2_2root[1024];
1122  Cname_2.ReplaceAll(".root",".png");
1123  sprintf(Cname2_2,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2.Data());
1124  sprintf(Cname2_2root,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2root.Data());
1125  c2->SaveAs(Cname2_2);
1126  c2->Print(Cname2_2root);
1127 
1128 //CARTELLA CCi_RHO_ANN_Plots
1129 }
1130 
char ofile[1024]
double rho
void graph(TString ifile)
Float_t * rho
biased null statistics
Definition: netevent.hh:112
TCanvas * c2
wavearray< double > a(hp.size())
par [0] name
int n
Definition: cwb_net.C:28
TString("c")
ofstream out
Definition: cwb_merge.C:214
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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
Int_t GetEntry(Int_t)
Definition: netevent.cc:409
CWB::Toolbox TB
STL namespace.
int j
Definition: cwb_net.C:28
i drho i
cout<< fileout<< endl;out.open(fileout, ios::out);if(!out.good()) { cout<< "cwb_mkhtml_cbc.C : Error Opening File : "<< fileout<< endl;exit(1);} out<< "<br><br><< endl; MakePlotsHtmlCellTable(&out, "Sensitive distance for M< sub > total</sub > bins
wavearray< double > hh
Definition: Regression_H1.C:73
void Plots(TString ifile)
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
#define nth
TText * text
wavearray< double > h
Definition: Regression_H1.C:25
cout<< "Injected signals: "<< mdc.GetEntries()<< endl;cout<< "Injected signals in histogram factor_events_inj: "<< NEVTS<< endl;float myifar, ecor, m1, m2, netcc[3], neted, penalty;float rho[2];float chirp[6];float range[2];float frequency[2];float iSNR[3], sSNR[3];sim.SetBranchAddress("mass", mass);sim.SetBranchAddress("factor", &factor);sim.SetBranchAddress("range", range);sim.SetBranchAddress("chirp", chirp);sim.SetBranchAddress("rho", rho);sim.SetBranchAddress("netcc", netcc);sim.SetBranchAddress("neted", &neted);sim.SetBranchAddress("ifar", &myifar);sim.SetBranchAddress("ecor", &ecor);sim.SetBranchAddress("penalty", &penalty);sim.SetBranchAddress("time", mytime);sim.SetBranchAddress("iSNR", iSNR);sim.SetBranchAddress("sSNR", sSNR);sim.SetBranchAddress("spin", spin);sim.SetBranchAddress("frequency", frequency);float **volume=new float *[NBINS_mass1];float **volume_first_shell=new float *[NBINS_mass1];float **radius=new float *[NBINS_mass1];float **error_volume=new float *[NBINS_mass1];float **error_volume_first_shell=new float *[NBINS_mass1];float **error_radius=new float *[NBINS_mass1];for(int i=0;i< NBINS_mass1;i++) { volume[i]=new float[NBINS_mass2];volume_first_shell[i]=new float[NBINS_mass2];radius[i]=new float[NBINS_mass2];error_volume[i]=new float[NBINS_mass2];error_volume_first_shell[i]=new float[NBINS_mass2];error_radius[i]=new float[NBINS_mass2];for(int j=0;j< NBINS_mass2;j++) { volume[i][j]=0.;volume_first_shell[i][j]=0.;radius[i][j]=0.;error_volume[i][j]=0.;error_volume_first_shell[i][j]=0.;error_radius[i][j]=0.;} } float **spin_mtot_volume=new float *[NBINS_MTOT+1];float **spin_mtot_radius=new float *[NBINS_MTOT+1];float **error_spin_mtot_volume=new float *[NBINS_MTOT+1];float **error_spin_mtot_radius=new float *[NBINS_MTOT+1];for(int i=0;i< NBINS_MTOT+1;i++) { spin_mtot_volume[i]=new float[NBINS_SPIN+1];spin_mtot_radius[i]=new float[NBINS_SPIN+1];error_spin_mtot_volume[i]=new float[NBINS_SPIN+1];error_spin_mtot_radius[i]=new float[NBINS_SPIN+1];for(int j=0;j< NBINS_SPIN+1;j++) { spin_mtot_volume[i][j]=0.;error_spin_mtot_volume[i][j]=0.;spin_mtot_radius[i][j]=0.;error_spin_mtot_radius[i][j]=0.;} } char fname[1024];sprintf(fname, "%s/recovered_signals.txt", netdir);ofstream fev;fev.open(fname, std::ofstream::out);sprintf(line, "#GPS@L1 FAR[Hz] eFAR[Hz] Pval " "ePval factor rho frequency iSNR sSNR \");fev<< line<< endl;ofstream *fev_single=new ofstream[nfactor];for(int l=1;l< nfactor+1;l++) { sprintf(fname, "%s/recovered_signals_%d.txt", netdir, l);fev_single[l - 1].open(fname, std::ofstream::out);fev_single[l - 1]<< line<< endl;} double Vrho[RHO_NBINS], eVrho[RHO_NBINS], Rrho[RHO_NBINS], eRrho[RHO_NBINS], Trho[RHO_NBINS];for(int i=0;i< RHO_NBINS;i++) { Vrho[i]=0.;eVrho[i]=0.;Rrho[i]=0.;eRrho[i]=0.;Trho[i]=RHO_MIN+i *RHO_BIN;} double dV, dV1, dV_spin_mtot, nevts, internal_volume;int nT;int countv=0;int cnt=0;int cnt2=0;int cntfreq=0;bool bcut=false;double liveTot=sim.GetMaximum("ifar");double BKG_LIVETIME_yr=liveTot/CYS;double BKG_LIVETIME_Myr=BKG_LIVETIME_yr/(1.e6);cout.precision(14);cout<< "Total live time ---> background
#define nIFO
char output[256]
int k
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:113
TFile * ifile
wavearray< double > yy
Definition: TestFrame5.C:12
s s
Definition: cwb_net.C:155
#define nRHO
void Annth(TString ifile)
#define deltarho
#define nCC
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
Int_t GetEntries()
Definition: netevent.cc:403
#define RHOth
void Test0_dANN_var(TString NN_FILE, TString TEST_FILE, int TS=0, int TB=0, int s=0, int b=0, int uf=0)
wavearray< double > y
Definition: Test10.C:31
#define nANN
#define deltacc
#define deltaANN
exit(0)