Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Average_FAout_moreGraphs_3.C
Go to the documentation of this file.
1 #define XIFO 4
2 
3 #pragma GCC system_header
4 
5 #include "cwb.hh"
6 #include "config.hh"
7 #include "network.hh"
8 #include "wavearray.hh"
9 #include "TString.h"
10 #include "TObjArray.h"
11 #include "TObjString.h"
12 #include "TPaletteAxis.h"
13 #include "TMultiLayerPerceptron.h"
14 #include "TMLPAnalyzer.h"
15 #include "TRandom.h"
16 #include "TComplex.h"
17 #include "TMath.h"
18 #include "mdc.hh"
19 #include "watplot.hh"
20 #include <vector>
21 
22 
23 //#define nINP 64
24 #define Ncc3D 50
25 #define maxCC 1
26 #define minCC 0
27 #define nIFO 3
28 #define nRHO 3
29 #define nANN 10
30 #define deltaANN 0.1
31 #define deltacc 0.05
32 #define deltarho 1
33 #define nCC 4
34 #define NPIX 20
35 #define RHOth 5.5
36 #define nth 10
37 #define rhomin 6
38 #define ccmin 0.6
39 #define iMc 1
40 #define NMc_inj 50
41 #define minMc_inj 0
42 #define maxMc_inj 100
43 #define NMc_r 50
44 #define maxMc_r 40
45 #define minMc_r -5
46 #define NANN_av 50
47 #define minANN_av -0.2
48 #define maxANN_av 1.5
49 
50 //definition parameters for SNR_plots
51 #define N_ANNth 50
52 #define N_RHOth 50
53 #define ANNth_max 1.5
54 #define ANNth_min -0.5
55 #define RHOth_max 10.0
56 #define RHOth_min 4.5
57 
58 using namespace std;
59 void graph(TString ifile);
60 void Annth(TString ifile, int &FAth_n, int &FDth_n, double &thresFA);
61 void Mcth(TString ifile, int &McFAth_n, int &McFDth_n, double &McthresFA);
62 void PlotsAv_cc(TString ifile);
63 void PlotsAv_Mc(TString ifile);
64 void SNR_plots(TString ifile);
65 void CoG_plots(TString ifile);
66 void duration_plot(TString ifile);
67 //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){
68 
69 //NN_FILE-NN_FILE8: name of the files where the interested ANNs are saved (!!they must become from a directory whose name ends in "NN").
70 //TEST_FILE: name of the file which contains the events you are going to analyse (!!they must become from a directory whose name ends in "nnTREE").
71 //TS, TB: numbers of events to test for each class.
72 //s, b: index values of the first events considered for both the classes. If b < #SIG_events, b is newly set to b=b+#SIG_events
73 //uf !=0 means the Test file of events is the same used to train the network, in this case ni==0. It means that the events used for the training are analysed only by the other ANNs, in this way we test the events with at least NNn-1 networks. If uf==0 the Test and Train files are different, in this case ni!=0 and all the events can beanalysed by all the nNN networks.
74 //consider_all=1(0): the events not used for the training procedure are (aren't) considered for the test.
75 //ni: =0 excludes from the test the training set of events
76 void Average_FAout_moreGraphs_3(double FAdes, TString NN_FILE,TString NN_FILE2, TString NN_FILE3, TString NN_FILE4, TString NN_FILE5, TString NN_FILE6, TString NN_FILE7, TString NN_FILE8, TString TEST_FILE, TString NOMEtot_S, int TS=0, int TB=0, int s=0, int b=0, int uf=1, int consider_all=0, int av_on_nNN=0){
77 int ni=0;
78 if (uf==0&&av_on_nNN!=0) cout<<"all events were not used for the training procedure"<<endl;
79 if(uf==0) {ni=1; consider_all=1; av_on_nNN=0;}
80 else ni=0; //events used for the training not considered for the ANN trained on them
81 
82 int TS0=0;
83 TS0=TS;
84 int TB0=0;
85 TB0=TB;
86 //if(uf!=0&&consider_all==0) ni=1;
87 
88 //COUNT THE ANNs
89 int n_NN=0;
90 if (NN_FILE.CompareTo("")){
91  n_NN=n_NN+1;
92  }
93 if (NN_FILE2.CompareTo("")){
94  n_NN=n_NN+1;
95  }
96 if (NN_FILE3.CompareTo("")){
97  n_NN=n_NN+1;
98  }
99 if (NN_FILE4.CompareTo("")){
100  n_NN=n_NN+1;
101  }
102 if (NN_FILE5.CompareTo("")){
103  n_NN=n_NN+1;
104  }
105 if (NN_FILE6.CompareTo("")){
106  n_NN=n_NN+1;
107  }
108 if (NN_FILE7.CompareTo("")){
109  n_NN=n_NN+1;
110  }
111 if (NN_FILE8.CompareTo("")){
112  n_NN=n_NN+1;
113  }
114 int const nNN=n_NN;
115 
116 // int p=0;
117  char NNi[nNN][1024];
118  char NNi2[nNN][1024];
119  sprintf(NNi2[0],"%s",NN_FILE.Data());
120  if(nNN>1) sprintf(NNi2[1],"%s",NN_FILE2.Data());
121  if(nNN>2) sprintf(NNi2[2],"%s",NN_FILE3.Data());
122  if(nNN>3) sprintf(NNi2[3],"%s",NN_FILE4.Data());
123  if(nNN>4) sprintf(NNi2[4],"%s",NN_FILE5.Data());
124  if(nNN>5) sprintf(NNi2[5],"%s",NN_FILE6.Data());
125  if(nNN>6) sprintf(NNi2[6],"%s",NN_FILE7.Data());
126  if(nNN>7) sprintf(NNi2[7],"%s",NN_FILE8.Data());
127 
128 if(nNN==1) {av_on_nNN=1; consider_all=0; ni=0;}
129 //saving the ANN names
130  for (int u=0;u<nNN;u++){
131  int p=0;
132  while (NNi2[u][p]){
133  if (NNi2[u][p]=='N'){
134  if((NNi2[u][p+1]=='N')&&(NNi2[u][p+2]=='\/')) {
135  int hh=p+3;
136  while (NNi2[u][hh]!='\0'){NNi[u][hh-p-3]=NNi2[u][hh];hh=hh+1;}
137  break;
138  }
139  }
140  p=p+1;
141  }
142 }
143 
144 //saving the Test-file name
145 int q=0;
146  char Filei[1024];
147  char Filei2[1024];
148  sprintf(Filei2,"%s",TEST_FILE.Data());
149  while (Filei2[q]){
150  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]=='\/')) {
151  int hh=q+7;
152  while (Filei2[hh]!='\0') {Filei[hh-q-7]=Filei2[hh];hh=hh+1;
153  }
154  for (int h0=hh-q-7;h0<1024;h0++) Filei[h0]='\0';
155  break;
156  }
157  q=q+1;
158  }
159 
160  cout<<Filei<<" original String: "<<Filei2<<endl;
161 
162 
163  TString NNi0[nNN];
164  NNi0[0]=NNi[0];
165  if(nNN>1) NNi0[1]=NNi[1];
166  if(nNN>2) NNi0[2]=NNi[2];
167  if(nNN>3) NNi0[3]=NNi[3];
168  if(nNN>4) NNi0[4]=NNi[4];
169  if(nNN>5) NNi0[5]=NNi[5];
170  if(nNN>6) NNi0[6]=NNi[6];
171  if(nNN>7) NNi0[7]=NNi[7];
172 
173  TString Filei0(Filei);
174 
175 //removing the extensions
176 for (int u=0;u<nNN;u++){
177  NNi0[u].ReplaceAll(".root","");
178  }
179  Filei0.ReplaceAll(".root","");
180 
181 //Definitions of some parameters
182  TFile *fnet[nNN];
183  TMultiLayerPerceptron* mlp[nNN];
184  TTree* infot[nNN];
185  int TotBg[nNN];
186  int FD0[nNN]; //FD=False Dismissal
187  int FA0[nNN]; //FA=False Alarm
188 
189  for (int yy=0; yy<nNN; yy++){TotBg[yy]=0;FA0[yy]=0;FD0[yy]=0;}
190 
191  int NNs[nNN]; //first SIG-event considered for ANN trainings
192  int NNb[nNN]; //first BKG-event considered for ANN trainings
193  int NNnTS[nNN]; //number of SIG-events considered for ANN trainings
194  int NNnTB[nNN]; //number of BKG-events considered for ANN trainings
195 
196 int min_NNb=5000000;
197 int min_NNs=5000000;
198 int max_NNb=0;
199 int max_NNs=0;
200 
201 //Extracting ANNs and information from the files
202 //if (uf!=0&&ni==0&&(TS!=0||TB!=0)){
203  for (int u=0;u<nNN;u++){
204  fnet[u]=TFile::Open(NNi2[u]);
205  mlp[u]=(TMultiLayerPerceptron*)fnet[u]->Get("TMultiLayerPerceptron");
206 
207  cout<<"dopo TMLP"<<endl;
208  if(mlp[u]==NULL) {cout << "Error getting mlp" << endl;exit(1);}
209  infot[u]=(TTree*)fnet[u]->Get("info");
210 
211  infot[u]->SetBranchAddress("Rand_start_Sig",&NNs[u]);
212  infot[u]->SetBranchAddress("Rand_start_Bg",&NNb[u]);
213  infot[u]->SetBranchAddress("#trainSig",&NNnTS[u]);
214  infot[u]->SetBranchAddress("#trainBg",&NNnTB[u]);
215  infot[u]->GetEntry(0);
216  cout<<"n "<<u<<" b "<<NNb[u]<<" min_NNb "<<min_NNb<<endl;
217 
218  if(NNs[u]<min_NNs) min_NNs=NNs[u];
219  if(NNb[u]<min_NNb) min_NNb=NNb[u];
220  if(NNs[u]>max_NNs) max_NNs=NNs[u];
221  if(NNb[u]>max_NNb) max_NNb=NNb[u];
222  cout<<"n "<<u<<" s "<<NNs[u]<<" b "<<NNb[u]<<" s fin "<<NNs[u]+NNnTS[u]<<" b fin "<<NNb[u]+NNnTB[u]<<" nTS "<<NNnTS[u]<<" nTB "<<NNnTB[u]<<endl;
223  }
224 
225 //Changing b value to its sum to the first event considered for the training
226 //if (uf!=0&&ni==0&&(b<min_NNb)&&(TS!=0||TB!=0)){ b+=min_NNb;cout<<" change in b value to have BKG events: "<<b<<endl;}
227 /*if (uf!=0&&ni==0&&(b<min_NNb)){ b+=min_NNb;cout<<" change in b value to have BKG events: "<<b<<endl;}
228 if (uf!=0&&ni==0&&(s<min_NNs)){ s+=min_NNs;cout<<" change in s value to have SIG events: "<<s<<endl;}*/
229 //if (uf!=0&&consider_all==0&&av_on_nNN==0&&(b<min_NNb)){ b+=min_NNb;cout<<" change in b value to have BKG events: "<<b<<endl;}
230 //if (uf!=0&&consider_all==0&&av_on_nNN==0&&(s<min_NNs)){ s+=min_NNs;cout<<" change in s value to have SIG events: "<<s<<endl;}
231 
232 //if (uf!=0&&av_on_nNN!=0&&(b>min_NNb)&&(b<))
233 cout<<" min_NNb "<<min_NNb<<" b "<<b<<endl;
234 cout<<"Def. tree e mlp"<<endl;
235 
236 //opening the Test-file
237  TFile* fTEST =TFile::Open(TEST_FILE.Data());
238 cout<<"Def. tree e mlp"<<endl;
239  TTree* NNTree=(TTree*)fTEST->Get("nnTree");
240 cout<<"Def. tree e mlp"<<endl;
241  int entries=NNTree->GetEntries();
242 // int entries=0;
243 
244  cout<<"entries: "<<entries<<endl;
245  cout<<" NNTree->GetEntries() "<<NNTree->GetEntries()<<endl;
246 //exit(0);
247 //Extract information necessary for the tree which defines the mlp
248  int ndim;
249  int ninp;
250  int y;
251  int entriesTot;
252  int bg_entries;
253  int sig_entries;
254 
255  NNTree->SetBranchAddress("#Entries_type",&entriesTot);
256  NNTree->SetBranchAddress("Matrix_dim",&ndim);
257  NNTree->SetBranchAddress("#inputs",&ninp);
258  NNTree->SetBranchAddress("amplitude_mode",&y);
259 
260  NNTree->GetEntry(0);
261  sig_entries=entriesTot;
262  NNTree->GetEntry(entries-1);
263  bg_entries=entriesTot;
264  cout<<"sig "<<sig_entries<<" bg "<<bg_entries<<endl;
265  int const NDIM=ndim;
266  int const nINP=ninp;
267  cout<<"NDIM: "<<NDIM<<endl;
268  cout<<"nINP: "<<nINP<<endl;
269  cout<<"sig e: "<<sig_entries<<endl;
270  cout<<"bg e: "<<bg_entries<<endl;
271 
272 //defining the sample which has less events
273  int minevents=0;
274  if (sig_entries>bg_entries) minevents=bg_entries;
275  else minevents=sig_entries;
276 
277  //if(b==0) b=sig_entries;
278  if(TB==0 && TS==0){
279  TS=minevents;
280  TB=minevents;
281  }
282 cout<<"TS "<<TS<<" TB "<<TB<<" minevents "<<minevents<<endl;
283 //defining compatible parameters
284  if (b<sig_entries){
285  int a =b;
286  b+=sig_entries;
287  cout<<"Error: Bg index<sig_entries: new set of b parameter-> b="<<b<<" instead of b="<<a<<endl;
288  //exit(0);
289 }
290  if(s>sig_entries){
291  int a=s;
292  s-=sig_entries;
293  cout<<"Error: Sig index>sig_entries: new set of s parameter-> s="<<b<<" instead of s="<<a<<endl;
294 }
295 if((TS>sig_entries-s||TB>(bg_entries-(b-sig_entries)))&&(TS==TB)) {TS=minevents-s;TB=minevents-(b-sig_entries);
296  if (TS<TB) TB=TS;
297  else TS=TB;
298  cout<<"Error:S>sig_entries or TB>bg_entries, to maintain ugual number of analysed events TS=TB="<<TB<<endl;
299  }
300  if((TS>sig_entries-s||TB>bg_entries-(b-sig_entries))&&(TS!=TB)) {
301  if(TS>sig_entries) TS=sig_entries-s;
302  if (TB>bg_entries) TB=bg_entries-(b-sig_entries);
303  cout<<"Error:S>sig_entries or TB>bg_entries, new TS and TB values are thus define-> TS="<<TS<<" TB="<<TB<<endl;
304  }
305 
306 
307  char NOMEtot[1024];
308  sprintf(NOMEtot,"%s",NOMEtot_S.Data());
309  cout<<"output name: "<<NOMEtot<<endl;
310 
311 //extracting original files
312  TString SIG_FILE;
313  TString BG_FILE;
314  char FILE_NAME[516];
315  NNTree->SetBranchAddress("Files_name",&FILE_NAME);
316  NNTree->GetEntry(0);
317  SIG_FILE=FILE_NAME;
318  NNTree->GetEntry(entries-1);
319  BG_FILE=FILE_NAME;
320  cout<<"fine ifdef RHO_CC"<<endl;
321 
322 //extracting other information
323  TChain sigTree("waveburst");//search the Tree "waveburst" in files
324  sigTree.Add(SIG_FILE.Data());
325  netevent signal(&sigTree,nIFO);
326  int sig_entries2 = signal.GetEntries();
327  cout << "sig entries2 : " << sig_entries2 << endl;
328 
329  TChain bgTree("waveburst");
330  bgTree.Add(BG_FILE.Data());
331  netevent background(&bgTree,nIFO);
332  int bg_entries2 = background.GetEntries();
333  cout << "bg entries2 : " << bg_entries2 << endl;
334 
335  cout<<"b: "<<b<<endl;
336  cout<<"s: "<<s<<endl;
337 
338  // add leaf
339  Float_t x[nINP];
340  for (int jj=0; jj<nINP;jj++) x[jj]=0.;
341  char ilabel[nINP][16];
342 
343  // define a branch for each input
344  for(int i=0;i<nINP;i++) {
345  sprintf(ilabel[i],"x%i",i+1);
346  NNTree->SetBranchAddress(ilabel[i], &x[i]);
347  }
348 
349 char ofile[1024];
350 sprintf(ofile,"average_file/%s.root",NOMEtot);
351 TFile*f =new TFile(ofile,"RECREATE");
352  TTree* NNTree2=new TTree("Parameters","Parameters");
353  NNTree2->SetDirectory(f);
354 double DtTree;
355  //NNTree->SetBranchAddress("t_duration", &DtTree);
356  NNTree->SetBranchAddress("duration", &DtTree);
357 // defining useful information
358  double average=0.;
359  double out[nNN];
360  for(int y=0; y<nNN; y++) out[y]=0.;
361  double Mc=0.;
362  double NNcc=0.;
363  double NNrho=0.;
364  double Mc_inj=0.;
365  double Dt=0.;
366  //Float_t Dt=0.;
367  double std=0.;
368  int ev_ind=0;
369  int CoG=0;
370  // int ev_ind;
371  int index_ev=0;
372  char TestFile[1024];
373 //adding new information to the original tree
374  NNTree2->Branch("Center_of_Gravity",&CoG,"Center_of_Gravity/I");
375  NNTree2->Branch("duration",&Dt,"duration/D");
376  //NNTree2->Branch("duration",&Dt,"duration/F");
377  NNTree2->Branch("Average",&average,"Average/D");
378  NNTree2->Branch("StandardDevaition",&std,"StandardDeviation/D");
379  NNTree2->Branch("cc",&NNcc,"cc/D");
380  NNTree2->Branch("Mchirp",&Mc,"Mchirp/D");
381  NNTree2->Branch("Mchirp_injected",&Mc_inj,"Mchirp_injected/D");
382  NNTree2->Branch("rho",&NNrho,"rho/D");
383  NNTree2->Branch("index",&index_ev,"index/I");
384 
385 // NNTree2->Branch("event_index",&ev_ind,"event_index/I");
386  NNTree2->Branch("Test_file",&TestFile,"Test_file/C");
387 
388  for(int u=0;u<nNN;u++){
389  char NNoutl[512];
390  char NNoutl2[512];
391  char NNnamel[512];
392  char NNnamel2[512];
393  sprintf(NNoutl,"NNout%i",u);
394  sprintf(NNoutl2,"NNout%i/D",u);
395  sprintf(NNnamel,"NNname%i",u);
396  sprintf(NNnamel2,"NNname%i/C",u);
397  NNTree2->Branch(NNoutl,&out[u],NNoutl2);
398  NNTree2->Branch(NNnamel,&NNi[u],NNnamel2);
399  }
400 
401  char Testf[1024];
402  NNTree2->Branch("TestFile",&Testf,"TestFile/C");
403  sprintf(Testf,"%s",TEST_FILE.Data());
404  int nTestS=0;
405 
406 
407 //defining how many examples to use
408  int scount=0;
409  if(uf==0) nTestS=TS;
410  else {
411  //for(int n=s;n<s+TS;n++) {
412  for(int n=s;n<sig_entries;n++) {
413  int countNN=0;
414  for(int u=0;u<nNN;u++){
415  if (uf!=0&&n>NNs[u]&&n<(NNs[u]+NNnTS[u])) countNN=countNN+1;
416  }
417 
418  // if(consider_all==0&&av_on_nNN==0&&countNN==0)continue; //skip if no network is trained on the considered event
419  // if(consider_all==0&&av_on_nNN!=0&&countNN!=0) {sigskip=sigskip+1; continue; }
420 
421  if(countNN==0&&av_on_nNN!=0) scount=scount+1;
422  else{
423  //if(nNN==1&&countNN==1&&av_on_nNN!=0) scount=scount+1;
424  if(countNN==1&&consider_all==0&&av_on_nNN==0) scount=scount+1;
425  if(countNN<2&&consider_all!=0) scount=scount+1; //if you are indifferently intereseted in average on nNN or nNN-1
426  if(countNN>1){cout<<"Error: training non independent"<<n<<endl; exit(0);} //exit if the ANNs are training on dependent set of events
427  cout<<"n "<<n<<"; countNN "<<countNN<<"; scount "<< scount<<"; consider_all "<<consider_all<<"; av_on_nNN "<<av_on_nNN<<endl;
428  if(scount==TS) break;}
429  }
430  nTestS=scount;
431 }
432 cout<<" nTestS "<< nTestS<<" TS "<<TS<<endl;
433 //exit(0);
434 int B_eff=0;
435 int FA=0;
436 int bg_05=0;
437 int cont_su=0;
438 int cont_su5=0;
439 int cont_su4=0;
440 int cont_su3=0;
441 
442 int nTestB=0;
443  int bcount=0;
444  if(uf==0) nTestB=TB;
445  else {
446  //for(int n=b;n<b+TB;n++) {
447  for(int n=b;n<sig_entries+bg_entries;n++) {
448  int countNN=0;
449  for(int u=0;u<nNN;u++){
450  if (uf!=0&&n>NNb[u]&&n<(NNb[u]+NNnTB[u])) countNN=countNN+1;
451  }
452  if(countNN==0&&av_on_nNN!=0) bcount=bcount+1;
453  //if(nNN==1&&countNN==1&&av_on_nNN!=0) bcount=bcount+1;
454  if(countNN==1&&consider_all==0&&av_on_nNN==0) bcount=bcount+1;
455  if(countNN<2&&consider_all!=0) bcount=bcount+1; //if you are indifferently intereseted in average on nNN or nNN-1
456  if(countNN>1){cout<<"Error: training non independent"<<n<<endl; exit(0);} //exit if the ANNs are training on dependent set of events
457  if(bcount==TB) break;
458  }
459  nTestB=bcount;
460 }
461 
462 if(uf!=0) nTestS=nTestS-1;
463 if(uf!=0) nTestB=nTestB-1;
464 //TB=TB-1;
465 //TS=TS-1;
466 cout<<"TS "<<TS<<" TB "<<TB<<" nTestS "<<nTestS<<" nTestB "<<nTestB<<endl;
467 if(TS0==0&&TB0==0) {
468  if(nTestB>nTestS) {nTestB=nTestS; TB=nTestS; TS=nTestS;}
469  else {nTestS=nTestB; TB=nTestB; TS=nTestB;}
470  }
471 
472 if(nTestS<TS) TS= nTestS;
473 else nTestS=TS;
474 
475 // NNTree2->Branch("#TestSig",&TS,"#TestSig/I");
476  NNTree2->Branch("#TestSig",&TS,"#TestSig/I");
477  //NNTree2->Branch("#TestSig",&nTestS,"#TestSig/I");
478  cout<<"nTestS: "<<nTestS<<" TS: "<<TS<<endl;
479  cout<<"dopo def tree"<<endl;
480 cout<<"s test "<<s<< " s+TS "<<s+TS<<" b "<<b<<" b+ TB "<<b+TB<<" nTestS "<<nTestS<<" TB "<<TB<<" nTestB "<<nTestB<<endl;
481 double params[nINP];
482  int sig_05=0;
483  for (int i=0; i<nINP; i++) params[i]=0.;
484 
485  int S_eff=0;
486  int FD=0;
487 int sigskip=0;
488 
489  TString txt_originaldata(NOMEtot_S);
490  txt_originaldata.ReplaceAll(".root","_originalData.txt");
491  //txt_originaldata.ReplaceAll(".root","_originalData.txt");
492  char txt_originaldata_c[1024];
493  sprintf(txt_originaldata_c,"txt_files/%s.txt",txt_originaldata.Data());
494  ofstream od_txt(txt_originaldata_c);
495 
496 
497  //for(int n=s;n<s+TS;n++) {
498  for(int n=s;n<sig_entries;n++) {
499 double ti_0=0;
500 double fi_0=0;
501 
502  index_ev=n;
503  average=0.;
504  int indNN=-1;
505  int countNN=0;
506  for(int u=0;u<nNN;u++){
507  if (uf!=0&&n>=NNs[u]&&n<(NNs[u]+NNnTS[u])) { indNN=u;countNN=countNN+1;}
508  }
509  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
510  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
511  if(consider_all==0&&av_on_nNN==0&&countNN==0)continue; //skip if no network is trained on the considered event
512  if(consider_all==0&&av_on_nNN!=0&&countNN!=0) {sigskip=sigskip+1; continue; }
513  NNTree->GetEntry(n);
514  signal.GetEntry(n);
515  char line_od[1024];
516 
517  // exracting parameters
518  NNcc=(double)signal.netcc[0];
519  NNrho=(double)signal.rho[1];
520  Mc=(double)signal.chirp[iMc];
521  Mc_inj=(double)signal.chirp[0];
522  Dt=(double)signal.duration[0];
523  //Dt=DtTree;
524  cout<<"DtTree "<<DtTree<<endl;
525  //Dt=(double)signal.duration[0];
526  sprintf(TestFile,"%s",TEST_FILE.Data());
527 // cout<<"background.duration[0] "<<signal.duration[0]<<" background.duration[1] "<<signal.duration[1]<<" background.duration[2] "<<signal.duration[2]<<endl;
528  //exit(0);
529  sprintf(line_od,"%i \n",n);
530  od_txt<<line_od;
531  char line_od2[1024];
532  // evaluating the event
533  for (int i=0; i<nINP; i++){
534  params[i]=x[i];
535 
536  int ti_i=0;
537  int fi_i=0;
538  ti_i=i/NDIM;
539  fi_i=i-ti_i*NDIM;
540  ti_0+=(double)x[i]*ti_i;
541  fi_0+=(double)x[i]*fi_i;
542  sprintf(line_od2,"%s %1.5f",line_od2, x[i]);
543  }
544  sprintf(line_od2,"%s \n",line_od2);
545 
546  if (ti_0<fi_0) CoG=1;
547  else CoG=0;
548 
549 //cout<<" ti_0 "<<ti_0<<" fi_0 "<<fi_0<<" CoG "<<CoG<<endl;
550  od_txt<<line_od2;
551  for (int u=0;u<nNN;u++) {
552 // cout<<" u "<<u<<" nNN "<<nNN<<endl;
553  double output=mlp[u]->Evaluate(0,params);
554 // cout<<output<<endl;
555  out[u]=output;
556  }
557 
558 // cout<<"Dopo param"<<endl;
559 
560 //calculating the average and the standard deviation
561 
562  int c_nNN0=0;
563  for (int u=0;u<nNN;u++){
564  if((u!=indNN&&ni==0)||uf==0) {average=out[u]+average;std+=out[u]*out[u];c_nNN0=c_nNN0+1; if(out[u]<0.6) FD0[u]=FD0[u]+1;}
565  }
566 
567  average=average/(c_nNN0);
568  if((c_nNN0)!=1)std=pow((std/(nNN-1-countNN)-average*average),0.5);
569  else std=-1;
570 
571  char line_od3[1024];
572  sprintf(line_od3,"ev_ind %i average %1.5f Mc %1.5f Mc_inj %1.5f n %i \n",index_ev, average,Mc,Mc_inj,n);
573  od_txt<<line_od3;
574 
575  if (average<0.6) FD=FD+1;
576  NNTree2->Fill();
577  S_eff=S_eff+1;
578  if(S_eff==1) s=n;
579  if(S_eff==TS) break;
580  cout<<" cc "<<NNcc<<" NNrho "<<NNrho<<" Mc_inj "<< Mc_inj<<endl;
581 }
582 //exit(0);
583 od_txt.close();
584 //exit(0);
585 cout<<txt_originaldata_c<<endl;
586 //exit(0);
587 //closing the cycle on SIG events
588 cout<<" S_eff "<<S_eff<<" s "<<s<<" b "<<b<<" sigskip "<<sigskip<<endl;
589 //exit(0);
590 //defining BKG parameters
591 
592 int bgskip=0;
593 //opening the cycle on BKG events
594  for(int n=b;n<sig_entries+bg_entries;n++) {
595 double ti_0=0;
596 double fi_0=0;
597  index_ev=0;
598  index_ev=n;
599  average=0.;
600  int indNN=-1;
601  int countNN=0;
602  for(int u=0;u<nNN;u++){
603  cout<<" uf "<<uf<<" n "<<n<<" u "<<u<<" NN b[u] "<<NNb[u]<<" NNb[u]+NNnTB[u]"<<NNb[u]+NNnTB[u]<<endl;
604  if (uf!=0&&n>=NNb[u]&&n<(NNb[u]+NNnTB[u])) { indNN=u;countNN=countNN+1; }
605  }
606  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
607  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
608  cout<<"consider_all "<<consider_all<<" av_on_nNN "<<av_on_nNN<<" countNN "<<countNN<<" indNN "<<indNN<<endl;
609  if(consider_all==0&&av_on_nNN==0&&countNN==0)continue; //skip if no network is trained on the considered event
610  if(consider_all==0&&av_on_nNN!=0&&countNN!=0) { cout<<" n skipped "<<n<<endl; bgskip=bgskip+1; continue;}
611  NNTree->GetEntry(n);
612  cout<<"BKG->n: "<<n<<"Bg index"<<(n-sig_entries)<<endl;
613  background.GetEntry(n-sig_entries);
614  //NNcc=(double)background.netcc[1];
615  NNcc=(double)background.netcc[0];
616  //NNrho=(double)background.rho[0];
617  NNrho=(double)background.rho[1];
618  Mc=(double)background.chirp[iMc];
619  Mc_inj=(double)background.chirp[0];
620  //Dt=(double)background.t_duration;
621  //Dt=DtTree;
622  Dt=(double)background.duration[0];
623  cout<<"background.duration[0] "<<background.duration[0]<<" background.duration[1] "<<background.duration[1]<<" background.duration[2] "<<background.duration[2]<<endl;
624  /*int hn=0;
625  ev_ind=0;
626  hn=n;
627  ev_ind=hn;*/
628  sprintf(TestFile,"%s",TEST_FILE.Data());
629 
630 for (int i=0; i<nINP; i++){
631  params[i]=x[i];
632  int ti_i=0;
633  int fi_i=0;
634  ti_i=i/NDIM;
635  fi_i=i-ti_i*NDIM;
636  ti_0+=(double)x[i]*ti_i;
637  fi_0+=(double)x[i]*fi_i;
638 
639  }
640  if (ti_0<fi_0) CoG=1;
641  else CoG=0;
642 
643  for(int u=0;u<nNN;u++) {
644  double output=mlp[u]->Evaluate(0,params);
645  cout<<output<<endl;
646  out[u]=output;
647  }
648 
649  int c_nNN0=0;
650  int nnsu=0;
651  for (int u=0;u<nNN;u++){
652  if((u!=indNN&&ni==0)||uf==0) {average=out[u]+average;std+=out[u]*out[u];c_nNN0=c_nNN0+1; TotBg[u]=TotBg[u]+1;if(out[u]>0.6) {FA0[u]=FA0[u]+1;nnsu=nnsu+1;}}
653  }
654 
655 
656  average=average/(c_nNN0);
657  if((c_nNN0)!=1)std=pow((std/(nNN-1-countNN)-average*average),0.5);
658  else std=-1;
659 
660  if(nnsu==6) cont_su=cont_su+1;
661  if(nnsu>4) cont_su5=cont_su5+1;
662  if(nnsu>3) cont_su4=cont_su4+1;
663  if(nnsu>2) cont_su3=cont_su3+1;
664 
665 //calculating the final average value
666 
667 
668 //saving false alarms
669  if(average>0.6) FA=FA+1;
670 
671  NNTree2->Fill();
672  B_eff=B_eff+1;
673  if(B_eff==1) b=n;
674  if(B_eff==TB) break;
675  }
676 //cout<<"bkg filled"<<endl;
677 
678 NNTree2->Write();
679 f->Close();
680 cout<<"B_eff "<<B_eff<<" TB "<<TB<<endl;
681 //cout<<"closed file"<<endl;
682 cout<<ofile<<endl;
683 
684 int cont=0;
685 
686 double freq_c[nNN];
687 double meanf=0.;
688 double dev=0.;
689 
690 for (int yy=0; yy<nNN; yy++){
691  freq_c[yy]=0.;
692 // cout<<" FA0[yy] "<<FA0[yy]<<" FD0[yy] "<<FD0[yy]<<" FD "<<FD<<" TotBg[yy] "<<TotBg[yy]<<endl;
693  if(TotBg[yy]!=0){freq_c[yy]=(double)FA0[yy]/TotBg[yy];
694  //cout<<" freq_c[yy] "<<freq_c[yy]<<endl;
695  }
696  if(freq_c[yy]!=0) {meanf=freq_c[yy]+meanf;cont=cont+1;
697  //cout<<" cont "<<cont<<endl;
698  }
699  dev=freq_c[yy]*freq_c[yy]+dev;
700 }
701 
702 if(cont==0) cout<<"Error cont==0"<<endl;
703 if(cont!=0) meanf=meanf/cont;
704 dev=0.;
705 for (int yy=0; yy<nNN; yy++){
706 //cout<<" dev "<<dev<<endl;
707 //cout<<"freq_c[yy] "<<freq_c[yy]<<" meanf "<<meanf<<endl;
708 if(freq_c[yy]!=0)dev+=(freq_c[yy]-meanf)*(freq_c[yy]-meanf);
709 }
710 double McFAdes=FAdes;
711 int FAth_n=FAdes*B_eff;
712 int McFAth_n=FAdes*B_eff;
713 cout<<"FAth_n "<<FAth_n<<" FAdes*B_eff "<<FAdes*B_eff<<endl;
714 int FDth_n=-1;
715 int McFDth_n=-1;
716 double thresFA=-1;
717 double McthresFA=-1;
718 Annth(ofile, FAth_n, FDth_n, thresFA);
719 Mcth(ofile, McFAth_n, McFDth_n, McthresFA);
720 graph(ofile);
721 PlotsAv_cc(ofile);
722 PlotsAv_Mc(ofile);
723 SNR_plots(ofile);
724 CoG_plots(ofile);
725 duration_plot(ofile);
726 
727 if(cont!=0) dev=pow(dev/(cont-1),0.5);
728 cout<<"meanf "<<meanf<<" dev "<<dev<<endl;
729 cout<<" FA average "<<FA<<endl;
730 cout<<" cont_su "<<cont_su<<" cont_su5 "<<cont_su5<<" cont_su4 "<<cont_su4<<" cont_su3 "<<cont_su3<<endl;
731 
732 cout<<" S_eff "<<S_eff<<" B_eff "<<B_eff<<" s "<<s<<" b "<<b<<" bgskip "<<bgskip<<" sigskip "<<sigskip<<endl;
733 cout<<" FA0[1] "<<FA0[1]<<" TotBg[1] "<<TotBg[1]<<endl;
734 cout<<" FD0[1] "<<FD0[1]<<" TotBg[1] "<<TotBg[1]<<endl;
735 //double S_eff2=(double
736 double FDdes=(double)FDth_n/S_eff;
737 double McFDdes=(double)McFDth_n/S_eff;
738 FAdes=(double)FAth_n/B_eff;
739 McFAdes=(double)McFAth_n/B_eff;
740 cout<<"FAth_n "<<FAth_n<<" FAdes: "<<FAdes<<" ANN_av_FDdes "<<FDdes<<" ANN_av_FDth_n "<<FDth_n<<" ANN_av_thresFA "<<thresFA<<endl;
741 cout<<"McFAth_n "<<McFAth_n<<" McFAdes: "<<McFAdes<<" Mc_FDdes "<<McFDdes<<" McFDth_n "<<McFDth_n<<" McthresFA "<<McthresFA<<endl;
742 }
743 
744 void graph(TString ifile){
745  TString name(ifile);
746  name.ReplaceAll("average_file/","");
747  TFile* fTEST =TFile::Open(ifile.Data());
748  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
749  double av;
750  double cc;
751  double rho;
752  int nSi;
753  NNTree2->SetBranchAddress("Average",&av);
754  NNTree2->SetBranchAddress("cc",&cc);
755  NNTree2->SetBranchAddress("rho",&rho);
756  NNTree2->SetBranchAddress("#TestSig",&nSi);
757  NNTree2->GetEntry(0);
758  int const nSig=nSi;
759  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
760  cout<<" BG: "<<NNTree2->GetEntries()-nSi<<endl;
761  int const ncurve=nANN*nCC;
762  cout<<"dentro funzione dopodef"<<endl;
763 //LOG(NBg)vs RHO------------------------------------------
764  double* rhoSig[ncurve];
765  for (int i=0;i<ncurve;i++) rhoSig[i]=new double[nSig];
766  cout<<"dopo def rhoSig"<<endl;
767  //double rhoSig[20][10];
768  int NSig[ncurve];
769 // int nBg=0;
770  int const nBg=NNTree2->GetEntries()-nSig;
771  for (int i=0;i<ncurve;i++) {
772  NSig[i]=0;
773  for (int j=0;j<nSig;j++) rhoSig[i][j]=0.;
774  }
775  cout<<"dopo def rhoSig"<<endl;
776  double* rhoBg[ncurve];
777  for (int i=0;i<ncurve;i++) rhoBg[i]=new double[nBg];
778  //double rhoBg[20][10];
779  int NBg[ncurve];
780  for (int i=0;i<ncurve;i++) {
781  NBg[i]=0;
782  for (int j=0;j<nBg;j++) rhoBg[i][j]=0.;
783  }
784  cout<<"dopo def rhoBg"<<endl;
785  double ccTh[nCC];
786  for (int i=0;i<nCC;i++) ccTh[i]=0.;
787  double NNTh[nANN];
788  for (int i=0;i<nANN;i++) NNTh[i]=0.;
789  //double deltacc=0.;
790 // double deltaANN=0.;
791 // deltacc=0.2/nCC;
792  //deltaANN=0.6/(nANN-1);
793 
794  cout<<NNTree2->GetEntries()<<endl;
795  for(int n=0;n<NNTree2->GetEntries();n++){
796  cout<<n<<endl;
797  NNTree2->GetEntry(n);
798  cout<<"rho "<<rho<<" cc "<<cc<<" av "<<av<<endl;
799  for(int i=0; i<nCC;i++){
800  ccTh[i]=0.5+i*deltacc;
801  if(cc<ccTh[i]) continue;
802  for(int j=0; j<nANN;j++){
803  if(j==0) NNTh[j]=-1000.;
804  //else NNTh[j]=0.5+(j-1)*deltaANN;
805  //else NNTh[j]=0.1+(j-1)*deltaANN;
806  else NNTh[j]=0.+(j)*deltaANN;
807  //else NNTh[j]=0.+(j-1)*deltaANN/1000.;
808  //else NNTh[j]=1.;
809  if(av<NNTh[j]) continue;
810  int ni=0;
811  if(n>=nSig) {
812 
813  NBg[i*nANN+j]= NBg[i*nANN+j]+1;
814  while(rhoBg[i*nANN+j][ni]!=0)ni=ni+1;
815  rhoBg[i*nANN+j][ni]=rho;
816  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
817 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
818  }
819  else {
820  NSig[i*nANN+j]= NSig[i*nANN+j]+1;
821  while(rhoSig[i*nANN+j][ni]!=0)ni=ni+1;
822  rhoSig[i*nANN+j][ni]=rho;
823  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
824 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
825  }
826  // }
827  }
828  }
829  }
830  cout<<"dopo riempimento variabili"<<endl;
831  int* indexS[ncurve];
832  for (int i=0;i<ncurve;i++) indexS[i]=new int[nSig];
833 
834  TGraph * gS[ncurve];
835  for (int y=0;y<ncurve;y++) {
836  int igS=0;
837  int igS_p=0;
838  /* int indexS[nSig];
839  double rhoS[nSig];
840  for(int i=0;i<nSig;i++) {
841  rhoS[i]=0.;
842  indexS[i]=0;
843  }*/
844  // ig=1;
845  // for(int i=0;i<nSig;i++) rhoS[i]=rhoSig[y][i];
846  gS[y]=new TGraph();
847 // gS[y]->SetMarkerStyle(7);
848  TMath::Sort(nSig,rhoSig[y],indexS[y],false);
849 // cout<<"dopo Sort "<<y<<endl;
850  for (int k=0;k<nSig;k++) {
851  int ii=indexS[y][k];
852  int yy=0;
853  if (k>0){
854 // cout<<"k "<<k<<endl;
855  int ij=indexS[y][k-1];
856  //if(rhoSig[y][ii]!=0&&rhoSig[y][ii]!=rhoSig[y][ij]) {
857  if(rhoSig[y][ii]!=0) {
858  yy=NSig[y]-igS;
859  //gS[y]->SetPoint(igS,rhoSig[y][ii],yy);
860  if(rhoSig[y][ii]!=rhoSig[y][ij]) gS[y]->SetPoint(igS_p++,rhoSig[y][ii],yy);
861  cout<<"igS"<<igS<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
862  igS=igS+1;
863  }
864  }
865  else {
866  if(rhoSig[y][ii]!=0){
867  yy=NSig[y]-igS;
868  gS[y]->SetPoint(0,rhoSig[y][ii],yy);
869  igS=igS+1;
870 // cout<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
871  }
872 
873  }
874  }
875  }
876  // cout<<"dopo inserimento puntiiS"<<endl;
877 // cout<<ncurve<<endl;
878  int* indexB[ncurve];
879  for (int i=0;i<ncurve;i++) indexB[i]=new int[nBg];
880 
881  TGraph * gB[ncurve];
882  for (int y=0;y<ncurve;y++) {
883  int igB=0;
884  int igB_p=0;
885  gB[y]=new TGraph();
886  TMath::Sort(nBg,rhoBg[y],indexB[y],false);
887 // cout<<"dopo Sort "<<y<<endl;
888  for (int k=0;k<nBg;k++) {
889  int ii=indexB[y][k];
890  int yy=0;
891  if (k>0){
892  int ij=indexB[y][k-1];
893  //if(rhoBg[y][ii]!=0&&rhoBg[y][ii]!=rhoBg[y][ij]) {
894  if(rhoBg[y][ii]!=0) {
895  yy=NBg[y]-igB;
896  //gB[y]->SetPoint(igB,rhoBg[y][ii],yy);
897  if(rhoBg[y][ii]!=rhoBg[y][ij]) gB[y]->SetPoint(igB_p++,rhoBg[y][ii],yy);
898  igB=igB+1;
899 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
900  }
901  }
902  else {
903  if(rhoBg[y][ii]!=0) {
904  yy=NBg[y]-igB;
905  gB[y]->SetPoint(0,rhoBg[y][ii],yy);
906  igB=igB+1;
907 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
908  }
909  }
910  }
911  }
912  cout<<"dopo inserimento puntiB"<<endl;
913  //gB[1]->SetMarkerStyle(7);
914  //gB[1]->SetPoint(1,1,2);
915  TCanvas* cS=new TCanvas("Efficiency_vs_rho","Efficiency_vs_rho",0,0,1200,700);
916  cS->Divide(2,2);
917  //cS->cd(1)->SetLogy();
918  cS->cd(1);
919  TMultiGraph* mg1=new TMultiGraph();
920 // gS[0]->SetMarkerStyle(7);
921  gS[0]->SetMarkerColor(2);
922  gS[0]->SetLineColor(2);
923  mg1->SetTitle("cc=0.5;rho;#Events");
924  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
925 // gS[0]->Draw("apl");
926  for (int h=1;h<nANN;h++){
927  gS[h]->SetMarkerColor(3);
928  gS[h]->SetLineColor(3);
929 // gS[h]->SetMarkerStyle(h+1);
930  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
931  // gS[h]->Draw("apl,same");
932  }
933  mg1->Draw("apl");
934  //cS->cd(2)->SetLogy();
935  cS->cd(2);
936  TMultiGraph* mg2=new TMultiGraph();
937 // gS[nANN]->SetMarkerStyle(7);
938  gS[nANN]->SetMarkerColor(2);
939  gS[nANN]->SetLineColor(2);
940  mg2->SetTitle("cc=0.55;rho;#Events");
941  if(gS[nANN]->GetN()!=0) mg2->Add(gS[nANN]);
942  for (int h=1;h<nANN;h++){
943  gS[nANN+h]->SetMarkerColor(3);
944  gS[nANN+h]->SetLineColor(3);
945 // gS[nANN+h]->SetMarkerStyle(h+1);
946  if(gS[nANN+h]->GetN()!=0) mg2->Add(gS[nANN+h]);
947  // gS[nANN+h]->Draw("apl,same");
948  }
949  mg2->Draw("apl");
950  //cS->cd(3)->SetLogy();
951  cS->cd(3);
952  TMultiGraph* mg3=new TMultiGraph();
953 // gS[nANN*2]->SetMarkerStyle(7);
954  gS[nANN*2]->SetMarkerColor(2);
955  gS[nANN*2]->SetLineColor(2);
956  mg3->SetTitle("cc=0.6;rho;#Events");
957  if(gS[nANN*2]->GetN()!=0) mg3->Add(gS[nANN*2]);
958  for (int h=1;h<nANN;h++){
959  gS[2*nANN+h]->SetMarkerColor(3);
960  gS[2*nANN+h]->SetLineColor(3);
961 // gS[2*nANN+h]->SetMarkerStyle(h+1);
962  if(gS[2*nANN+h]->GetN()!=0) mg3->Add(gS[2*nANN+h]);
963  // gS[2*nANN+h]->Draw("apl,same");
964  }
965  mg3->Draw("apl");
966  //cS->cd(4)->SetLogy();
967  cS->cd(4);
968  TMultiGraph* mg4=new TMultiGraph();
969 // gS[nANN*3]->SetMarkerStyle(7);
970  gS[nANN*3]->SetMarkerColor(2);
971  gS[nANN*3]->SetLineColor(2);
972  mg4->SetTitle("cc=0.65;rho;#Events");
973  if(gS[nANN*3]->GetN()!=0) mg4->Add(gS[nANN*3]);
974  // gS[nANN*3]->Draw("apl");
975  for (int h=1;h<nANN;h++){
976  gS[3*nANN+h]->SetMarkerColor(3);
977  gS[3*nANN+h]->SetLineColor(3);
978 // gS[3*nANN+h]->SetMarkerStyle(h+1);
979  if(gS[3*nANN+h]->GetN()!=0) mg4->Add(gS[3*nANN+h]);
980  // gS[3*nANN+h]->Draw("apl,same");
981  }
982  mg4->Draw("apl");
983  cout<<"nuovo canv"<<endl;
984  TCanvas* cB=new TCanvas("Number_vs_rho","Number_vs_rho",0,0,1200,700);
985  cB->Divide(2,2);
986  cB->cd(1)->SetLogy();
987  TMultiGraph* mg1B=new TMultiGraph();
988 // gB[0]->SetMarkerStyle(7);
989  gB[0]->SetMarkerColor(2);
990  gB[0]->SetLineColor(2);
991  mg1B->SetTitle("cc=0.5;rho;#Events");
992  if(gB[0]->GetN()!=0) mg1B->Add(gB[0]);
993  for (int h=1;h<nANN;h++){
994  gB[h]->SetMarkerColor(3);
995  gB[h]->SetLineColor(3);
996  // gB[h]>SetMarkerStyle(h+1);
997  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
998  //gB[h]->Draw("apl,same");
999  }
1000  mg1B->Draw("apl");
1001  cB->cd(2)->SetLogy();
1002  TMultiGraph* mg2B=new TMultiGraph();
1003  //gB[nANN]->SetMarkerStyle(7);
1004  gB[nANN]->SetMarkerColor(2);
1005  gB[nANN]->SetLineColor(2);
1006  mg2B->SetTitle("cc=0.55;rho;#Events");
1007  if(gB[nANN]->GetN()!=0) mg2B->Add(gB[nANN]);
1008  for (int h=1;h<nANN;h++){
1009  gB[nANN+h]->SetMarkerColor(3);
1010  gB[nANN+h]->SetLineColor(3);
1011  //gB[nANN+h]>SetMarkerStyle(h+1);
1012  if(gB[nANN+h]->GetN()!=0) mg2B->Add(gB[nANN+h]);
1013  //gB[h]->Draw("apl,same");
1014  }
1015  mg2B->Draw("apl");
1016  cB->cd(3)->SetLogy();
1017  TMultiGraph* mg3B=new TMultiGraph();
1018 // gB[2*nANN]->SetMarkerStyle(7);
1019  gB[2*nANN]->SetMarkerColor(2);
1020  gB[2*nANN]->SetLineColor(2);
1021  mg3B->SetTitle("cc=0.6;rho;#Events");
1022  if(gB[2*nANN]->GetN()!=0) mg3B->Add(gB[2*nANN]);
1023  for (int h=1;h<nANN;h++){
1024  gB[2*nANN+h]->SetMarkerColor(3);
1025  gB[2*nANN+h]->SetLineColor(3);
1026 // gB[2*nANN+h]>SetMarkerStyle(h+1);
1027  if(gB[2*nANN+h]->GetN()!=0) mg3B->Add(gB[2*nANN+h]);
1028  //gB[h]->Draw("apl,same");
1029  }
1030  mg3B->Draw("apl");
1031  cB->cd(4)->SetLogy();
1032  TMultiGraph* mg4B=new TMultiGraph();
1033 // gB[3*nANN]->SetMarkerStyle(7);
1034  gB[3*nANN]->SetMarkerColor(2);
1035  gB[3*nANN]->SetLineColor(2);
1036  mg4B->SetTitle("cc=0.65;rho;#Events");
1037  if(gB[3*nANN]->GetN()!=0) mg4B->Add(gB[3*nANN]);
1038  for (int h=1;h<nANN;h++){
1039  gB[3*nANN+h]->SetMarkerColor(3);
1040  gB[3*nANN+h]->SetLineColor(3);
1041 // gB[3*nANN+h]>SetMarkerStyle(h+1);
1042  if(gB[3*nANN+h]->GetN()!=0) mg4B->Add(gB[3*nANN+h]);
1043  //gB[h]->Draw("apl,same");
1044  }
1045  mg4B->Draw("apl");
1046 
1047 // cout<<"dopo Draw()"<<endl;
1048  TString CnameS(name);
1049  TString CnameB(name);
1050  TString CnameSroot(name);
1051  TString CnameBroot(name);
1052  char CnameS2[1024];
1053  char CnameB2[1024];
1054  char CnameS2root[1024];
1055  char CnameB2root[1024];
1056  CnameS.ReplaceAll(".root",".png");
1057  CnameB.ReplaceAll(".root",".png");
1058  sprintf(CnameS2,"logN_rho_av/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameS.Data());
1059  sprintf(CnameB2,"logN_rho_av/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameB.Data());
1060  sprintf(CnameS2root,"logN_rho_av/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameSroot.Data());
1061  sprintf(CnameB2root,"logN_rho_av/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameBroot.Data());
1062  cS->SaveAs(CnameS2);
1063  cB->SaveAs(CnameB2);
1064  cS->Print(CnameS2root);
1065  cB->Print(CnameB2root);
1066 // cout<<"fine"<<endl;
1067 
1068 }
1069 
1070 
1071 void Annth(TString ifile, int &FAth_n, int &FDth_n, double &thresFA){
1072  TString name(ifile);
1073  name.ReplaceAll("average_file/","");
1074  TFile* fTEST =TFile::Open(ifile.Data());
1075  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1076  double av;
1077  double cc;
1078  double rho;
1079  int nSi;
1080  NNTree2->SetBranchAddress("Average",&av);
1081  NNTree2->SetBranchAddress("cc",&cc);
1082  NNTree2->SetBranchAddress("rho",&rho);
1083  NNTree2->SetBranchAddress("#TestSig",&nSi);
1084  NNTree2->GetEntry(0);
1085  int const nSig=nSi;
1086  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1087  int const ncurve2=nRHO*nCC;
1088 
1089 
1090  double* ANNSig[ncurve2];
1091  for (int i=0;i<ncurve2;i++) ANNSig[i]=new double[nSig];
1092  //double rhoSig[20][10];
1093  int NSig[ncurve2];
1094 // int nBg=0;
1095  int const nBg=NNTree2->GetEntries()-nSig;
1096  for (int i=0;i<ncurve2;i++) {
1097  NSig[i]=0;
1098  for (int j=0;j<nSig;j++) ANNSig[i][j]=0.;
1099  }
1100  double* ANNBg[ncurve2];
1101  for (int i=0;i<ncurve2;i++) ANNBg[i]=new double[nBg];
1102 
1103  //double rhoBg[20][10];
1104  int NBg[ncurve2];
1105  for (int i=0;i<ncurve2;i++) {
1106  NBg[i]=0;
1107  for (int j=0;j<nBg;j++) ANNBg[i][j]=0.;
1108  }
1109  double ccTh[nCC];
1110  for (int i=0;i<nCC;i++) ccTh[i]=0.;
1111  double rhoTh[nRHO];
1112  for (int i=0;i<nRHO;i++) rhoTh[i]=0.;
1113  // double deltacc=0.;
1114  // double deltarho=0.;
1115  // deltacc=0.2/nCC;
1116 // deltarho=1./(nRHO);
1117 
1118 //int thresFA=-1;
1119  for(int n=0;n<NNTree2->GetEntries();n++){
1120  NNTree2->GetEntry(n);
1121  cout<<"rho "<<rho<<" cc "<<cc<<" av "<<av<<endl;
1122  for(int i=0; i<nCC;i++){
1123  ccTh[i]=0.5+i*deltacc;
1124  ccTh[0]=0.;
1125  if(cc<ccTh[i]) continue;
1126  for(int j=0; j<nRHO;j++){
1127  rhoTh[j]=5+j*deltarho;
1128  rhoTh[0]=0.;
1129  if(rho<rhoTh[j]) continue;
1130  int ni=0;
1131  if(n>=nSig) {
1132  NBg[i*nRHO+j]= NBg[i*nRHO+j]+1;
1133  if(i*nRHO+j==0) cout<<" NBg[i*nRHO+j] "<<NBg[i*nRHO+j]<<endl;
1134  while(ANNBg[i*nRHO+j][ni]!=0)ni=ni+1;
1135  ANNBg[i*nRHO+j][ni]=av;
1136  // if (NBg[0]==FAth_n) thresFA=av;
1137  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1138  }
1139  else {
1140  NSig[i*nRHO+j]= NSig[i*nRHO+j]+1;
1141  if(i*nRHO+j==0) cout<<" NSig[i*nRHO+j] "<<NSig[i*nRHO+j]<<endl;
1142  while(ANNSig[i*nRHO+j][ni]!=0)ni=ni+1;
1143  ANNSig[i*nRHO+j][ni]=av;
1144  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1145  }
1146  }
1147  }
1148  }
1149  int* indexS[ncurve2];
1150  for (int i=0;i<ncurve2;i++) indexS[i]=new int[nSig];
1151 
1152  TGraph * gS[ncurve2];
1153  for (int y=0;y<ncurve2;y++) {
1154  int igS=0;
1155  int igS_p=0;
1156  gS[y]=new TGraph();
1157  TMath::Sort(nSig,ANNSig[y],indexS[y],false);
1158  for (int k=0;k<nSig;k++) {
1159  int ii=indexS[y][k];
1160  int yy=0;
1161  if (k>0){
1162  //cout<<"k "<<k<<endl;
1163  int ij=indexS[y][k-1];
1164  //if(ANNSig[y][ii]!=0&&ANNSig[y][ii]!=ANNSig[y][ij]) {
1165  if(ANNSig[y][ii]!=0) {
1166  yy=NSig[y]-igS;
1167  //gS[y]->SetPoint(igS,ANNSig[y][ii],yy);
1168  if(ANNSig[y][ii]!=ANNSig[y][ij]) gS[y]->SetPoint(igS_p++,ANNSig[y][ii],yy);
1169  // cout<<"igS"<<igS<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1170  igS=igS+1;
1171 
1172  }
1173  }
1174  else {
1175  if(ANNSig[y][ii]!=0){
1176  yy=NSig[y]-igS;
1177  gS[y]->SetPoint(0,ANNSig[y][ii],yy);
1178  igS=igS+1;
1179  // cout<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1180  }
1181 
1182  }
1183  }
1184  }
1185 
1186 
1187 int thres_ind=0;
1188  int* indexB[ncurve2];
1189  for (int i=0;i<ncurve2;i++) indexB[i]=new int[nBg];
1190 
1191  TGraph * gB[ncurve2];
1192  for (int y=0;y<ncurve2;y++) {
1193  int igB=0;
1194  int igB_p=0;
1195  gB[y]=new TGraph();
1196  TMath::Sort(nBg,ANNBg[y],indexB[y],false);
1197  // cout<<"dopo Sort "<<y<<endl;
1198  for (int k=0;k<nBg;k++) {
1199  int ii=indexB[y][k];
1200  int ik=indexB[y][k-1];
1201  int yy=0;
1202  if (k>0){
1203  int ij=indexB[y][k-1];
1204  //if(ANNBg[y][ii]!=0&&ANNBg[y][ii]!=ANNBg[y][ij]) {
1205  if(ANNBg[y][ii]!=0) {
1206  yy=NBg[y]-igB;
1207  if(y==0 && yy==FAth_n+1) {thresFA=ANNBg[y][ii];}
1208  if (y==0 && yy<FAth_n+1) if(ANNBg[y][ii]==ANNBg[y][ik]){FAth_n=FAth_n-1;}
1209  if(y==0) cout<<"yy "<<yy<<" y "<<y<<" FAth_n "<<FAth_n<<" thresFA "<<thresFA<<" ANNBg[y][ii] "<<ANNBg[y][ii]<<endl;
1210  //gB[y]->SetPoint(igB,ANNBg[y][ii],yy);
1211  if(ANNBg[y][ii]!=ANNBg[y][ij]) gB[y]->SetPoint(igB_p++,ANNBg[y][ii],yy);
1212  igB=igB+1;
1213  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1214  }
1215  }
1216  else {
1217  if(ANNBg[y][ii]!=0) {
1218  yy=NBg[y]-igB;
1219  gB[y]->SetPoint(0,ANNBg[y][ii],yy);
1220  igB=igB+1;
1221  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1222  }
1223  }
1224  }
1225  }
1226 //int FDth_n=0;
1227 //for(int i=0; i<nCC;i++){
1228 // for(int j=0; j<nRHO;j++){
1229 FDth_n=0;
1230  for(int ni=1;ni<nSig;ni++){
1231  cout<<"FAth_n "<<FAth_n<<" thresFA "<<thresFA<<" ANNSig[0][indexS[0][ni]] "<<ANNSig[0][indexS[0][ni]]<<" FDth_n "<<FDth_n<<" NSig[0]-(ni-1 )"<<NSig[0]-(ni-1)<<" NSig[0] "<<NSig[0]<<endl;
1232 // if (ANNSig[0][indexS[0][ni]]==thresFA || (ANNSig[0][indexS[0][ni-1]]<thresFA && ANNSig[0][indexS[0][ni]]>thresFA)) FDth_n=NSig[0]-(ni-1);
1233  if (ANNSig[0][indexS[0][ni]]<=thresFA) FDth_n=FDth_n+1;
1234  }
1235 // }
1236 //}
1237 //exit(0);
1238  TCanvas* cS=new TCanvas("Efficiency_vs_ANN","Efficiency_vs_ANN",0,0,1200,700);
1239  cS->Divide(2,2);
1240  //cS->cd(1)->SetLogy();
1241  cS->cd(1);
1242  TMultiGraph* mg1=new TMultiGraph();
1243  mg1->SetTitle("cc: no_cut;ANN;#Events");
1244  for (int h=0;h<nRHO;h++){
1245  gS[h]->SetLineColor(4);
1246  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
1247  }
1248  mg1->Draw("al");
1249  cS->cd(2);
1250  // cS->cd(2)->SetLogy();
1251  TMultiGraph* mg2=new TMultiGraph();
1252  mg2->SetTitle("cc=0.55;ANN;#Events");
1253  for (int h=0;h<nRHO;h++){
1254  gS[nRHO+h]->SetLineColor(4);
1255  if(gS[nRHO+h]->GetN()!=0) mg2->Add(gS[nRHO+h]);
1256  }
1257  mg2->Draw("al");
1258 
1259  //cS->cd(3)->SetLogy();
1260  cS->cd(3);
1261  TMultiGraph* mg3=new TMultiGraph();
1262  mg3->SetTitle("cc=0.6;ANN;#Events");
1263  for (int h=0;h<nRHO;h++){
1264  gS[2*nRHO+h]->SetLineColor(4);
1265  if(gS[2*nRHO+h]->GetN()!=0) mg3->Add(gS[2*nRHO+h]);
1266  }
1267  mg3->Draw("al");
1268  //cS->cd(4)->SetLogy();
1269  cS->cd(4);
1270  TMultiGraph* mg4=new TMultiGraph();
1271  mg4->SetTitle("cc=0.65;ANN;#Events");
1272  for (int h=0;h<nRHO;h++){
1273  gS[3*nRHO+h]->SetLineColor(4);
1274  if(gS[3*nRHO+h]->GetN()!=0) mg4->Add(gS[3*nRHO+h]);
1275  }
1276  mg4->Draw("al");
1277 
1278  TCanvas* cB=new TCanvas("Number_vs_ANN","Number_vs_ANN",0,0,1200,700);
1279  cB->Divide(2,2);
1280  cB->cd(1)->SetLogy();
1281  TMultiGraph* mg1B=new TMultiGraph();
1282  mg1B->SetTitle("cc: no_cut;ANN;#Events");
1283  for (int h=0;h<nRHO;h++){
1284  gB[h]->SetLineColor(4);
1285  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
1286  }
1287  mg1B->Draw("al");
1288  cB->cd(2)->SetLogy();
1289  TMultiGraph* mg2B=new TMultiGraph();
1290  mg2B->SetTitle("cc=0.55;ANN;#Events");
1291  for (int h=0;h<nRHO;h++){
1292  gB[nRHO+h]->SetLineColor(4);
1293  if(gB[nRHO+h]->GetN()!=0) mg2B->Add(gB[nRHO+h]);
1294  }
1295  mg2B->Draw("al");
1296  cB->cd(3)->SetLogy();
1297  TMultiGraph* mg3B=new TMultiGraph();
1298  mg3B->SetTitle("cc=0.6;ANN;#Events");
1299  for (int h=0;h<nRHO;h++){
1300  gB[2*nRHO+h]->SetLineColor(4);
1301  if(gB[2*nRHO+h]->GetN()!=0) mg3B->Add(gB[2*nRHO+h]);
1302  }
1303  mg3B->Draw("al");
1304  cB->cd(4)->SetLogy();
1305  TMultiGraph* mg4B=new TMultiGraph();
1306  mg4B->SetTitle("cc=0.65;ANN;#Events");
1307  for (int h=0;h<nRHO;h++){
1308  gB[3*nRHO+h]->SetLineColor(4);
1309  if(gB[3*nRHO+h]->GetN()!=0) mg4B->Add(gB[3*nRHO+h]);
1310  }
1311  mg4B->Draw("al");
1312 
1313 
1314  //cout<<"dopo Draw()"<<endl;
1315  TString CnameS(name);
1316  TString CnameB(name);
1317  TString CnameSroot(name);
1318  TString CnameBroot(name);
1319  char CnameS2[1024];
1320  char CnameB2[1024];
1321  char CnameS2root[1024];
1322  char CnameB2root[1024];
1323  CnameS.ReplaceAll(".root",".png");
1324  CnameB.ReplaceAll(".root",".png");
1325  sprintf(CnameS2,"ANNthres_av/N_ANN_S_%s",CnameS.Data());
1326  sprintf(CnameB2,"ANNthres_av/N_ANN_B_%s",CnameB.Data());
1327  sprintf(CnameS2root,"ANNthres_av/N_ANN_S_%s",CnameSroot.Data());
1328  sprintf(CnameB2root,"ANNthres_av/N_ANN_B_%s",CnameBroot.Data());
1329  cS->SaveAs(CnameS2);
1330  cB->SaveAs(CnameB2);
1331  cS->Print(CnameS2root);
1332  cB->Print(CnameB2root);
1333  //cout<<"fine"<<endl;
1334 
1335 //CARTELLA Annth
1336 
1337 
1338 
1339 }
1340 
1341 void Mcth(TString ifile, int &McFAth_n, int &McFDth_n, double &McthresFA){
1342 TString name(ifile);
1343  name.ReplaceAll("outfile/","");
1344  name.ReplaceAll("average_file/","");
1345  TFile* fTEST =TFile::Open(ifile.Data());
1346  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1347  //double outbis;
1348  double av;
1349  double out;
1350  double cc;
1351  double Mc;
1352  double Mc_inj;
1353  double rho;
1354  int nSi;
1355 // NNTree2->SetBranchAddress("ANNout",&out);
1356  NNTree2->SetBranchAddress("Mchirp",&Mc);
1357  NNTree2->SetBranchAddress("Average",&av);
1358  NNTree2->SetBranchAddress("cc",&cc);
1359  NNTree2->SetBranchAddress("rho",&rho);
1360  NNTree2->SetBranchAddress("Mchirp_injected",&Mc_inj);
1361  NNTree2->SetBranchAddress("#TestSig",&nSi);
1362  NNTree2->GetEntry(0);
1363  int const nSig=nSi;
1364  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1365  int const ncurve2=nRHO*nCC;
1366 
1367  double* McSig[ncurve2];
1368  for (int i=0;i<ncurve2;i++) McSig[i]=new double[nSig];
1369  //double rhoSig[20][10];
1370  int NSig[ncurve2];
1371 // int nBg=0;
1372  int const nBg=NNTree2->GetEntries()-nSig;
1373  for (int i=0;i<ncurve2;i++) {
1374  NSig[i]=0;
1375  for (int j=0;j<nSig;j++) McSig[i][j]=0.;
1376  }
1377  double* McBg[ncurve2];
1378  for (int i=0;i<ncurve2;i++) McBg[i]=new double[nBg];
1379 
1380  //double rhoBg[20][10];
1381  int NBg[ncurve2];
1382  for (int i=0;i<ncurve2;i++) {
1383  NBg[i]=0;
1384  for (int j=0;j<nBg;j++) McBg[i][j]=0.;
1385  }
1386  double ccTh[nCC];
1387  for (int i=0;i<nCC;i++) ccTh[i]=0.;
1388  double rhoTh[nRHO];
1389  for (int i=0;i<nRHO;i++) rhoTh[i]=0.;
1390 // double deltacc=0.;
1391 // double deltarho=0.;
1392 // deltacc=0.2/nCC;
1393 // deltarho=1./(nRHO);
1394 
1395  for(int n=0;n<NNTree2->GetEntries();n++){
1396  NNTree2->GetEntry(n);
1397  cout<<"rho "<<rho<<" cc "<<cc<<" av "<<av<<endl;
1398  for(int i=0; i<nCC;i++){
1399  ccTh[i]=0.5+i*deltacc;
1400  ccTh[0]=0.;
1401  if(cc<ccTh[i]) continue;
1402  for(int j=0; j<nRHO;j++){
1403  rhoTh[j]=5+j*deltarho;
1404  rhoTh[0]=0.;
1405  if(rho<rhoTh[j]) continue;
1406  int ni=0;
1407  if(n>=nSig) {
1408  NBg[i*nRHO+j]= NBg[i*nRHO+j]+1;
1409  if(i*nRHO+j==0) cout<<" NBg[i*nRHO+j] "<<NBg[i*nRHO+j]<<endl;
1410  while(McBg[i*nRHO+j][ni]!=0)ni=ni+1;
1411  McBg[i*nRHO+j][ni]=Mc;
1412  // if (NBg[0]==FAth_n) thresFA=av;
1413  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1414  }
1415  else {
1416  NSig[i*nRHO+j]= NSig[i*nRHO+j]+1;
1417  if(i*nRHO+j==0) cout<<" NSig[i*nRHO+j] "<<NSig[i*nRHO+j]<<endl;
1418  while(McSig[i*nRHO+j][ni]!=0)ni=ni+1;
1419  McSig[i*nRHO+j][ni]=Mc;
1420  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1421  }
1422  }
1423  }
1424  }
1425 
1426  int* indexS[ncurve2];
1427  for (int i=0;i<ncurve2;i++) indexS[i]=new int[nSig];
1428 
1429  TGraph * gS[ncurve2];
1430  for (int y=0;y<ncurve2;y++) {
1431  int igS=0;
1432  int igS_p=0;
1433  gS[y]=new TGraph();
1434  TMath::Sort(nSig,McSig[y],indexS[y],false);
1435  for (int k=0;k<nSig;k++) {
1436  int ii=indexS[y][k];
1437  int yy=0;
1438  if (k>0){
1439  //cout<<"k "<<k<<endl;
1440  int ij=indexS[y][k-1];
1441  //if(ANNSig[y][ii]!=0&&ANNSig[y][ii]!=ANNSig[y][ij]) {
1442  if(McSig[y][ii]!=0) {
1443  yy=NSig[y]-igS;
1444  //gS[y]->SetPoint(igS,ANNSig[y][ii],yy);
1445  if(McSig[y][ii]!=McSig[y][ij]) gS[y]->SetPoint(igS_p++,McSig[y][ii],yy);
1446  // cout<<"igS"<<igS<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1447  igS=igS+1;
1448 
1449  }
1450  }
1451  else {
1452  if(McSig[y][ii]!=0){
1453  yy=NSig[y]-igS;
1454  gS[y]->SetPoint(0,McSig[y][ii],yy);
1455  igS=igS+1;
1456  // cout<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1457  }
1458 
1459  }
1460  }
1461  }
1462 
1463 int thres_ind=0;
1464  int* indexB[ncurve2];
1465  for (int i=0;i<ncurve2;i++) indexB[i]=new int[nBg];
1466 
1467  TGraph * gB[ncurve2];
1468  for (int y=0;y<ncurve2;y++) {
1469  int igB=0;
1470  int igB_p=0;
1471  gB[y]=new TGraph();
1472  TMath::Sort(nBg,McBg[y],indexB[y],false);
1473  // cout<<"dopo Sort "<<y<<endl;
1474  for (int k=0;k<nBg;k++) {
1475  int ii=indexB[y][k];
1476  int ik=indexB[y][k-1];
1477  int yy=0;
1478  if (k>0){
1479  int ij=indexB[y][k-1];
1480  //if(ANNBg[y][ii]!=0&&ANNBg[y][ii]!=ANNBg[y][ij]) {
1481  if(McBg[y][ii]!=0) {
1482  yy=NBg[y]-igB;
1483  if(y==0 && yy==McFAth_n+1) {McthresFA=McBg[y][ii];
1484  }
1485  if (y==0 && yy<McFAth_n+1) { cout<<"dentro primo if"<<endl;
1486  //if(McBg[y][ii]==McBg[y][ik]){McFAth_n=McFAth_n-1; cout<<"dentro secondo if"<<endl;}
1487  if(McBg[y][ii]>McBg[y][ik]){McFAth_n=McFAth_n;cout<<"maggiore ii "<<" diff"<<(McBg[y][ii]-McBg[y][ik])<<endl;}
1488  if(McBg[y][ii]<McBg[y][ik]){McFAth_n=McFAth_n;cout<<"minore ii"<<endl;}
1489  if(McBg[y][ii]==McBg[y][ik]){McFAth_n=McFAth_n-1;cout<<"uguale"<<endl;}
1490  //else { McFAth_n=McFAth_n-1;cout<<"dentro secondo if"<<endl;}
1491 }
1492  if(y==0) cout<<"ii "<<ii<<" yy "<<yy<<" y "<<y<<" McFAth_n "<<McFAth_n<<"McthresFA "<<McthresFA<<" McBg[y][ii] "<<McBg[y][ii]<<" McBg[y][indexB[y][k-1] "<<McBg[y][indexB[y][k-1]]<<endl;
1493  //gB[y]->SetPoint(igB,ANNBg[y][ii],yy);
1494  if(McBg[y][ii]!=McBg[y][ij]) gB[y]->SetPoint(igB_p++,McBg[y][ii],yy);
1495  igB=igB+1;
1496  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1497  }
1498  }
1499  else {
1500  if(McBg[y][ii]!=0) {
1501  yy=NBg[y]-igB;
1502  gB[y]->SetPoint(0,McBg[y][ii],yy);
1503  igB=igB+1;
1504  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1505  }
1506  }
1507  }
1508  }
1509 
1510 McFDth_n=0;
1511  for(int ni=1;ni<nSig;ni++){
1512  cout<<"McFAth_n "<<McFAth_n<<" McthresFA "<<McthresFA<<" McSig[0][indexS[0][ni]] "<<McSig[0][indexS[0][ni]]<<" McFDth_n "<<McFDth_n<<" NSig[0]-(ni-1 )"<<NSig[0]-(ni-1)<<" NSig[0] "<<NSig[0]<<endl;
1513 // if (ANNSig[0][indexS[0][ni]]==thresFA || (ANNSig[0][indexS[0][ni-1]]<thresFA && ANNSig[0][indexS[0][ni]]>thresFA)) FDth_n=NSig[0]-(ni-1);
1514  if (McSig[0][indexS[0][ni]]<=McthresFA) McFDth_n=McFDth_n+1;
1515  }
1516 
1517 
1518  TCanvas* cS=new TCanvas("Efficiency_vs_Mchirp","Efficiency_vs_Mchirp",0,0,1200,700);
1519  cS->Divide(2,2);
1520  //cS->cd(1)->SetLogy();
1521  cS->cd(1);
1522  TMultiGraph* mg1=new TMultiGraph();
1523  mg1->SetTitle("cc: no_cut;Mchirp;#Events");
1524  for (int h=0;h<nRHO;h++){
1525  gS[h]->SetLineColor(4);
1526  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
1527  }
1528  mg1->Draw("al");
1529  cS->cd(2);
1530  // cS->cd(2)->SetLogy();
1531  TMultiGraph* mg2=new TMultiGraph();
1532  mg2->SetTitle("cc=0.55;Mchirp;#Events");
1533  for (int h=0;h<nRHO;h++){
1534  gS[nRHO+h]->SetLineColor(4);
1535  if(gS[nRHO+h]->GetN()!=0) mg2->Add(gS[nRHO+h]);
1536  }
1537  mg2->Draw("al");
1538 
1539  //cS->cd(3)->SetLogy();
1540  cS->cd(3);
1541  TMultiGraph* mg3=new TMultiGraph();
1542  mg3->SetTitle("cc=0.6;Mchirp;#Events");
1543  for (int h=0;h<nRHO;h++){
1544  gS[2*nRHO+h]->SetLineColor(4);
1545  if(gS[2*nRHO+h]->GetN()!=0) mg3->Add(gS[2*nRHO+h]);
1546  }
1547  mg3->Draw("al");
1548  //cS->cd(4)->SetLogy();
1549  cS->cd(4);
1550  TMultiGraph* mg4=new TMultiGraph();
1551  mg4->SetTitle("cc=0.65;Mchirp;#Events");
1552  for (int h=0;h<nRHO;h++){
1553  gS[3*nRHO+h]->SetLineColor(4);
1554  if(gS[3*nRHO+h]->GetN()!=0) mg4->Add(gS[3*nRHO+h]);
1555  }
1556  mg4->Draw("al");
1557 
1558  TCanvas* cB=new TCanvas("Number_vs_Mchirp","Number_vs_Mchirp",0,0,1200,700);
1559  cB->Divide(2,2);
1560  cB->cd(1)->SetLogy();
1561  TMultiGraph* mg1B=new TMultiGraph();
1562  mg1B->SetTitle("cc: no_cut;Mchirp;#Events");
1563  for (int h=0;h<nRHO;h++){
1564  gB[h]->SetLineColor(4);
1565  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
1566  }
1567  mg1B->Draw("al");
1568  cB->cd(2)->SetLogy();
1569  TMultiGraph* mg2B=new TMultiGraph();
1570  mg2B->SetTitle("cc=0.55;Mchirp;#Events");
1571  for (int h=0;h<nRHO;h++){
1572  gB[nRHO+h]->SetLineColor(4);
1573  if(gB[nRHO+h]->GetN()!=0) mg2B->Add(gB[nRHO+h]);
1574  }
1575  mg2B->Draw("al");
1576  cB->cd(3)->SetLogy();
1577  TMultiGraph* mg3B=new TMultiGraph();
1578  mg3B->SetTitle("cc=0.6;Mchirp;#Events");
1579  for (int h=0;h<nRHO;h++){
1580  gB[2*nRHO+h]->SetLineColor(4);
1581  if(gB[2*nRHO+h]->GetN()!=0) mg3B->Add(gB[2*nRHO+h]);
1582  }
1583  mg3B->Draw("al");
1584  cB->cd(4)->SetLogy();
1585  TMultiGraph* mg4B=new TMultiGraph();
1586  mg4B->SetTitle("cc=0.6;Mchirp;#Events");
1587  for (int h=0;h<nRHO;h++){
1588  gB[3*nRHO+h]->SetLineColor(4);
1589  if(gB[3*nRHO+h]->GetN()!=0) mg4B->Add(gB[3*nRHO+h]);
1590  }
1591  mg4B->Draw("al");
1592 
1593 
1594  TString CnameS(name);
1595  TString CnameB(name);
1596  TString CnameSroot(name);
1597  TString CnameBroot(name);
1598  char CnameS2[1024];
1599  char CnameB2[1024];
1600  char CnameS2root[1024];
1601  char CnameB2root[1024];
1602  CnameS.ReplaceAll(".root",".png");
1603  CnameB.ReplaceAll(".root",".png");
1604  sprintf(CnameS2,"Mcthres/N_Mc_S_%s",CnameS.Data());
1605  sprintf(CnameB2,"Mcthres/N_Mc_B_%s",CnameB.Data());
1606  sprintf(CnameS2root,"Mcthres/N_Mc_S_%s",CnameSroot.Data());
1607  sprintf(CnameB2root,"Mcthres/N_Mc_B_%s",CnameBroot.Data());
1608  cS->SaveAs(CnameS2);
1609  cB->SaveAs(CnameB2);
1610  cS->Print(CnameS2root);
1611  cB->Print(CnameB2root);
1612 
1613 
1614 }
1615 
1616 void PlotsAv_cc(TString ifile){
1617  TString name(ifile);
1618  name.ReplaceAll("outfile/","");
1619  name.ReplaceAll("average_file/","");
1620  TFile* fTEST =TFile::Open(ifile.Data());
1621  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1622  //double outbis;
1623  double av;
1624  double out;
1625  double cc;
1626  double rho;
1627  int nSi;
1628 // NNTree2->SetBranchAddress("ANNout",&out);
1629  NNTree2->SetBranchAddress("Average",&av);
1630  NNTree2->SetBranchAddress("cc",&cc);
1631  NNTree2->SetBranchAddress("rho",&rho);
1632  NNTree2->SetBranchAddress("#TestSig",&nSi);
1633  NNTree2->GetEntry(0);
1634  int const nSig=nSi;
1635  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1636  int const nBg=NNTree2->GetEntries()-nSig;
1637  cout<<"nBg: "<<nBg<<" Entries: "<<NNTree2->GetEntries()<<endl;
1638 
1639  TGraph * gS[3];
1640  TGraph * gB[3];
1641  gB[0]=new TGraph();
1642  gS[0]=new TGraph();
1643 // cout<<"nSig: "<<nSig<<endl;
1644  for (int n = 0; n <NNTree2->GetEntries(); n++){
1645  NNTree2->GetEntry(n);
1646  if(n<nSig) {
1647  gS[0]->SetPoint(n,cc,av);
1648  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1649  }
1650  else {
1651  gB[0]->SetPoint(n-nSig,cc,av);
1652  cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1653  }
1654  }
1655 
1656 
1657  /*for (int a=0;a<nBg;a++){
1658  if(aRHOB[a]>=RHOth){
1659  for (int b=0;b<nth;b++){
1660  if(aCCB[a]>=cc1th[b]){
1661  for(int c=0;c<nth+1;c++){
1662  if(aANNB[a]>=ANN1th[c]) Z[b*nth+c]=Z[b*nth+c]+1;
1663  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1664  }
1665  }
1666  }
1667  }
1668  }*/
1669  //gS[0]->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1670  //gB[0]->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1671  //gS[0]->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1672  //gB[0]->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1673  gS[0]->SetMarkerColor(2);
1674  gB[0]->SetMarkerColor(4);
1675  gS[0]->SetMarkerStyle(6);
1676  gB[0]->SetMarkerStyle(7);
1677 
1678  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1679  c->Divide(1,2);
1680  c->cd(1);
1681  TMultiGraph* mg1=new TMultiGraph();
1682  mg1->SetTitle("Av_cc");
1683  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1684  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1685  mg1->Draw("ap");
1686  mg1->GetHistogram()->GetXaxis()->SetTitle("cc");
1687  mg1->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1688 // mg1->Draw("ap");
1689  c->cd(2);
1690  TMultiGraph* mg2=new TMultiGraph();
1691  mg2->SetTitle("Av_cc");
1692  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1693  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1694  mg2->Draw("ap");
1695  mg2->GetHistogram()->GetXaxis()->SetTitle("cc");
1696  mg2->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1697 
1698  cout<<" name "<<name<<endl;
1699  TString Cname(name);
1700  TString Cnameroot(name);
1701  char Cname2[1024];
1702  char Cname2root[1024];
1703  Cname.ReplaceAll(".root",".png");
1704  sprintf(Cname2,"average_png/out_Plots_%s",Cname.Data());
1705  sprintf(Cname2root,"average_png/out_Plots_%s",Cnameroot.Data());
1706  c->SaveAs(Cname2);
1707  c->Print(Cname2root);
1708 /*
1709  TCanvas* c2=new TCanvas("Plots_Bkg_on_Sig","Plots_Bkg_on_Sig",0,0,1200,700);
1710  c2->Divide(2,2);
1711  c2->cd(1);
1712  gS[0]->Draw("ap");
1713  gB[0]->Draw("p,same");
1714 
1715  c2->cd(2);
1716  gS[1]->Draw("ap");
1717  gB[1]->Draw("p,same");
1718  c2->cd(3);
1719  gS[2]->Draw("ap");
1720  gB[2]->Draw("p,same");
1721  c2->cd(4);
1722  TText* text2=new TText(0.37,0.0,"no cuts on ANN");
1723  hglitch->SetStats(0);
1724  hglitch->GetXaxis()->SetTitle("cc");
1725  hglitch->GetYaxis()->SetTitle("ANNoutput");
1726  hglitch->GetZaxis()->SetTitle("count");
1727  hglitch->Draw("colz");
1728  text2->Draw();
1729  gPad->SetLogz();
1730 
1731  TString Cname_2(name);
1732  TString Cname_2root(name);
1733  char Cname2_2[1024];
1734  char Cname2_2root[1024];
1735  Cname_2.ReplaceAll(".root",".png");
1736  sprintf(Cname2_2,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2.Data());
1737  sprintf(Cname2_2root,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2root.Data());
1738  c2->SaveAs(Cname2_2);
1739  c2->Print(Cname2_2root);
1740 
1741 //CARTELLA CCi_RHO_ANN_Plots */
1742 }
1743 void PlotsAv_Mc(TString ifile){
1744  TString name(ifile);
1745  name.ReplaceAll("outfile/","");
1746  name.ReplaceAll("average_file/","");
1747  TFile* fTEST =TFile::Open(ifile.Data());
1748  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1749  //double outbis;
1750  double av;
1751  double out;
1752  double cc;
1753  double Mc;
1754  double Mc_inj;
1755  double rho;
1756  int nSi;
1757 // NNTree2->SetBranchAddress("ANNout",&out);
1758  NNTree2->SetBranchAddress("Mchirp",&Mc);
1759  NNTree2->SetBranchAddress("Average",&av);
1760  NNTree2->SetBranchAddress("cc",&cc);
1761  NNTree2->SetBranchAddress("rho",&rho);
1762  NNTree2->SetBranchAddress("Mchirp_injected",&Mc_inj);
1763  NNTree2->SetBranchAddress("#TestSig",&nSi);
1764  NNTree2->GetEntry(0);
1765  int const nSig=nSi;
1766  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1767  //exit(0);
1768  int const nBg=NNTree2->GetEntries()-nSig;
1769 
1770  TGraph * gS[3];
1771  TGraph * gB[3];
1772  gB[0]=new TGraph();
1773  gS[0]=new TGraph();
1774 
1775 
1776  TCanvas* ANN_Mc_c3D=new TCanvas("ANN_average_vs_Mc_reconstructed","ANN_average_vs_Mc_reconstructed",0,0,1200,700);
1777  ANN_Mc_c3D->Divide(1,2);
1778  TH2D* ANN_Mc_SIG3D=new TH2D("ANN_average_vs_Mc_recostructed_sig","ANN_average_vs_Mc_recostructed_sig",NMc_r,minMc_r,maxMc_r,NANN_av,minANN_av,maxANN_av);
1779  TH2D* ANN_Mc_BKG3D=new TH2D("ANN_average_vs_Mc_recostructed_bkg","ANN_average_vs_Mc_recostructed_bkg",NMc_r,minMc_r,maxMc_r,NANN_av,minANN_av,maxANN_av);
1780  // Mc_inj_g->SetPoint(0,rhoSig[y][ii],yy);
1781  ANN_Mc_BKG3D->SetDirectory(0);
1782  ANN_Mc_BKG3D->SetStats(0);
1783  ANN_Mc_SIG3D->SetDirectory(0);
1784  ANN_Mc_SIG3D->SetStats(0);
1785 
1786  TCanvas* ANN_cc_c3D=new TCanvas("ANN_average_vs_cc","ANN_average_vs_cc",0,0,1200,700);
1787  ANN_cc_c3D->Divide(1,2);
1788  TH2D* ANN_cc_SIG3D=new TH2D("ANN_average_vs_cc_SIG","ANN_average_vs_cc_SIG",Ncc3D, minCC, maxCC,NANN_av,minANN_av,maxANN_av);
1789  TH2D* ANN_cc_BKG3D=new TH2D("ANN_average_vs_cc_BKG","ANN_average_vs_cc_BKG",Ncc3D, minCC, maxCC,NANN_av,minANN_av,maxANN_av);
1790  ANN_cc_SIG3D->SetDirectory(0);
1791  ANN_cc_SIG3D->SetStats(0);
1792  ANN_cc_BKG3D->SetDirectory(0);
1793  ANN_cc_BKG3D->SetStats(0);
1794  double deltacc3D=0;
1795  deltacc3D=(maxCC-minCC)/(Ncc3D);
1796 
1797 
1798  TCanvas* Mc_inj_c=new TCanvas("ANN_vs_Mc_injected","ANN_vs_Mc_injected",0,0,1200,700);
1799  Mc_inj_c->Divide(1,2);
1800  TH2D* Mc_inj_gMc=new TH2D("Mc_injected_Mc_recostructed","Mc_injected_recostructed",NMc_inj,minMc_inj,maxMc_inj,NMc_r,minMc_r,maxMc_r);
1801  TH2D* Mc_inj_g=new TH2D("Mc_injected_comparison","Mc_injected_comparison",NMc_inj,minMc_inj,maxMc_inj,NANN_av,minANN_av,maxANN_av);
1802  // Mc_inj_g->SetPoint(0,rhoSig[y][ii],yy);
1803  Mc_inj_gMc->SetDirectory(0);
1804  Mc_inj_gMc->SetStats(0);
1805  Mc_inj_g->SetDirectory(0);
1806  Mc_inj_g->SetStats(0);
1807  double deltaMc_inj=0;
1808  deltaMc_inj=(maxMc_inj-minMc_inj)/(NMc_inj);
1809  double deltaANN_av_M=0;
1810  deltaANN_av_M=(maxANN_av-minANN_av)/(NANN_av);
1811  double deltaMc_r=0;
1812  deltaMc_r=(maxMc_r-minMc_r)/(NMc_r);
1813 
1814 
1815  TString txt_name0(name);
1816  txt_name0.ReplaceAll(".root",".txt");
1817  char txt_name[1024];
1818  sprintf(txt_name,"txt_files/%s",txt_name0.Data());
1819  ofstream file_txt_out(txt_name);
1820 
1821 
1822 // cout<<"nSig: "<<nSig<<endl;
1823  int aa1=0;
1824  int bb1=0;
1825  for (int n = 0; n <NNTree2->GetEntries(); n++){
1826  NNTree2->GetEntry(n);
1827  char line[2048];
1828  if(n<nSig&&n!=nSig) {
1829  Mc_inj_gMc->Fill(Mc_inj,Mc);
1830  Mc_inj_g->Fill(Mc_inj,av);
1831  gS[0]->SetPoint(n,Mc,av);
1832  cout<<" n "<<n<<" Mc "<<Mc<<" av "<<av<<" M_inj "<<Mc_inj<<endl;
1833  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1834  aa1=aa1+1;
1835  ANN_Mc_SIG3D->Fill(Mc,av);
1836  ANN_cc_SIG3D->Fill(cc,av);
1837  sprintf(line,"SIG:Average_%1.5f; Mc_rec_%1.5f; Mc_inj_%1.5f; event %i; count_%i; nSig_%i\n",av,Mc,Mc_inj,n,aa1,nSig);
1838  file_txt_out<<line;
1839  cout<<" n "<<n<<" nSig "<<nSig<<endl;
1840  }
1841  else {
1842  bb1=bb1+1;
1843  gB[0]->SetPoint(n-nSig,Mc,av);
1844  ANN_Mc_BKG3D->Fill(Mc,av);
1845  ANN_cc_BKG3D->Fill(cc,av);
1846  // cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1847  sprintf(line,"BKG:Average %1.5f; Mc_rec %1.5f; Mc_inj_%1.5f; event %i; count_%i\n",av,Mc,Mc_inj,n,bb1);
1848  file_txt_out<<line;
1849 
1850  }
1851  }
1852  //exit(0);
1853 
1854 
1855  gS[0]->SetMarkerColor(2);
1856  gB[0]->SetMarkerColor(4);
1857  gS[0]->SetMarkerStyle(6);
1858  gB[0]->SetMarkerStyle(7);
1859 
1860  ANN_Mc_c3D->cd(1);
1861  ANN_Mc_SIG3D->GetXaxis()->SetTitle("SIG:Mchirp_reconstructed");
1862  ANN_Mc_SIG3D->GetYaxis()->SetTitle("SIG:ANN_average");
1863  ANN_Mc_SIG3D->GetZaxis()->SetTitle("count");
1864  ANN_Mc_SIG3D->Draw("colz");
1865  ANN_Mc_c3D->cd(2);
1866  ANN_Mc_BKG3D->GetXaxis()->SetTitle("BKG:Mchirp_reconstructed");
1867  ANN_Mc_BKG3D->GetYaxis()->SetTitle("BKG:ANN_average");
1868  ANN_Mc_BKG3D->GetZaxis()->SetTitle("count");
1869  ANN_Mc_BKG3D->Draw("colz");
1870 
1871 
1872  TString ANN_Mc_name3D(name);
1873  TString ANN_Mc_rootname3D(name);
1874  char ANN_Mc_cname23D[1024];
1875  char ANN_Mc_crootname23D[1024];
1876  ANN_Mc_name3D.ReplaceAll(".root",".png");
1877  sprintf(ANN_Mc_cname23D,"average_png/nSig%i_nBg%i_3DMc_ANNav_Plots_%s",nSig,nBg,ANN_Mc_name3D.Data());
1878  sprintf(ANN_Mc_crootname23D,"average_png/nSig%i_nBg%i_3DMc_ANNav_Plots_%s",nSig,nBg,ANN_Mc_rootname3D.Data());
1879  ANN_Mc_c3D->SaveAs(ANN_Mc_cname23D);
1880  ANN_Mc_c3D->Print(ANN_Mc_crootname23D);
1881 
1882 
1883  ANN_cc_c3D->cd(1);
1884  ANN_cc_SIG3D->GetXaxis()->SetTitle("SIG:cc");
1885  ANN_cc_SIG3D->GetYaxis()->SetTitle("SIG:ANN_average");
1886  ANN_cc_SIG3D->GetZaxis()->SetTitle("count");
1887  ANN_cc_SIG3D->Draw("colz");
1888  ANN_cc_c3D->cd(2);
1889  ANN_cc_BKG3D->GetXaxis()->SetTitle("BKG:cc");
1890  ANN_cc_BKG3D->GetYaxis()->SetTitle("BKG:ANN_average");
1891  ANN_cc_BKG3D->GetZaxis()->SetTitle("count");
1892  ANN_cc_BKG3D->Draw("colz");
1893 
1894 
1895  TString ANN_cc_cname3D(name);
1896  TString ANN_cc_crootname3D(name);
1897  char ANN_cc_cname23D[1024];
1898  char ANN_cc_crootname23D[1024];
1899  ANN_cc_cname3D.ReplaceAll(".root",".png");
1900  sprintf(ANN_cc_cname23D,"average_png/nSig%i_nBg%i_3DANNav_cc_Plots_%s",nSig,nBg,ANN_cc_cname3D.Data());
1901  sprintf(ANN_cc_crootname23D,"average_png/nSig%i_nBg%i_3DANNav_cc_Plots_%s",nSig,nBg,ANN_cc_crootname3D.Data());
1902  ANN_cc_c3D->SaveAs(ANN_cc_cname23D);
1903  ANN_cc_c3D->Print(ANN_cc_crootname23D);
1904 
1905  Mc_inj_c->cd(1);
1906  Mc_inj_g->GetXaxis()->SetTitle("Mchirp_injected");
1907  Mc_inj_g->GetYaxis()->SetTitle("ANN_average");
1908  Mc_inj_g->GetZaxis()->SetTitle("count");
1909  Mc_inj_g->Draw("colz");
1910  Mc_inj_c->cd(2);
1911  Mc_inj_gMc->GetXaxis()->SetTitle("Mchirp_injected");
1912  Mc_inj_gMc->GetYaxis()->SetTitle("Mchirp_reconstructed");
1913  Mc_inj_gMc->GetZaxis()->SetTitle("count");
1914  Mc_inj_gMc->Draw("colz");
1915 
1916 
1917  TString Mc_inj_cname(name);
1918  TString Mc_inj_crootname(name);
1919  char Mc_inj_cname2[1024];
1920  char Mc_inj_crootname2[1024];
1921  Mc_inj_cname.ReplaceAll(".root",".png");
1922  sprintf(Mc_inj_cname2,"average_png/nSig%i_nBg%i_Mc_inj_Plots_%s",nSig,nBg,Mc_inj_cname.Data());
1923  sprintf(Mc_inj_crootname2,"average_png/nSig%i_nBg%i_Mc_inj_Plots_%s",nSig,nBg,Mc_inj_crootname.Data());
1924  Mc_inj_c->SaveAs(Mc_inj_cname2);
1925  Mc_inj_c->Print(Mc_inj_crootname2);
1926 
1927 
1928  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1929  c->Divide(1,2);
1930  c->cd(1);
1931  TMultiGraph* mg1=new TMultiGraph();
1932  mg1->SetTitle("Av_Mc");
1933  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1934  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1935  mg1->Draw("ap");
1936  mg1->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1937  mg1->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1938  c->cd(2);
1939  TMultiGraph* mg2=new TMultiGraph();
1940  mg2->SetTitle("Av_Mc");
1941  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1942  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1943  mg2->Draw("ap");
1944  mg2->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1945  mg2->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1946 
1947  cout<<" name "<<name<<endl;
1948  TString Cname(name);
1949  TString Cnameroot(name);
1950  char Cname2[1024];
1951  char Cname2root[1024];
1952  Cname.ReplaceAll(".root",".png");
1953  sprintf(Cname2,"average_png/nSig%i_nBg%i_Mc_Plots_%s",nSig,nBg,Cname.Data());
1954  sprintf(Cname2root,"average_png/nSig%i_nBg%i_Mc_Plots_%s",nSig,nBg,Cnameroot.Data());
1955  cout<<"Cname2 "<<Cname2<<endl;
1956  c->SaveAs(Cname2);
1957  c->Print(Cname2root);
1958  cout<<" gS[0]->GetN() "<<gS[0]->GetN()<<endl;
1959  //cout<<"entries "<<mg1->GetHistogram()->GetEntries()<<endl;
1960  cout<<" nSig "<<nSig<<" nBg "<<nBg<<" aa1 "<<aa1<<endl;
1961 }
1962 
1963 
1964 void SNR_plots(TString ifile){
1965 
1966  TString name(ifile);
1967  name.ReplaceAll("outfile/","");
1968  name.ReplaceAll("average_file/","");
1969  TFile* fTEST =TFile::Open(ifile.Data());
1970  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1971  double av;
1972  double out;
1973  double cc;
1974  double Mc;
1975  double Mc_inj;
1976  double rho;
1977  int nSi;
1978  NNTree2->SetBranchAddress("Mchirp",&Mc);
1979  NNTree2->SetBranchAddress("Average",&av);
1980  NNTree2->SetBranchAddress("cc",&cc);
1981  NNTree2->SetBranchAddress("rho",&rho);
1982  NNTree2->SetBranchAddress("Mchirp_injected",&Mc_inj);
1983  NNTree2->SetBranchAddress("#TestSig",&nSi);
1984  NNTree2->GetEntry(0);
1985  int const nSig=nSi;
1986  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1987  int const nBg=NNTree2->GetEntries()-nSig;
1988  int const nTot=NNTree2->GetEntries();
1989  int const n_points=N_ANNth*N_RHOth;
1990  int nBgSur[3][n_points];
1991  int nSigSur[3][n_points];
1992  double zax[3][n_points];
1993 
1994  for (int i=0; i<3; i++) for (int u=0; u<n_points; u++) {nBgSur[i][u]=0; nSigSur[i][u]=0;zax[i][u]=0.;}
1995 
1996  double pANNth[N_ANNth];
1997  double pRHOth[N_RHOth];
1998 
1999  double pCCth[3];
2000 
2001  pCCth[0]=0.;
2002  pCCth[1]=0.;
2003  pCCth[2]=0.;
2004 
2005  pCCth[0]=0.5;
2006  pCCth[1]=0.6;
2007  pCCth[2]=0.7;
2008 
2009  double deltaANNth=0.;
2010  double deltaRHOth=0.;
2011 
2012  deltaANNth=(double)(ANNth_max-ANNth_min)/N_ANNth;
2013  deltaRHOth=(double)(RHOth_max-RHOth_min)/N_RHOth;
2014 
2015 
2016  TCanvas* SNR_c=new TCanvas("SIG-BKG_reduction","SIG-BKG_reduction",0,0,1200,700);
2017  SNR_c->Divide(1,3);
2018 
2019  TH2D* SNR3Dcc0=new TH2D("cc=0.5","cc=0.5",N_ANNth,ANNth_min,ANNth_max,N_RHOth,RHOth_min,RHOth_max);
2020  TH2D* SNR3Dcc1=new TH2D("cc=0.6","cc=0.6",N_ANNth,ANNth_min,ANNth_max,N_RHOth,RHOth_min,RHOth_max);
2021  TH2D* SNR3Dcc2=new TH2D("cc=0.7","cc=0.7",N_ANNth,ANNth_min,ANNth_max,N_RHOth,RHOth_min,RHOth_max);
2022  SNR3Dcc0->SetDirectory(0);
2023  SNR3Dcc0->SetStats(0);
2024  SNR3Dcc1->SetDirectory(0);
2025  SNR3Dcc1->SetStats(0);
2026  SNR3Dcc2->SetDirectory(0);
2027  SNR3Dcc2->SetStats(0);
2028 
2029  for (int u=0; u<N_ANNth; u++) { pANNth[u]=0.; pANNth[u]=ANNth_min+u*deltaANNth; cout<<" pANNth[u] "<<pANNth[u]<<" u "<<u<<endl;}
2030  for (int u=0; u<N_RHOth; u++) { pRHOth[u]=0.; pRHOth[u]=RHOth_min+u*deltaRHOth; cout<<" pRHOth[u] "<<pRHOth[u]<<" u "<<u<<endl;}
2031 
2032 cout<<"fino def cose"<<endl;
2033 /*
2034  for (int u=0; u<nTot; u++){
2035  NNTree2->GetEntry(u);
2036  cout<<" cc "<<cc<<" av "<<av<<" rho "<< rho<<endl;
2037  for (int y=2; y>-1; y--){
2038  // cout<<" cc "<<cc<<" pCCth "<<pCCth[y]<<" y "<<y<<endl;
2039  if(cc>pCCth[y]) {
2040  for (int i=N_ANNth-1; i>-1; i--){
2041  // cout<<" av "<<av<<" pANNth[i] "<<pANNth[i]<<" i "<<i<<endl;
2042  if(av>pANNth[i]){
2043  for (int j=N_RHOth-1; j>-1; j++){
2044  //cout<<" rho "<<rho<<" pRHOth[j] "<<pRHOth[j]<<" j "<<j<<" nSig "<<nSig<<" u "<<u<<endl;
2045  if(rho>pRHOth[j]){
2046  for(int yy=y; yy>-1;yy--) for (int ii=i; ii>-1; ii--) for (int jj=j; jj>-1;jj--){
2047  if(u<nSig) nSigSur[yy][jj+N_ANNth*ii]=nSigSur[yy][jj+N_ANNth*ii]+1;
2048  else nBgSur[yy][jj+N_ANNth*ii]=nBgSur[yy][jj+N_ANNth*ii]+1;
2049  // cout<<" yy "<<yy<<" ii "<<ii<<" jj "<<jj<<" nBgSur[yy][jj+N_ANNth*ii] "<<nBgSur[yy][jj+N_ANNth*ii]<<" nSigSur[yy][jj+N_ANNth*ii] "<<nSigSur[yy][jj+N_ANNth*ii]<<endl;
2050  }
2051  break;
2052  }
2053  else continue;
2054  }
2055  break;
2056  }
2057  else continue;
2058  }
2059  break;
2060  }
2061  else continue;
2062  }
2063  }
2064 //exit(0);
2065 cout<<"dopo ciclo"<<endl;
2066 */
2067 for (int u=0; u<nTot; u++){
2068  NNTree2->GetEntry(u);
2069  for (int y=0; y<3; y++){
2070  if(cc>pCCth[y]){
2071  for (int i=0; i<N_ANNth; i++){
2072  if(av>pANNth[i]) {
2073  for (int j=0; j<N_RHOth; j++){
2074  if(rho>pRHOth[j]){
2075  if(u<nSig) nSigSur[y][j+N_ANNth*i]=nSigSur[y][j+N_ANNth*i]+1;
2076  else nBgSur[y][j+N_ANNth*i]=nBgSur[y][j+N_ANNth*i]+1;
2077  }
2078  else continue;
2079  }
2080  }
2081  else continue;
2082  }
2083  }
2084  else continue;
2085  }
2086  }
2087 
2088 
2089 //exit(0);
2090 int nBgSur0=0;
2091 nBgSur0=nBgSur[0][0];
2092 int nSigSur0=0;
2093 nSigSur0=nSigSur[0][0];
2094 if(nBgSur0==0){nBgSur0=nBg;}
2095 if(nSigSur0==0){nSigSur0=nSig;}
2096 int count_zax=0;
2097 cout<<"dopo ciclo"<<" nSigSur0 "<<nSigSur0<<" nBgSur0 "<<nBgSur0<<endl;
2098 for (int u=0; u<3; u++) for (int j=0; j<N_RHOth; j++) for (int i=0; i<N_ANNth; i++) {
2099  //cout<<" u "<<u<<" j "<<j<<" i "<<i<<" nSigSur[u][i+j*N_ANNth] "<<nSigSur[u][i+j*N_ANNth]<<endl;
2100  //if(nSigSur[u][i+j*N_ANNth]!=0){zax[u][i+j*N_ANNth]=(nBgSur[u][i+j*N_ANNth]/nBgSur0)/(nSigSur[u][i+j*N_ANNth]/nSigSur0);};
2101  if(nSigSur[u][i+j*N_ANNth]==0) {
2102  if (nSigSur[u][i+j*N_ANNth]==0 && nBgSur[u][i+j*N_ANNth]==0) zax[u][i+j*N_ANNth]=-1;
2103  else zax[u][i+j*N_ANNth]=-2;
2104  }
2105  else {
2106  double num=0.;
2107  double den=0.;
2108  num=(double)nBgSur[u][i+j*N_ANNth]/nBgSur0;
2109  den=(double)nSigSur0/nSigSur[u][i+j*N_ANNth];
2110  zax[u][i+j*N_ANNth]=(double)num/den; count_zax=count_zax+1;
2111  cout<<(double)(nBgSur[u][i+j*N_ANNth]/nBgSur0)*(nSigSur0/nSigSur[u][i+j*N_ANNth])<<endl;
2112  }
2113 // cout<<" nBgSur[yy][jj+N_ANNth*ii] "<<nBgSur[u][j+N_ANNth*i]<<" nSigSur[y][j+N_ANNth*i] "<<nSigSur[u][j+N_ANNth*i]<<" zax[u][i+j*N_ANNth] "<<zax[u][i+j*N_ANNth]<<endl;
2114  //else zax[u][i+j*N_ANNth]=(nBgSur[u][i+j*N_ANNth]/nBgSur0)/(nSigSur[u][i+j*N_ANNth]/nSigSur0);
2115  if(u==2){
2116  /*cout<<" zax[0][i+j*N_ANNth] "<<zax[0][i+j*N_ANNth]<<endl;
2117  cout<<" zax[1][i+j*N_ANNth] "<<zax[1][i+j*N_ANNth]<<endl;
2118  cout<<" zax[2][i+j*N_ANNth] "<<zax[2][i+j*N_ANNth]<<endl;*/
2119  SNR3Dcc0->Fill(pANNth[i]+ deltaANNth/2.,pRHOth[j]+ deltaRHOth/2.,zax[0][i+j*N_ANNth]);
2120  SNR3Dcc1->Fill(pANNth[i]+ deltaANNth/2.,pRHOth[j]+ deltaRHOth/2.,zax[1][i+j*N_ANNth]);
2121  SNR3Dcc2->Fill(pANNth[i]+ deltaANNth/2.,pRHOth[j]+ deltaRHOth/2.,zax[2][i+j*N_ANNth]);
2122  }
2123 }
2124 for(int yy=0; yy<3;yy++) {
2125  cout<<" pCCth "<<pCCth[yy]<<endl;
2126  for (int ii=N_ANNth-1; ii>-1; ii--) {
2127  cout<<" pANNth[i] "<<pANNth[ii]<<endl;
2128  for (int jj=N_RHOth-1; jj>-1;jj--){
2129  cout<<" pRHOth[jj] "<<pRHOth[jj]<<endl;
2130  cout<<" nBgSur[yy][ii+N_ANNth*jj] "<<nBgSur[yy][ii+N_ANNth*jj]<<endl;
2131  cout<<" nSigSur[yy][ii+N_ANNth*jj] "<<nSigSur[yy][ii+N_ANNth*jj]<<endl;
2132  cout<<" zax[0][i+j*N_ANNth]"<<zax[yy][ii+jj*N_ANNth]<<endl;
2133  cout<<" (double)(nBgSur[u][i+j*N_ANNth]/nBgSur0)"<<(double)nBgSur[yy][ii+jj*N_ANNth]/nBgSur0<<endl;
2134  //cout<<" yy "<<yy<<" ii "<<ii<<" jj "<<jj<<" nBgSur[yy][jj+N_ANNth*ii] "<<nBgSur[yy][jj+N_ANNth*ii]<<" nSigSur[yy][jj+N_ANNth*ii] "<<nSigSur[yy][jj+N_ANNth*ii]<<endl;
2135  }
2136  }
2137  }
2138 cout<<"coutn_zax "<<count_zax<<endl;
2139  SNR_c->cd(1);
2140  SNR3Dcc0->GetXaxis()->SetTitle("ANN_average_thresholds");
2141  SNR3Dcc0->GetYaxis()->SetTitle("RHO_thresholds");
2142  SNR3Dcc0->GetZaxis()->SetTitle("normBKG/normSIG");
2143  SNR3Dcc0->Draw("colz");
2144 
2145  SNR_c->cd(2);
2146  SNR3Dcc1->GetXaxis()->SetTitle("ANN_average_thresholds");
2147  SNR3Dcc1->GetYaxis()->SetTitle("RHO_thresholds");
2148  SNR3Dcc1->GetZaxis()->SetTitle("normBKG/normSIG");
2149  SNR3Dcc1->Draw("colz");
2150 
2151  SNR_c->cd(3);
2152  SNR3Dcc2->GetXaxis()->SetTitle("ANN_average_thresholds");
2153  SNR3Dcc2->GetYaxis()->SetTitle("RHO_thresholds");
2154  SNR3Dcc2->GetZaxis()->SetTitle("normBKG/normSIG");
2155  SNR3Dcc2->Draw("colz");
2156 
2157  TString SNR_n(name);
2158  TString SNR_nroot(name);
2159  char SNR_c_n[1024];
2160  char SNR_c_nroot[1024];
2161  SNR_n.ReplaceAll(".root",".png");
2162  sprintf(SNR_c_n,"SNR_plots/nSig%i_nBg%i_invSNR_Plots_%s",nSig,nBg,SNR_n.Data());
2163  sprintf(SNR_c_nroot,"SNR_plots/nSig%i_nBg%i_invSNR_Plots_%s",nSig,nBg,SNR_nroot.Data());
2164  SNR_c->SaveAs(SNR_c_n);
2165  SNR_c->Print(SNR_c_nroot);
2166 
2167 
2168 }
2169 
2170 void CoG_plots(TString ifile){
2171 
2172  TString name(ifile);
2173  name.ReplaceAll("outfile/","");
2174  name.ReplaceAll("average_file/","");
2175  TFile* fTEST =TFile::Open(ifile.Data());
2176  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
2177  double av;
2178  double out;
2179  double cc;
2180  double Mc;
2181  double Mc_inj;
2182  double rho;
2183  int nSi;
2184  int CoG;
2185  NNTree2->SetBranchAddress("Center_of_Gravity",&CoG);
2186  NNTree2->SetBranchAddress("Mchirp",&Mc);
2187  NNTree2->SetBranchAddress("Average",&av);
2188  NNTree2->SetBranchAddress("cc",&cc);
2189  NNTree2->SetBranchAddress("rho",&rho);
2190  NNTree2->SetBranchAddress("Mchirp_injected",&Mc_inj);
2191  NNTree2->SetBranchAddress("#TestSig",&nSi);
2192  NNTree2->GetEntry(0);
2193  int const nSig=nSi;
2194  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
2195  int const nTot=NNTree2->GetEntries();
2196  int const nBg=nTot-nSig;
2197 
2198  TCanvas* CoG_c=new TCanvas("SIG-BKG_reduction","SIG-BKG_reduction",0,0,1200,700);
2199 
2200  TH1D* imp=new TH1D("Center_of_Gravity_VS_ANN","Center_of_Gravity_VS_ANN",N_ANNth,ANNth_min,ANNth_max);
2201  TH1D* CoG1Sig=new TH1D("Center_of_Gravity_VS_ANN","Center_of_Gravity_VS_ANN",N_ANNth,ANNth_min,ANNth_max);
2202  TH1D* CoG1Bg=new TH1D("Center_of_Gravity_VS_ANN","Center_of_Gravity_VS_ANN",N_ANNth,ANNth_min,ANNth_max);
2203  TH1D* CoG0Sig=new TH1D("Center_of_Gravity_VS_ANN","Center_of_Gravity_VS_ANN",N_ANNth,ANNth_min,ANNth_max);
2204  TH1D* CoG0Bg=new TH1D("Center_of_Gravity_VS_ANN","Center_of_Gravity_VS_ANN",N_ANNth,ANNth_min,ANNth_max);
2205  TH1D* ibkgCC=new TH1D("Number_of_events_vs_CC","Number_of_events_vs_CC",50,0.,1.);
2206  TH1D* isigCC=new TH1D("Number_of_events_vs_CC","Number_of_events_vs_CC",50,0.,1.);
2207  TH1D* ibkgRHO=new TH1D("Number_of_events_vs_RHO","Number_of_events_vs_RHO",50,0.,50.);
2208  TH1D* isigRHO=new TH1D("Number_of_events_vs_RHO","Number_of_events_vs_RHO",50,0.,50.);
2209  TH1D* ibkgANN=new TH1D("Number_of_events_vs_ANN","Number_of_events_vs_ANN",50,-0.2,1.2);
2210  TH1D* isigANN=new TH1D("Number_of_events_vs_ANN","Number_of_events_vs_ANN",50,-0.2,1.2);
2211 
2212  CoG1Sig->SetDirectory(0);
2213  CoG1Sig->SetStats(0);
2214  CoG0Sig->SetDirectory(0);
2215  CoG0Sig->SetStats(0);
2216  CoG1Bg->SetDirectory(0);
2217  CoG1Bg->SetStats(0);
2218  CoG0Bg->SetDirectory(0);
2219  CoG0Bg->SetStats(0);
2220  ibkgCC->SetDirectory(0);
2221  ibkgCC->SetStats(0);
2222  isigRHO->SetDirectory(0);
2223  isigRHO->SetStats(0);
2224  ibkgANN->SetDirectory(0);
2225  ibkgANN->SetStats(0);
2226  isigCC->SetDirectory(0);
2227  isigCC->SetStats(0);
2228  ibkgRHO->SetDirectory(0);
2229  ibkgRHO->SetStats(0);
2230  isigANN->SetDirectory(0);
2231  isigANN->SetStats(0);
2232 
2233  for (int n=0; n<nTot; n++){
2234  NNTree2->GetEntry(n);
2235  if (n<nSig) {
2236  if(CoG==0) {CoG0Sig->Fill(av);}
2237  else {CoG1Sig->Fill(av);}
2238  isigCC->Fill(cc);
2239  isigRHO->Fill(rho);
2240  isigANN->Fill(av);
2241  }
2242  else {
2243  if(CoG==0) {CoG0Bg->Fill(av);}
2244  else {CoG1Bg->Fill(av);}
2245  ibkgCC->Fill(cc);
2246  ibkgRHO->Fill(rho);
2247  ibkgANN->Fill(av);
2248 
2249  }
2250  cout<<" CoG "<<CoG<<endl;
2251  }
2252 
2253  imp->SetDirectory(0);
2254  imp->SetStats(0);
2255  int max_n=0;
2256  // if(nBg>nSig) max_n=nBg;
2257  // else max_n=nSig;
2258 
2259  int max[4];
2260  for(int i=0; i<4; i++) max[i]=0;
2261 
2262  max[0]=CoG1Sig->GetMaximum();
2263  max[1]=CoG0Sig->GetMaximum();
2264  max[2]=CoG1Bg->GetMaximum();
2265  max[3]=CoG0Bg->GetMaximum();
2266 
2267 
2268 max_n=max[0];
2269  for(int i=0; i<4; i++) if(max_n<max[i]) max_n=max[i] ;
2270 
2271  imp->Fill(ANNth_min,max_n);
2272  imp->SetLineColor(10);
2273  imp->SetBins(N_ANNth,ANNth_min,ANNth_max);
2274  imp->GetXaxis()->SetTitle("ANN");
2275  imp->GetYaxis()->SetTitle("Count");
2276  imp->Draw();
2277 
2278  isigCC->SetFillStyle(3018);
2279  isigRHO->SetFillStyle(3018);
2280  isigANN->SetFillStyle(3018);
2281  ibkgCC->SetFillStyle(3004);
2282  ibkgRHO->SetFillStyle(3004);
2283  ibkgANN->SetFillStyle(3004);
2284  CoG1Sig->SetFillStyle(3018);
2285  CoG0Sig->SetFillStyle(3018);
2286  CoG0Bg->SetFillStyle(3004);
2287  CoG1Bg->SetFillStyle(3004);
2288  ibkgCC->SetLineColor(4);//blu
2289  ibkgANN->SetLineColor(4);//blu
2290  ibkgRHO->SetLineColor(4);//blu
2291  isigCC->SetLineColor(2);//rosso
2292  isigANN->SetLineColor(2);//rosso
2293  isigRHO->SetLineColor(2);//rosso
2294  ibkgCC->SetFillColor(4);//blu
2295  ibkgANN->SetFillColor(4);//blu
2296  ibkgRHO->SetFillColor(4);//blu
2297  isigCC->SetFillColor(2);//rosso
2298  isigANN->SetFillColor(2);//rosso
2299  isigRHO->SetFillColor(2);//rosso
2300  CoG1Sig->SetFillColor(2);//rosso
2301  CoG0Sig->SetFillColor(6);//fucsia
2302  CoG0Bg->SetFillColor(4);//blu
2303  CoG1Bg->SetFillColor(9);//violaceo
2304  CoG1Sig->SetLineColor(2);
2305  CoG0Sig->SetLineColor(6);
2306  CoG0Bg->SetLineColor(4);
2307  CoG1Bg->SetLineColor(9);
2308  // CoG1Sig->SetBins(N_ANNth,ANNth_min,ANNth_max);
2309  // CoG1Sig->GetXaxis()->SetTitle("ANN");
2310  // CoG1Sig->GetYaxis()->SetTitle("Count");
2311 CoG_c->cd(1);
2312  CoG1Sig->Draw("same");
2313  CoG0Sig->Draw("same");
2314  CoG0Bg->Draw("same");
2315  CoG1Bg->Draw("same");
2316 
2317 
2318 
2319  TString CoG_n(name);
2320  TString CoG_nroot(name);
2321  char CoG_c_n[1024];
2322  char CoG_c_nroot[1024];
2323  CoG_n.ReplaceAll(".root",".png");
2324  sprintf(CoG_c_n,"CoG_plots/nSig%i_nBg%i_CoG_Plots_%s",nSig,nBg,CoG_n.Data());
2325  sprintf(CoG_c_nroot,"CoG_plots/nSig%i_nBg%i_CoG_Plots_%s",nSig,nBg,CoG_nroot.Data());
2326  CoG_c->SaveAs(CoG_c_n);
2327  CoG_c->Print(CoG_c_nroot);
2328 
2329 TCanvas* Istogram=new TCanvas("SIG-BKG_reduction","SIG-BKG_reduction",0,0,1200,400);
2330 
2331 Istogram->Divide(3,1);
2332 Istogram->cd(1);
2333  if (ibkgCC->GetMaximum()>isigCC->GetMaximum()) {ibkgCC->GetXaxis()->SetTitle("Network Correlation Coefficient"); ibkgCC->GetYaxis()->SetTitle("Count"); ibkgCC->Draw(); isigCC->Draw("same");}
2334  else {isigCC->GetXaxis()->SetTitle("Network Correlation Coefficient"); isigCC->GetXaxis()->SetTitle("Count"); isigCC->Draw();ibkgCC->Draw("same");}
2335 Istogram->cd(2);
2336  if (ibkgRHO->GetMaximum()>isigRHO->GetMaximum()) {ibkgRHO->GetXaxis()->SetTitle("Effective Correlated SNR"); ibkgRHO->GetYaxis()->SetTitle("Count"); ibkgCC->Draw(); ibkgRHO->Draw(); isigRHO->Draw("same");}
2337  else {isigRHO->GetXaxis()->SetTitle("Effective Correlated SNR"); isigRHO->GetXaxis()->SetTitle("Count"); isigRHO->Draw();ibkgRHO->Draw("same");}
2338 Istogram->cd(3);
2339  if (ibkgANN->GetMaximum()>isigANN->GetMaximum()) {ibkgANN->GetXaxis()->SetTitle("ANN parameter"); ibkgANN->GetYaxis()->SetTitle("Count"); ibkgANN->Draw(); isigANN->Draw("same");}
2340  else {isigANN->GetXaxis()->SetTitle("ANN parameter"); isigANN->GetYaxis()->SetTitle("Count"); isigANN->Draw();ibkgANN->Draw("same");}
2341 
2342  TString isto_n(name);
2343  TString isto_nroot(name);
2344  char isto_c_n[1024];
2345  char isto_c_nroot[1024];
2346  isto_n.ReplaceAll(".root",".png");
2347  sprintf(isto_c_n,"Istograms/nSig%i_nBg%i_istograms_%s",nSig,nBg,isto_n.Data());
2348  sprintf(isto_c_nroot,"Istograms/nSig%i_nBg%i_istograms_%s",nSig,nBg,isto_nroot.Data());
2349  Istogram->SaveAs(isto_c_n);
2350  Istogram->Print(isto_c_nroot);
2351 
2352 }
2353 
2354 void duration_plot(TString ifile){
2355 TString name(ifile);
2356  name.ReplaceAll("outfile/","");
2357  name.ReplaceAll("average_file/","");
2358  TFile* fTEST =TFile::Open(ifile.Data());
2359  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
2360  double av;
2361  double out;
2362  double cc;
2363  double Mc;
2364  double Mc_inj;
2365  double rho;
2366  int nSi;
2367  int CoG;
2368  double Dt;
2369 // float Dt;
2370 
2371  NNTree2->SetBranchAddress("Center_of_Gravity",&CoG);
2372  NNTree2->SetBranchAddress("duration",&Dt);
2373  NNTree2->SetBranchAddress("Mchirp",&Mc);
2374  NNTree2->SetBranchAddress("Average",&av);
2375  NNTree2->SetBranchAddress("cc",&cc);
2376  NNTree2->SetBranchAddress("rho",&rho);
2377  NNTree2->SetBranchAddress("Mchirp_injected",&Mc_inj);
2378  NNTree2->SetBranchAddress("#TestSig",&nSi);
2379  NNTree2->GetEntry(0);
2380  int const nSig=nSi;
2381  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
2382  int const nTot=NNTree2->GetEntries();
2383  int const nBg=nTot-nSig;
2384 
2385  TCanvas* Dt_c=new TCanvas("Duration","Duration",0,0,1200,700);
2386  double maxDt=0.;
2387  double minDt=0.;
2388  maxDt=Dt;
2389  minDt=Dt;
2390  cout<<"maxDt "<<maxDt<<" minDt "<<minDt<<endl;
2391 
2392  for( int i=0; i<nTot; i++){
2393  NNTree2->GetEntry(i);
2394  if(Dt>maxDt) maxDt=Dt;
2395  if(Dt<minDt) minDt=Dt;
2396  }
2397 
2398  /*TH1D* imp=new TH1D("Duration","Duration",100,minDt,maxDt);
2399  TH1D* Dt_hB=new TH1D("Duration","Duration",100,minDt,maxDt);
2400  TH1D* Dt_hS=new TH1D("Duration","Duration",100,minDt,maxDt);*/
2401  TH1D* imp=new TH1D("Duration","Duration",100,0.,2.);
2402  TH1D* Dt_hB=new TH1D("Duration","Duration",100,0.,2.);
2403  TH1D* Dt_hS=new TH1D("Duration","Duration",100,0.,2.);
2404 
2405  imp->SetDirectory(0);
2406  imp->SetStats(0);
2407  Dt_hS->SetDirectory(0);
2408  Dt_hS->SetStats(0);
2409  Dt_hB->SetDirectory(0);
2410  Dt_hB->SetStats(0);
2411 
2412  for (int n=0; n<nTot; n++){
2413  NNTree2->GetEntry(n);
2414  if (n<nSig) Dt_hS->Fill(Dt);
2415  else Dt_hB->Fill(Dt);
2416  }
2417 
2418  Dt_hS->SetFillStyle(3018);
2419  Dt_hB->SetFillStyle(3004);
2420  Dt_hS->SetFillColor(2);
2421  Dt_hB->SetFillColor(4);
2422  Dt_hS->SetLineColor(2);
2423  Dt_hB->SetLineColor(4);
2424 
2425  int max_n=0;
2426  int max[2];
2427  for(int i=0; i<2; i++) max[i]=0;
2428 
2429  max[0]=Dt_hS->GetMaximum();
2430  max[1]=Dt_hB->GetMaximum();
2431 
2432  max_n=max[0];
2433  for(int i=0; i<2; i++) if(max_n<max[i]) max_n=max[i] ;
2434 
2435  imp->Fill(minDt,max_n);
2436  imp->SetLineColor(10);
2437 // imp->SetBins(N_ANNth,ANNth_min,ANNth_max);
2438  imp->GetXaxis()->SetTitle("Duration");
2439  imp->GetYaxis()->SetTitle("Count");
2440  imp->Draw();
2441  Dt_hS->Draw("same");
2442  Dt_hB->Draw("same");
2443 
2444  TString CoG_n(name);
2445  TString CoG_nroot(name);
2446  char CoG_c_n[1024];
2447  char CoG_c_nroot[1024];
2448  CoG_n.ReplaceAll(".root",".png");
2449  sprintf(CoG_c_n,"Duration_plots/nSig%i_nBg%i_invSNR_Plots_%s",nSig,nBg,CoG_n.Data());
2450  sprintf(CoG_c_nroot,"Duration_plots/nSig%i_nBg%i_invSNR_Plots_%s",nSig,nBg,CoG_nroot.Data());
2451  Dt_c->SaveAs(CoG_c_n);
2452  Dt_c->Print(CoG_c_nroot);
2453 
2454 }
nCC
#define nCC
Definition: Average_FAout_moreGraphs_3.C:33
minCC
#define minCC
Definition: Average_FAout_moreGraphs_3.C:26
rho
double rho
Definition: cwb_report_prod_2.C:149
NMc_r
#define NMc_r
Definition: Average_FAout_moreGraphs_3.C:43
Ncc3D
#define Ncc3D
Definition: Average_FAout_moreGraphs_3.C:24
nRHO
#define nRHO
Definition: Average_FAout_moreGraphs_3.C:28
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:93
cwb_online.f
tuple f
Definition: cwb_online.py:91
netevent::duration
Float_t * duration
max cluster time relative to segment start
Definition: netevent.hh:79
a
wavearray< double > a(hp.size())
name
par[0] name
Definition: CWB_Plugin_HEN_SIM_Config.C:87
n
int n
Definition: cwb_net.C:10
Mcth
void Mcth(TString ifile, int &McFAth_n, int &McFDth_n, double &McthresFA)
Definition: Average_FAout_moreGraphs_3.C:1341
TString
TString("c")
Definition: cwb_report_skymap.C:111
x
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
Definition: BaseFunctionWDM.C:51
netevent::GetEntry
Int_t GetEntry(Int_t)
Definition: netevent.cc:387
Get
cout<<"Comparison bkg events above threshold: "<< entriesx<< endl;double *rhox=cnet.GetV1();comph=(TH1D *) gDirectory-> Get("hcomp")
Definition: compare_bkg.C:201
deltaANN
#define deltaANN
Definition: Average_FAout_moreGraphs_3.C:30
SNR_plots
void SNR_plots(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:1964
std
STL namespace.
Dt
TH1F * Dt
Definition: cbc_plots.C:398
ofile
char ofile[512]
Definition: cwb_condor_rescue.C:40
j
int j
Definition: cwb_net.C:10
i
i drho i
Definition: cwb_epparameters.C:85
graph
void graph(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:744
TB
CWB::Toolbox TB
Definition: ComputeSNR.C:5
out
ofstream out
Definition: cwb_merge.C:196
NDIM
#define NDIM
Definition: ClusterToFrameNew01_dtmin_approx_cent.C:31
maxMc_inj
#define maxMc_inj
Definition: Average_FAout_moreGraphs_3.C:42
h
wavearray< double > h
Definition: Regression_H1.C:25
RHOth_min
#define RHOth_min
Definition: Average_FAout_moreGraphs_3.C:56
nANN
#define nANN
Definition: Average_FAout_moreGraphs_3.C:29
PlotsAv_Mc
void PlotsAv_Mc(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:1743
iMc
#define iMc
Definition: Average_FAout_moreGraphs_3.C:39
background
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 \n");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
Definition: cbc_plots_sim4_ifar.C:1070
ANNth_max
#define ANNth_max
Definition: Average_FAout_moreGraphs_3.C:53
output
char output[256]
Definition: cwb_report_prod_2.C:876
k
int k
Definition: cwb_report_prod_2.C:150
minMc_inj
#define minMc_inj
Definition: Average_FAout_moreGraphs_3.C:41
netevent::netcc
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:94
CoG_plots
void CoG_plots(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:2170
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
num
int num
Definition: check_slag_bkg.C:71
yy
wavearray< double > yy
Definition: TestFrame5.C:12
maxMc_r
#define maxMc_r
Definition: Average_FAout_moreGraphs_3.C:44
s
s s
Definition: cwb_net.C:137
line
TLine * line
Definition: compare_bkg.C:482
minANN_av
#define minANN_av
Definition: Average_FAout_moreGraphs_3.C:47
NANN_av
#define NANN_av
Definition: Average_FAout_moreGraphs_3.C:46
Annth
void Annth(TString ifile, int &FAth_n, int &FDth_n, double &thresFA)
Definition: Average_FAout_moreGraphs_3.C:1071
deltarho
#define deltarho
Definition: Average_FAout_moreGraphs_3.C:32
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
nIFO
#define nIFO
Definition: Average_FAout_moreGraphs_3.C:27
RHOth_max
#define RHOth_max
Definition: Average_FAout_moreGraphs_3.C:55
minMc_r
#define minMc_r
Definition: Average_FAout_moreGraphs_3.C:45
PlotsAv_cc
void PlotsAv_cc(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:1616
netevent::GetEntries
Int_t GetEntries()
Definition: netevent.cc:381
N_ANNth
#define N_ANNth
Definition: Average_FAout_moreGraphs_3.C:51
N_RHOth
#define N_RHOth
Definition: Average_FAout_moreGraphs_3.C:52
Average_FAout_moreGraphs_3
void Average_FAout_moreGraphs_3(double FAdes, TString NN_FILE, TString NN_FILE2, TString NN_FILE3, TString NN_FILE4, TString NN_FILE5, TString NN_FILE6, TString NN_FILE7, TString NN_FILE8, TString TEST_FILE, TString NOMEtot_S, int TS=0, int TB=0, int s=0, int b=0, int uf=1, int consider_all=0, int av_on_nNN=0)
Definition: Average_FAout_moreGraphs_3.C:76
y
wavearray< double > y
Definition: Test10.C:31
ANNth_min
#define ANNth_min
Definition: Average_FAout_moreGraphs_3.C:54
NMc_inj
#define NMc_inj
Definition: Average_FAout_moreGraphs_3.C:40
deltacc
#define deltacc
Definition: Average_FAout_moreGraphs_3.C:31
netevent::chirp
Float_t * chirp
range to source: [0/1]-rec/inj
Definition: netevent.hh:109
duration_plot
void duration_plot(TString ifile)
Definition: Average_FAout_moreGraphs_3.C:2354
exit
exit(0)
maxANN_av
#define maxANN_av
Definition: Average_FAout_moreGraphs_3.C:48
maxCC
#define maxCC
Definition: Average_FAout_moreGraphs_3.C:25