Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Average.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 nIFO 3
25 #define nRHO 3
26 #define nANN 10
27 #define deltaANN 0.1
28 #define nCC 4
29 #define NPIX 20
30 #define RHOth 5.5
31 #define nth 10
32 #define rhomin 6
33 #define ccmin 0.6
34 #define iMc 1
35 using namespace std;
36 void graph(TString ifile);
37 void Annth(TString ifile);
38 void PlotsAv_cc(TString ifile);
39 void PlotsAv_Mc(TString ifile);
40 //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){
41 
42 //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").
43 //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").
44 //TS, TB: numbers of events to test for each class.
45 //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
46 //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.
47 void Average(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=1){
48 int ni=0;
49 if(uf==0) ni=1;
50 else ni=0;
51 
52 if (uf!=0&&consider_all==0) ni=1;
53 
54 //COUNT THE ANNs
55 int n_NN=0;
56 if (NN_FILE.CompareTo("")){
57  n_NN=n_NN+1;
58  }
59 if (NN_FILE2.CompareTo("")){
60  n_NN=n_NN+1;
61  }
62 if (NN_FILE3.CompareTo("")){
63  n_NN=n_NN+1;
64  }
65 if (NN_FILE4.CompareTo("")){
66  n_NN=n_NN+1;
67  }
68 if (NN_FILE5.CompareTo("")){
69  n_NN=n_NN+1;
70  }
71 if (NN_FILE6.CompareTo("")){
72  n_NN=n_NN+1;
73  }
74 if (NN_FILE7.CompareTo("")){
75  n_NN=n_NN+1;
76  }
77 if (NN_FILE8.CompareTo("")){
78  n_NN=n_NN+1;
79  }
80 int const nNN=n_NN;
81 
82 // int p=0;
83  char NNi[nNN][1024];
84  char NNi2[nNN][1024];
85  sprintf(NNi2[0],"%s",NN_FILE.Data());
86  if(nNN>1) sprintf(NNi2[1],"%s",NN_FILE2.Data());
87  if(nNN>2) sprintf(NNi2[2],"%s",NN_FILE3.Data());
88  if(nNN>3) sprintf(NNi2[3],"%s",NN_FILE4.Data());
89  if(nNN>4) sprintf(NNi2[4],"%s",NN_FILE5.Data());
90  if(nNN>5) sprintf(NNi2[5],"%s",NN_FILE6.Data());
91  if(nNN>6) sprintf(NNi2[6],"%s",NN_FILE7.Data());
92  if(nNN>7) sprintf(NNi2[7],"%s",NN_FILE8.Data());
93 
94 //saving the ANN names
95  for (int u=0;u<nNN;u++){
96  int p=0;
97  while (NNi2[u][p]){
98  if (NNi2[u][p]=='N'){
99  if((NNi2[u][p+1]=='N')&&(NNi2[u][p+2]=='\/')) {
100  int hh=p+3;
101  while (NNi2[u][hh]!='\0'){NNi[u][hh-p-3]=NNi2[u][hh];hh=hh+1;}
102  break;
103  }
104  }
105  p=p+1;
106  }
107 }
108 
109 //saving the Test-file name
110 int q=0;
111  char Filei[1024];
112  char Filei2[1024];
113  sprintf(Filei2,"%s",TEST_FILE.Data());
114  while (Filei2[q]){
115  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]=='\/')) {
116  int hh=q+7;
117  while (Filei2[hh]!='\0') {Filei[hh-q-7]=Filei2[hh];hh=hh+1;
118  }
119  for (int h0=hh-q-7;h0<1024;h0++) Filei[h0]='\0';
120  break;
121  }
122  q=q+1;
123  }
124 
125  cout<<Filei<<" original String: "<<Filei2<<endl;
126 
127 
128  TString NNi0[nNN];
129  NNi0[0]=NNi[0];
130  if(nNN>1) NNi0[1]=NNi[1];
131  if(nNN>2) NNi0[2]=NNi[2];
132  if(nNN>3) NNi0[3]=NNi[3];
133  if(nNN>4) NNi0[4]=NNi[4];
134  if(nNN>5) NNi0[5]=NNi[5];
135  if(nNN>6) NNi0[6]=NNi[6];
136  if(nNN>7) NNi0[7]=NNi[7];
137 
138  TString Filei0(Filei);
139 
140 //removing the extensions
141 for (int u=0;u<nNN;u++){
142  NNi0[u].ReplaceAll(".root","");
143  }
144  Filei0.ReplaceAll(".root","");
145 
146 //Definitions of some parameters
147  TFile *fnet[nNN];
148  TMultiLayerPerceptron* mlp[nNN];
149  TTree* infot[nNN];
150  int TotBg[nNN];
151  int FD0[nNN]; //FD=False Dismissal
152  int FA0[nNN]; //FA=False Alarm
153 
154  for (int yy=0; yy<nNN; yy++){TotBg[yy]=0;FA0[yy]=0;FD0[yy]=0;}
155 
156  int NNs[nNN]; //first SIG-event considered for ANN trainings
157  int NNb[nNN]; //first BKG-event considered for ANN trainings
158  int NNnTS[nNN]; //number of SIG-events considered for ANN trainings
159  int NNnTB[nNN]; //number of BKG-events considered for ANN trainings
160 
161 int min_NNb=5000000;
162 int min_NNs=5000000;
163 
164 //Extracting ANNs and information from the files
165 //if (uf!=0&&ni==0&&(TS!=0||TB!=0)){
166  for (int u=0;u<nNN;u++){
167  fnet[u]=TFile::Open(NNi2[u]);
168  mlp[u]=(TMultiLayerPerceptron*)fnet[u]->Get("TMultiLayerPerceptron");
169 
170  cout<<"dopo TMLP"<<endl;
171  if(mlp[u]==NULL) {cout << "Error getting mlp" << endl;exit(1);}
172  infot[u]=(TTree*)fnet[u]->Get("info");
173 
174  infot[u]->SetBranchAddress("Rand_start_Sig",&NNs[u]);
175  infot[u]->SetBranchAddress("Rand_start_Bg",&NNb[u]);
176  infot[u]->SetBranchAddress("#trainSig",&NNnTS[u]);
177  infot[u]->SetBranchAddress("#trainBg",&NNnTB[u]);
178  infot[u]->GetEntry(0);
179  cout<<"n "<<u<<" b "<<NNb[u]<<" min_NNb "<<min_NNb<<endl;
180 
181  if(NNs[u]<min_NNs) min_NNs=NNs[u];
182  if(NNb[u]<min_NNb) min_NNb=NNb[u];
183  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;
184  }
185 
186 //check if b and s have correct values.
187 int b_ok=0;
188 int s_ok=0;
189 for (int u=0;u<nNN;u++){
190  if((NNb[u]<b||NNb[u]==b)&&b>NNb[u]+ NNnTB[u]) b_ok=1;
191  if((NNs[u]<s||NNs[u]==s)&&b>NNs[u]+ NNnTS[u]) s_ok=1;
192 
193 }
194 
195 //Changing b value to its sum to the first event considered for the training
196 //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;}
197 if (uf!=0&&ni==0&&(b_ok==0)){ b=min_NNb;cout<<" change b value to have BKG events used for training procedures, b="<<b<<endl;}
198 if (uf!=0&&ni==0&&(s_ok==0)){ s=min_NNs;cout<<" change s value to have SIG events used for training procedures s="<<s<<endl;}
199 cout<<" min_NNb "<<min_NNb<<" b "<<b<<endl;
200 cout<<"Def. tree e mlp"<<endl;
201 
202 //opening the Test-file
203  TFile* fTEST =TFile::Open(TEST_FILE.Data());
204  TTree* NNTree=(TTree*)fTEST->Get("nnTree");
205  int entries=NNTree->GetEntries();
206  cout<<"entries: "<<entries<<endl;
207 
208 //Extract information necessary for the tree which defines the mlp
209  int ndim;
210  int ninp;
211  int y;
212  int entriesTot;
213  int bg_entries;
214  int sig_entries;
215 
216  NNTree->SetBranchAddress("#Entries_type",&entriesTot);
217  NNTree->SetBranchAddress("Matrix_dim",&ndim);
218  NNTree->SetBranchAddress("#inputs",&ninp);
219  NNTree->SetBranchAddress("amplitude_mode",&y);
220 
221  NNTree->GetEntry(0);
222  sig_entries=entriesTot;
223  NNTree->GetEntry(entries-1);
224  bg_entries=entriesTot;
225 
226  int const NDIM=ndim;
227  int const nINP=ninp;
228 
229  cout<<"NDIM: "<<NDIM<<endl;
230  cout<<"nINP: "<<nINP<<endl;
231  cout<<"sig e: "<<sig_entries<<endl;
232  cout<<"bg e: "<<bg_entries<<endl;
233 
234 //defining the sample which has less events
235  int minevents=0;
236  if (sig_entries>bg_entries) minevents=bg_entries;
237  else minevents=sig_entries;
238 
239  //if(b==0) b=sig_entries;
240  if(TB==0 && TS==0){
241  TS=minevents;
242  TB=minevents;
243  }
244 
245 //defining compatible parameters
246  if (b<sig_entries){
247  int a =b;
248  b+=sig_entries;
249  cout<<"Error: Bg index<sig_entries: new set of b parameter-> b="<<b<<" instead of b="<<a<<endl;
250  //exit(0);
251 }
252 
253 //set TS and TB to their maximum posible values
254 if((TS>sig_entries-s||TB>(bg_entries-(b-sig_entries)))&&(TS==TB)) {TS=minevents-s;TB=minevents-(b-sig_entries);
255  if (TS<TB) TB=TS;
256  else TS=TB;
257  cout<<"Error:S>sig_entries or TB>bg_entries, to maintain ugual number of analysed events TS=TB="<<TB<<endl;
258  }
259  if((TS>sig_entries-s||TB>bg_entries-(b-sig_entries))&&(TS!=TB)) {
260  if(TS>sig_entries) TS=sig_entries-s;
261  if (TB>bg_entries) TB=bg_entries-(b-sig_entries);
262  cout<<"Error:TS>sig_entries or TB>bg_entries, new TS and TB values are thus define-> TS="<<TS<<" TB="<<TB<<endl;
263  }
264 
265 
266  char NOMEtot[1024];
267  sprintf(NOMEtot,"%s",NOMEtot_S.Data());
268  cout<<"output name: "<<NOMEtot<<endl;
269 
270 //extracting original files
271  TString SIG_FILE;
272  TString BG_FILE;
273  char FILE_NAME[516];
274  NNTree->SetBranchAddress("Files_name",&FILE_NAME);
275  NNTree->GetEntry(0);
276  SIG_FILE=FILE_NAME;
277  NNTree->GetEntry(entries-1);
278  BG_FILE=FILE_NAME;
279  cout<<"fine ifdef RHO_CC"<<endl;
280 
281 //extracting other information
282  TChain sigTree("waveburst");//search the Tree "waveburst" in files
283  sigTree.Add(SIG_FILE.Data());
284  netevent signal(&sigTree,nIFO);
285  int sig_entries2 = signal.GetEntries();
286  cout << "sig entries2 : " << sig_entries2 << endl;
287 
288  TChain bgTree("waveburst");
289  bgTree.Add(BG_FILE.Data());
290  netevent background(&bgTree,nIFO);
291  int bg_entries2 = background.GetEntries();
292  cout << "bg entries2 : " << bg_entries2 << endl;
293 
294  cout<<"b: "<<b<<endl;
295  cout<<"s: "<<s<<endl;
296 
297  // add leaf
298  Float_t x[nINP];
299  for (int jj=0; jj<nINP;jj++) x[jj]=0.;
300  char ilabel[nINP][16];
301 
302  // define a branch for each input
303  for(int i=0;i<nINP;i++) {
304  sprintf(ilabel[i],"x%i",i+1);
305  NNTree->SetBranchAddress(ilabel[i], &x[i]);
306  }
307 
308 char ofile[1024];
309 sprintf(ofile,"average_file/%s.root",NOMEtot);
310 TFile*f =new TFile(ofile,"RECREATE");
311  TTree* NNTree2=new TTree("Parameters","Parameters");
312  NNTree2->SetDirectory(f);
313 
314 // defining useful information
315  double average=0.;
316  double out[nNN];
317  for(int y=0; y<nNN; y++) out[y]=0.;
318  double Mc=0.;
319  double NNcc=0.;
320  double NNrho=0.;
321  double std=0.;
322 
323 //adding new information to the original tree
324  NNTree2->Branch("Average",&average,"Average/D");
325  NNTree2->Branch("StandardDevaition",&std,"StandardDeviation/D");
326  NNTree2->Branch("cc",&NNcc,"cc/D");
327  NNTree2->Branch("Mchirp",&Mc,"Mchirp/D");
328  NNTree2->Branch("rho",&NNrho,"rho/D");
329 
330  for(int u=0;u<nNN;u++){
331  char NNoutl[512];
332  char NNoutl2[512];
333  char NNnamel[512];
334  char NNnamel2[512];
335  sprintf(NNoutl,"NNout%i",u);
336  sprintf(NNoutl2,"NNout%i/D",u);
337  sprintf(NNnamel,"NNname%i",u);
338  sprintf(NNnamel2,"NNname%i/C",u);
339  NNTree2->Branch(NNoutl,&out[u],NNoutl2);
340  NNTree2->Branch(NNnamel,&NNi[u],NNnamel2);
341  }
342 
343  char Testf[1024];
344  NNTree2->Branch("TestFile",&Testf,"TestFile/C");
345  sprintf(Testf,"%s",TEST_FILE.Data());
346  int nTestS=0;
347 
348  NNTree2->Branch("#TestSig",&nTestS,"#TestSig/I");
349  cout<<"nTestS: "<<nTestS<<" TS: "<<TS<<endl;
350  cout<<"dopo def tree"<<endl;
351 
352 //defining how many examples to use
353  int scount=0;
354  if(uf==0) nTestS=TS;
355  else {
356  for(int n=s;n<s+TS;n++) {
357  int countNN=0;
358  for(int u=0;u<nNN;u++){
359  if (uf!=0&&n>NNs[u]&&n<(NNs[u]+NNnTS[u])) countNN=countNN+1;
360  }
361  if(countNN==1&&ni==0) scount=scount+1;
362  if(countNN<2&&ni!=0) scount=scount+1; //if you are indifferently intereseted in average on nNN or nNN-1
363  if(countNN>1){cout<<"Error: training non independent"<<n<<endl; exit(0);} //exit if the ANNs are training on dependent set of events
364  }
365  nTestS=scount;
366 }
367 
368 cout<<"s test "<<s<< " s+TS "<<s+TS<<" b "<<b<<" b+ TB "<<b+TB<<" nTestS "<<nTestS<<endl;
369 
370  double params[nINP];
371  int sig_05=0;
372  for (int i=0; i<nINP; i++) params[i]=0.;
373 
374  int S_eff=0;
375  int FD=0;
376 
377  for(int n=s;n<s+TS;n++) {
378  average=0.;
379  int indNN=-1;
380  int countNN=0;
381  for(int u=0;u<nNN;u++){
382  if (uf!=0&&n>=NNs[u]&&n<(NNs[u]+NNnTS[u])) { indNN=u;countNN=countNN+1;}
383  }
384  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
385  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
386  if(ni==0&&countNN==0)continue; //skip if any network is trained on the considered event
387  NNTree->GetEntry(n);
388  signal.GetEntry(n);
389 
390  // exracting parameters
391  NNcc=(double)signal.netcc[1];
392  NNrho=(double)signal.rho[0];
393  Mc=(double)signal.chirp[iMc];
394 
395  // evaluating the event
396  for (int i=0; i<nINP; i++){
397  params[i]=x[i];
398  }
399  for (int u=0;u<nNN;u++) {
400  cout<<" u "<<u<<" nNN "<<nNN<<endl;
401  double output=mlp[u]->Evaluate(0,params);
402  cout<<output<<endl;
403  out[u]=output;
404  }
405 
406 
407  cout<<"Dopo param"<<endl;
408 
409 //calculating the average and the standard deviation
410 
411  //wrong!!!
412  /*for (int u=0;u<nNN;u++){
413  if((u!=indNN&&ni==0)||ni!=0) {average=out[u]+average;std+=out[u]*out[u];if(out[u]<0.6) FD0[u]=FD0[u]+1;}
414  }*/
415  for (int u=0;u<nNN;u++){
416  if((u!=indNN&&ni==0)||uf==0) {average=out[u]+average;std+=out[u]*out[u];if(out[u]<0.6) FD0[u]=FD0[u]+1;}
417  }
418 
419 
420  if(nNN!=1&&ni==0) {average=average/(nNN-countNN); if((nNN-1-countNN)!=0)std=pow((std/(nNN-1-countNN)-average*average),0.5);else std=pow((std/(nNN-countNN)-average*average),0.5);}
421 
422 
423  //wrong!!!
424  //if(nNN!=1&&ni!=0) {average=average/(nNN); if((nNN-1)!=0)std=pow((std/(nNN-1)-average*average),0.5);else std=pow((std/nNN-average*average),0.5);}
425  if(nNN!=1&&uf==0) {average=average/(nNN); if((nNN-1)!=0)std=pow((std/(nNN-1)-average*average),0.5);else std=pow((std/nNN-average*average),0.5);}
426  cout<<"average: "<<average<<endl;
427  cout<<"Mc "<<Mc<<endl;
428  if (average<0.6) FD=FD+1;
429  NNTree2->Fill();
430  S_eff=S_eff+1;
431 }
432 //closing the cycle on SIG events
433 
434 //defining BKG parameters
435 int B_eff=0;
436 int FA=0;
437 cout<<"riempito sig S_eff"<<S_eff<<endl;
438 int bg_05=0;
439 int cont_su=0;
440 int cont_su5=0;
441 int cont_su4=0;
442 int cont_su3=0;
443 
444 //opening the cycle on BKG events
445  for(int n=b;n<b+TB;n++) {
446  average=0.;
447  int indNN=-1;
448  int countNN=0;
449  for(int u=0;u<nNN;u++){
450  cout<<" uf "<<uf<<" n "<<n<<" u "<<u<<" NN b[u] "<<NNb[u]<<" NNb[u]+NNnTB[u]"<<NNb[u]+NNnTB[u]<<endl;
451  if (uf!=0&&n>=NNb[u]&&n<(NNb[u]+NNnTB[u])) { indNN=u;countNN=countNN+1; }
452  }
453 
454  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
455  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
456  if(ni==0&&countNN==0){ cout<<"countNN"<<countNN<<endl; continue;}
457 
458  NNTree->GetEntry(n);
459  cout<<"BKG->n: "<<n<<"Bg index"<<(n-sig_entries)<<endl;
460  background.GetEntry(n-sig_entries);
461  NNcc=(double)background.netcc[1];
462  NNrho=(double)background.rho[0];
463  Mc=(double)background.chirp[iMc];
464  for (int i=0; i<nINP; i++){
465  params[i]=x[i];
466  }
467 
468  for(int u=0;u<nNN;u++) {
469  double output=mlp[u]->Evaluate(0,params);
470  cout<<output<<endl;
471  out[u]=output;
472  }
473  int nnsu=0;
474  for (int u=0;u<nNN;u++){
475  if((u!=indNN&&ni==0)||ni!=0) { cout<<" u "<<u<<" n "<<n<<" FA0[u] "<<FA0[u]<<endl;average=out[u]+average;std=out[u]*out[u];TotBg[u]=TotBg[u]+1; if(out[u]>0.6) {FA0[u]=FA0[u]+1;nnsu=nnsu+1;}}
476  }
477  if(nnsu==6) cont_su=cont_su+1;
478  if(nnsu>4) cont_su5=cont_su5+1;
479  if(nnsu>3) cont_su4=cont_su4+1;
480  if(nnsu>2) cont_su3=cont_su3+1;
481 
482 //calculating the final average value
483  if(nNN!=1&&ni==0) {cout<<average<<endl;average=average/(nNN-countNN);
484  cout<<"ni==0"<<average<<endl; if((nNN-1-countNN)!=0)std=pow((std/(nNN-1-countNN)-average*average),0.5);else std=pow((std/(nNN-countNN)-average*average),0.5);}
485 
486  //if(nNN!=1&&ni!=0) {cout<<average<<endl;average=average/(nNN);if((nNN-1)!=0)std=pow((std/(nNN-1)-average*average),0.5);else std=pow((std/(nNN)-average*average)0,5);cout<<"ni!=0"<<average<<endl;}
487  //wrong!!!
488  //if(nNN!=1&&ni!=0) {cout<<average<<endl;average=average/(nNN);if((nNN-1)!=0)std=pow((std/(nNN-1)-average*average),0.5);else std=pow((std/nNN-average*average),0.5); cout<<"ni!=0"<<average<<endl;}
489 
490  if(nNN!=1&&uf==0) {cout<<average<<endl;average=average/(nNN);if((nNN-1)!=0)std=pow((std/(nNN-1)-average*average),0.5);else std=pow((std/nNN-average*average),0.5); cout<<"uf==0"<<average<<endl;}
491  cout<<"average: "<<average<<endl;
492 
493 //saving false alarms
494  if(average>0.6) FA=FA+1;
495 
496  NNTree2->Fill();
497  B_eff=B_eff+1;
498  }
499 cout<<"bkg filled"<<endl;
500 
501 NNTree2->Write();
502 f->Close();
503 
504 cout<<"closed file"<<endl;
505 cout<<ofile<<endl;
506 
507 int cont=0;
508 cout<<" S_eff "<<S_eff<<" B_eff "<<B_eff<<endl;
509 
510 double freq_c[nNN];
511 double meanf=0.;
512 double dev=0.;
513 
514 for (int yy=0; yy<nNN; yy++){
515  freq_c[yy]=0.;
516  cout<<" FA0[yy] "<<FA0[yy]<<" FD0[yy] "<<FD0[yy]<<" FD "<<FD<<" TotBg[yy] "<<TotBg[yy]<<endl;
517  if(TotBg[yy]!=0){freq_c[yy]=(double)FA0[yy]/TotBg[yy];cout<<" freq_c[yy] "<<freq_c[yy]<<endl;}
518  if(freq_c[yy]!=0) {meanf=freq_c[yy]+meanf;cont=cont+1;cout<<" cont "<<cont<<endl;}
519  dev=freq_c[yy]*freq_c[yy]+dev;
520 }
521 
522 if(cont==0) cout<<"Error cont==0"<<endl;
523 if(cont!=0) meanf=meanf/cont;
524 dev=0.;
525 for (int yy=0; yy<nNN; yy++){
526 cout<<" dev "<<dev<<endl;
527 cout<<"freq_c[yy] "<<freq_c[yy]<<" meanf "<<meanf<<endl;
528 if(freq_c[yy]!=0)dev+=(freq_c[yy]-meanf)*(freq_c[yy]-meanf);
529 }
530 Annth(ofile);
531 graph(ofile);
532 PlotsAv_cc(ofile);
533 PlotsAv_Mc(ofile);
534 if(cont!=0) dev=pow(dev/(cont-1),0.5);
535 cout<<"meanf "<<meanf<<" dev "<<dev<<endl;
536 cout<<" FA average "<<FA<<endl;
537 cout<<" cont_su "<<cont_su<<" cont_su5 "<<cont_su5<<" cont_su4 "<<cont_su4<<" cont_su3 "<<cont_su3<<endl;
538 
539 cout<<" FA0[1] "<<FA0[1]<<" TotBg[1] "<<TotBg[1]<<endl;
540 cout<<" FD0[1] "<<FD0[1]<<" TotBg[1] "<<TotBg[1]<<endl;
541 }
542 
543 void graph(TString ifile){
544  TString name(ifile);
545  name.ReplaceAll("average_file/","");
546  TFile* fTEST =TFile::Open(ifile.Data());
547  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
548  double av;
549  double cc;
550  double rho;
551  int nSi;
552  NNTree2->SetBranchAddress("Average",&av);
553  NNTree2->SetBranchAddress("cc",&cc);
554  NNTree2->SetBranchAddress("rho",&rho);
555  NNTree2->SetBranchAddress("#TestSig",&nSi);
556  NNTree2->GetEntry(0);
557  int const nSig=nSi;
558  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
559  cout<<" BG: "<<NNTree2->GetEntries()-nSi<<endl;
560  int const ncurve=nANN*nCC;
561  cout<<"dentro funzione dopodef"<<endl;
562 //LOG(NBg)vs RHO------------------------------------------
563  double* rhoSig[ncurve];
564  for (int i=0;i<ncurve;i++) rhoSig[i]=new double[nSig];
565  cout<<"dopo def rhoSig"<<endl;
566  //double rhoSig[20][10];
567  int NSig[ncurve];
568 // int nBg=0;
569  int const nBg=NNTree2->GetEntries()-nSig;
570  for (int i=0;i<ncurve;i++) {
571  NSig[i]=0;
572  for (int j=0;j<nSig;j++) rhoSig[i][j]=0.;
573  }
574  cout<<"dopo def rhoSig"<<endl;
575  double* rhoBg[ncurve];
576  for (int i=0;i<ncurve;i++) rhoBg[i]=new double[nBg];
577  //double rhoBg[20][10];
578  int NBg[ncurve];
579  for (int i=0;i<ncurve;i++) {
580  NBg[i]=0;
581  for (int j=0;j<nBg;j++) rhoBg[i][j]=0.;
582  }
583  cout<<"dopo def rhoBg"<<endl;
584  double ccTh[nCC];
585  for (int i=0;i<nCC;i++) ccTh[i]=0.;
586  double NNTh[nANN];
587  for (int i=0;i<nANN;i++) NNTh[i]=0.;
588  double deltacc=0.;
589 // double deltaANN=0.;
590  deltacc=0.2/nCC;
591  //deltaANN=0.6/(nANN-1);
592 
593  cout<<NNTree2->GetEntries()<<endl;
594  for(int n=0;n<NNTree2->GetEntries();n++){
595  cout<<n<<endl;
596  NNTree2->GetEntry(n);
597  cout<<"rho "<<rho<<" cc "<<cc<<" av "<<av<<endl;
598  for(int i=0; i<nCC;i++){
599  ccTh[i]=0.5+i*deltacc;
600  if(cc<ccTh[i]) continue;
601  for(int j=0; j<nANN;j++){
602  if(j==0) NNTh[j]=-1000.;
603  //else NNTh[j]=0.5+(j-1)*deltaANN;
604  //else NNTh[j]=0.1+(j-1)*deltaANN;
605  else NNTh[j]=0.+(j)*deltaANN;
606  //else NNTh[j]=0.+(j-1)*deltaANN/1000.;
607  //else NNTh[j]=1.;
608  if(av<NNTh[j]) continue;
609  int ni=0;
610  if(n>nSig) {
611 
612  NBg[i*nANN+j]= NBg[i*nANN+j]+1;
613  while(rhoBg[i*nANN+j][ni]!=0)ni=ni+1;
614  rhoBg[i*nANN+j][ni]=rho;
615  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
616 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
617  }
618  else {
619  NSig[i*nANN+j]= NSig[i*nANN+j]+1;
620  while(rhoSig[i*nANN+j][ni]!=0)ni=ni+1;
621  rhoSig[i*nANN+j][ni]=rho;
622  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
623 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
624  }
625  // }
626  }
627  }
628  }
629  cout<<"dopo riempimento variabili"<<endl;
630  int* indexS[ncurve];
631  for (int i=0;i<ncurve;i++) indexS[i]=new int[nSig];
632 
633  TGraph * gS[ncurve];
634  for (int y=0;y<ncurve;y++) {
635  int igS=0;
636  int igS_p=0;
637  /* int indexS[nSig];
638  double rhoS[nSig];
639  for(int i=0;i<nSig;i++) {
640  rhoS[i]=0.;
641  indexS[i]=0;
642  }*/
643  // ig=1;
644  // for(int i=0;i<nSig;i++) rhoS[i]=rhoSig[y][i];
645  gS[y]=new TGraph();
646 // gS[y]->SetMarkerStyle(7);
647  TMath::Sort(nSig,rhoSig[y],indexS[y],false);
648 // cout<<"dopo Sort "<<y<<endl;
649  for (int k=0;k<nSig;k++) {
650  int ii=indexS[y][k];
651  int yy=0;
652  if (k>0){
653 // cout<<"k "<<k<<endl;
654  int ij=indexS[y][k-1];
655  //if(rhoSig[y][ii]!=0&&rhoSig[y][ii]!=rhoSig[y][ij]) {
656  if(rhoSig[y][ii]!=0) {
657  yy=NSig[y]-igS;
658  //gS[y]->SetPoint(igS,rhoSig[y][ii],yy);
659  if(rhoSig[y][ii]!=rhoSig[y][ij]) gS[y]->SetPoint(igS_p++,rhoSig[y][ii],yy);
660  cout<<"igS"<<igS<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
661  igS=igS+1;
662  }
663  }
664  else {
665  if(rhoSig[y][ii]!=0){
666  yy=NSig[y]-igS;
667  gS[y]->SetPoint(0,rhoSig[y][ii],yy);
668  igS=igS+1;
669 // cout<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
670  }
671 
672  }
673  }
674  }
675  // cout<<"dopo inserimento puntiiS"<<endl;
676 // cout<<ncurve<<endl;
677  int* indexB[ncurve];
678  for (int i=0;i<ncurve;i++) indexB[i]=new int[nBg];
679 
680  TGraph * gB[ncurve];
681  for (int y=0;y<ncurve;y++) {
682  int igB=0;
683  int igB_p=0;
684  gB[y]=new TGraph();
685  TMath::Sort(nBg,rhoBg[y],indexB[y],false);
686 // cout<<"dopo Sort "<<y<<endl;
687  for (int k=0;k<nBg;k++) {
688  int ii=indexB[y][k];
689  int yy=0;
690  if (k>0){
691  int ij=indexB[y][k-1];
692  //if(rhoBg[y][ii]!=0&&rhoBg[y][ii]!=rhoBg[y][ij]) {
693  if(rhoBg[y][ii]!=0) {
694  yy=NBg[y]-igB;
695  //gB[y]->SetPoint(igB,rhoBg[y][ii],yy);
696  if(rhoBg[y][ii]!=rhoBg[y][ij]) gB[y]->SetPoint(igB_p++,rhoBg[y][ii],yy);
697  igB=igB+1;
698 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
699  }
700  }
701  else {
702  if(rhoBg[y][ii]!=0) {
703  yy=NBg[y]-igB;
704  gB[y]->SetPoint(0,rhoBg[y][ii],yy);
705  igB=igB+1;
706 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
707  }
708  }
709  }
710  }
711  cout<<"dopo inserimento puntiB"<<endl;
712  //gB[1]->SetMarkerStyle(7);
713  //gB[1]->SetPoint(1,1,2);
714  TCanvas* cS=new TCanvas("Efficiency_vs_rho","Efficiency_vs_rho",0,0,1200,700);
715  cS->Divide(2,2);
716  cS->cd(1)->SetLogy();
717  TMultiGraph* mg1=new TMultiGraph();
718 // gS[0]->SetMarkerStyle(7);
719  gS[0]->SetMarkerColor(2);
720  gS[0]->SetLineColor(2);
721  mg1->SetTitle("cc=0.5;rho;#Events");
722  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
723 // gS[0]->Draw("apl");
724  for (int h=1;h<nANN;h++){
725  gS[h]->SetMarkerColor(3);
726  gS[h]->SetLineColor(3);
727 // gS[h]->SetMarkerStyle(h+1);
728  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
729  // gS[h]->Draw("apl,same");
730  }
731  mg1->Draw("apl");
732  cS->cd(2)->SetLogy();
733  TMultiGraph* mg2=new TMultiGraph();
734 // gS[nANN]->SetMarkerStyle(7);
735  gS[nANN]->SetMarkerColor(2);
736  gS[nANN]->SetLineColor(2);
737  mg2->SetTitle("cc=0.55;rho;#Events");
738  if(gS[nANN]->GetN()!=0) mg2->Add(gS[nANN]);
739  for (int h=1;h<nANN;h++){
740  gS[nANN+h]->SetMarkerColor(3);
741  gS[nANN+h]->SetLineColor(3);
742 // gS[nANN+h]->SetMarkerStyle(h+1);
743  if(gS[nANN+h]->GetN()!=0) mg2->Add(gS[nANN+h]);
744  // gS[nANN+h]->Draw("apl,same");
745  }
746  mg2->Draw("apl");
747  cS->cd(3)->SetLogy();
748  TMultiGraph* mg3=new TMultiGraph();
749 // gS[nANN*2]->SetMarkerStyle(7);
750  gS[nANN*2]->SetMarkerColor(2);
751  gS[nANN*2]->SetLineColor(2);
752  mg3->SetTitle("cc=0.6;rho;#Events");
753  if(gS[nANN*2]->GetN()!=0) mg3->Add(gS[nANN*2]);
754  for (int h=1;h<nANN;h++){
755  gS[2*nANN+h]->SetMarkerColor(3);
756  gS[2*nANN+h]->SetLineColor(3);
757 // gS[2*nANN+h]->SetMarkerStyle(h+1);
758  if(gS[2*nANN+h]->GetN()!=0) mg3->Add(gS[2*nANN+h]);
759  // gS[2*nANN+h]->Draw("apl,same");
760  }
761  mg3->Draw("apl");
762  cS->cd(4)->SetLogy();
763  TMultiGraph* mg4=new TMultiGraph();
764 // gS[nANN*3]->SetMarkerStyle(7);
765  gS[nANN*3]->SetMarkerColor(2);
766  gS[nANN*3]->SetLineColor(2);
767  mg4->SetTitle("cc=0.65;rho;#Events");
768  if(gS[nANN*3]->GetN()!=0) mg4->Add(gS[nANN*3]);
769  // gS[nANN*3]->Draw("apl");
770  for (int h=1;h<nANN;h++){
771  gS[3*nANN+h]->SetMarkerColor(3);
772  gS[3*nANN+h]->SetLineColor(3);
773 // gS[3*nANN+h]->SetMarkerStyle(h+1);
774  if(gS[3*nANN+h]->GetN()!=0) mg4->Add(gS[3*nANN+h]);
775  // gS[3*nANN+h]->Draw("apl,same");
776  }
777  mg4->Draw("apl");
778  cout<<"nuovo canv"<<endl;
779  TCanvas* cB=new TCanvas("Number_vs_rho","Number_vs_rho",0,0,1200,700);
780  cB->Divide(2,2);
781  cB->cd(1)->SetLogy();
782  TMultiGraph* mg1B=new TMultiGraph();
783 // gB[0]->SetMarkerStyle(7);
784  gB[0]->SetMarkerColor(2);
785  gB[0]->SetLineColor(2);
786  mg1B->SetTitle("cc=0.5;rho;#Events");
787  if(gB[0]->GetN()!=0) mg1B->Add(gB[0]);
788  for (int h=1;h<nANN;h++){
789  gB[h]->SetMarkerColor(3);
790  gB[h]->SetLineColor(3);
791  // gB[h]>SetMarkerStyle(h+1);
792  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
793  //gB[h]->Draw("apl,same");
794  }
795  mg1B->Draw("apl");
796  cB->cd(2)->SetLogy();
797  TMultiGraph* mg2B=new TMultiGraph();
798  //gB[nANN]->SetMarkerStyle(7);
799  gB[nANN]->SetMarkerColor(2);
800  gB[nANN]->SetLineColor(2);
801  mg2B->SetTitle("cc=0.55;rho;#Events");
802  if(gB[nANN]->GetN()!=0) mg2B->Add(gB[nANN]);
803  for (int h=1;h<nANN;h++){
804  gB[nANN+h]->SetMarkerColor(3);
805  gB[nANN+h]->SetLineColor(3);
806  //gB[nANN+h]>SetMarkerStyle(h+1);
807  if(gB[nANN+h]->GetN()!=0) mg2B->Add(gB[nANN+h]);
808  //gB[h]->Draw("apl,same");
809  }
810  mg2B->Draw("apl");
811  cB->cd(3)->SetLogy();
812  TMultiGraph* mg3B=new TMultiGraph();
813 // gB[2*nANN]->SetMarkerStyle(7);
814  gB[2*nANN]->SetMarkerColor(2);
815  gB[2*nANN]->SetLineColor(2);
816  mg3B->SetTitle("cc=0.6;rho;#Events");
817  if(gB[2*nANN]->GetN()!=0) mg3B->Add(gB[2*nANN]);
818  for (int h=1;h<nANN;h++){
819  gB[2*nANN+h]->SetMarkerColor(3);
820  gB[2*nANN+h]->SetLineColor(3);
821 // gB[2*nANN+h]>SetMarkerStyle(h+1);
822  if(gB[2*nANN+h]->GetN()!=0) mg3B->Add(gB[2*nANN+h]);
823  //gB[h]->Draw("apl,same");
824  }
825  mg3B->Draw("apl");
826  cB->cd(4)->SetLogy();
827  TMultiGraph* mg4B=new TMultiGraph();
828 // gB[3*nANN]->SetMarkerStyle(7);
829  gB[3*nANN]->SetMarkerColor(2);
830  gB[3*nANN]->SetLineColor(2);
831  mg4B->SetTitle("cc=0.65;rho;#Events");
832  if(gB[3*nANN]->GetN()!=0) mg4B->Add(gB[3*nANN]);
833  for (int h=1;h<nANN;h++){
834  gB[3*nANN+h]->SetMarkerColor(3);
835  gB[3*nANN+h]->SetLineColor(3);
836 // gB[3*nANN+h]>SetMarkerStyle(h+1);
837  if(gB[3*nANN+h]->GetN()!=0) mg4B->Add(gB[3*nANN+h]);
838  //gB[h]->Draw("apl,same");
839  }
840  mg4B->Draw("apl");
841 
842 // cout<<"dopo Draw()"<<endl;
843  TString CnameS(name);
844  TString CnameB(name);
845  TString CnameSroot(name);
846  TString CnameBroot(name);
847  char CnameS2[1024];
848  char CnameB2[1024];
849  char CnameS2root[1024];
850  char CnameB2root[1024];
851  CnameS.ReplaceAll(".root",".png");
852  CnameB.ReplaceAll(".root",".png");
853  sprintf(CnameS2,"logN_rho_av/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameS.Data());
854  sprintf(CnameB2,"logN_rho_av/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameB.Data());
855  sprintf(CnameS2root,"logN_rho_av/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameSroot.Data());
856  sprintf(CnameB2root,"logN_rho_av/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameBroot.Data());
857  cS->SaveAs(CnameS2);
858  cB->SaveAs(CnameB2);
859  cS->Print(CnameS2root);
860  cB->Print(CnameB2root);
861 // cout<<"fine"<<endl;
862 
863 }
864 
865 
866 void Annth(TString ifile){
867  TString name(ifile);
868  name.ReplaceAll("average_file/","");
869  TFile* fTEST =TFile::Open(ifile.Data());
870  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
871  double av;
872  double cc;
873  double rho;
874  int nSi;
875  NNTree2->SetBranchAddress("Average",&av);
876  NNTree2->SetBranchAddress("cc",&cc);
877  NNTree2->SetBranchAddress("rho",&rho);
878  NNTree2->SetBranchAddress("#TestSig",&nSi);
879  NNTree2->GetEntry(0);
880  int const nSig=nSi;
881 // cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
882  int const ncurve2=nRHO*nCC;
883 
884 
885  double* ANNSig[ncurve2];
886  for (int i=0;i<ncurve2;i++) ANNSig[i]=new double[nSig];
887  //double rhoSig[20][10];
888  int NSig[ncurve2];
889 // int nBg=0;
890  int const nBg=NNTree2->GetEntries()-nSig;
891  for (int i=0;i<ncurve2;i++) {
892  NSig[i]=0;
893  for (int j=0;j<nSig;j++) ANNSig[i][j]=0.;
894  }
895  double* ANNBg[ncurve2];
896  for (int i=0;i<ncurve2;i++) ANNBg[i]=new double[nBg];
897 
898  //double rhoBg[20][10];
899  int NBg[ncurve2];
900  for (int i=0;i<ncurve2;i++) {
901  NBg[i]=0;
902  for (int j=0;j<nBg;j++) ANNBg[i][j]=0.;
903  }
904  double ccTh[nCC];
905  for (int i=0;i<nCC;i++) ccTh[i]=0.;
906  double rhoTh[nRHO];
907  for (int i=0;i<nRHO;i++) rhoTh[i]=0.;
908  double deltacc=0.;
909  double deltarho=0.;
910  deltacc=0.2/nCC;
911  deltarho=1./(nRHO);
912 
913 
914  for(int n=0;n<NNTree2->GetEntries();n++){
915  NNTree2->GetEntry(n);
916  cout<<"rho "<<rho<<" cc "<<cc<<" av "<<av<<endl;
917  for(int i=0; i<nCC;i++){
918  ccTh[i]=0.5+i*deltacc;
919  if(cc<ccTh[i]) continue;
920  for(int j=0; j<nRHO;j++){
921  rhoTh[j]=5+j*deltarho;
922  if(rho<rhoTh[j]) continue;
923  int ni=0;
924  if(n>nSig) {
925  NBg[i*nRHO+j]= NBg[i*nRHO+j]+1;
926  while(ANNBg[i*nRHO+j][ni]!=0)ni=ni+1;
927  ANNBg[i*nRHO+j][ni]=av;
928  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
929  }
930  else {
931  NSig[i*nRHO+j]= NSig[i*nRHO+j]+1;
932  while(ANNSig[i*nRHO+j][ni]!=0)ni=ni+1;
933  ANNSig[i*nRHO+j][ni]=av;
934  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
935  }
936  }
937  }
938  }
939 
940  int* indexS[ncurve2];
941  for (int i=0;i<ncurve2;i++) indexS[i]=new int[nSig];
942 
943  TGraph * gS[ncurve2];
944  for (int y=0;y<ncurve2;y++) {
945  int igS=0;
946  int igS_p=0;
947  gS[y]=new TGraph();
948  TMath::Sort(nSig,ANNSig[y],indexS[y],false);
949  for (int k=0;k<nSig;k++) {
950  int ii=indexS[y][k];
951  int yy=0;
952  if (k>0){
953  //cout<<"k "<<k<<endl;
954  int ij=indexS[y][k-1];
955  //if(ANNSig[y][ii]!=0&&ANNSig[y][ii]!=ANNSig[y][ij]) {
956  if(ANNSig[y][ii]!=0) {
957  yy=NSig[y]-igS;
958  //gS[y]->SetPoint(igS,ANNSig[y][ii],yy);
959  if(ANNSig[y][ii]!=ANNSig[y][ij]) gS[y]->SetPoint(igS_p++,ANNSig[y][ii],yy);
960  // cout<<"igS"<<igS<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
961  igS=igS+1;
962 
963  }
964  }
965  else {
966  if(ANNSig[y][ii]!=0){
967  yy=NSig[y]-igS;
968  gS[y]->SetPoint(0,ANNSig[y][ii],yy);
969  igS=igS+1;
970  // cout<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
971  }
972 
973  }
974  }
975  }
976 
977 
978 
979  int* indexB[ncurve2];
980  for (int i=0;i<ncurve2;i++) indexB[i]=new int[nBg];
981 
982  TGraph * gB[ncurve2];
983  for (int y=0;y<ncurve2;y++) {
984  int igB=0;
985  int igB_p=0;
986  gB[y]=new TGraph();
987  TMath::Sort(nBg,ANNBg[y],indexB[y],false);
988  // cout<<"dopo Sort "<<y<<endl;
989  for (int k=0;k<nBg;k++) {
990  int ii=indexB[y][k];
991  int yy=0;
992  if (k>0){
993  int ij=indexB[y][k-1];
994  //if(ANNBg[y][ii]!=0&&ANNBg[y][ii]!=ANNBg[y][ij]) {
995  if(ANNBg[y][ii]!=0) {
996  yy=NBg[y]-igB;
997  //gB[y]->SetPoint(igB,ANNBg[y][ii],yy);
998  if(ANNBg[y][ii]!=ANNBg[y][ij]) gB[y]->SetPoint(igB_p++,ANNBg[y][ii],yy);
999  igB=igB+1;
1000  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1001  }
1002  }
1003  else {
1004  if(ANNBg[y][ii]!=0) {
1005  yy=NBg[y]-igB;
1006  gB[y]->SetPoint(0,ANNBg[y][ii],yy);
1007  igB=igB+1;
1008  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1009  }
1010  }
1011  }
1012  }
1013 
1014 
1015  TCanvas* cS=new TCanvas("Efficiency_vs_ANN","Efficiency_vs_ANN",0,0,1200,700);
1016  cS->Divide(2,2);
1017  //cS->cd(1)->SetLogy();
1018  cS->cd(1);
1019  TMultiGraph* mg1=new TMultiGraph();
1020  mg1->SetTitle("cc=0.5;ANN;#Events");
1021  for (int h=0;h<nRHO;h++){
1022  gS[h]->SetLineColor(4);
1023  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
1024  }
1025  mg1->Draw("al");
1026  cS->cd(2);
1027  // cS->cd(2)->SetLogy();
1028  TMultiGraph* mg2=new TMultiGraph();
1029  mg2->SetTitle("cc=0.55;ANN;#Events");
1030  for (int h=0;h<nRHO;h++){
1031  gS[nRHO+h]->SetLineColor(4);
1032  if(gS[nRHO+h]->GetN()!=0) mg2->Add(gS[nRHO+h]);
1033  }
1034  mg2->Draw("al");
1035 
1036  //cS->cd(3)->SetLogy();
1037  cS->cd(3);
1038  TMultiGraph* mg3=new TMultiGraph();
1039  mg3->SetTitle("cc=0.6;ANN;#Events");
1040  for (int h=0;h<nRHO;h++){
1041  gS[2*nRHO+h]->SetLineColor(4);
1042  if(gS[2*nRHO+h]->GetN()!=0) mg3->Add(gS[2*nRHO+h]);
1043  }
1044  mg3->Draw("al");
1045  //cS->cd(4)->SetLogy();
1046  cS->cd(4);
1047  TMultiGraph* mg4=new TMultiGraph();
1048  mg4->SetTitle("cc=0.65;ANN;#Events");
1049  for (int h=0;h<nRHO;h++){
1050  gS[3*nRHO+h]->SetLineColor(4);
1051  if(gS[3*nRHO+h]->GetN()!=0) mg4->Add(gS[3*nRHO+h]);
1052  }
1053  mg4->Draw("al");
1054 
1055  TCanvas* cB=new TCanvas("Number_vs_ANN","Number_vs_ANN",0,0,1200,700);
1056  cB->Divide(2,2);
1057  cB->cd(1)->SetLogy();
1058  TMultiGraph* mg1B=new TMultiGraph();
1059  mg1B->SetTitle("cc=0.5;ANN;#Events");
1060  for (int h=0;h<nRHO;h++){
1061  gB[h]->SetLineColor(4);
1062  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
1063  }
1064  mg1B->Draw("al");
1065  cB->cd(2)->SetLogy();
1066  TMultiGraph* mg2B=new TMultiGraph();
1067  mg2B->SetTitle("cc=0.55;ANN;#Events");
1068  for (int h=0;h<nRHO;h++){
1069  gB[nRHO+h]->SetLineColor(4);
1070  if(gB[nRHO+h]->GetN()!=0) mg2B->Add(gB[nRHO+h]);
1071  }
1072  mg2B->Draw("al");
1073  cB->cd(3)->SetLogy();
1074  TMultiGraph* mg3B=new TMultiGraph();
1075  mg3B->SetTitle("cc=0.6;ANN;#Events");
1076  for (int h=0;h<nRHO;h++){
1077  gB[2*nRHO+h]->SetLineColor(4);
1078  if(gB[2*nRHO+h]->GetN()!=0) mg3B->Add(gB[2*nRHO+h]);
1079  }
1080  mg3B->Draw("al");
1081  cB->cd(4)->SetLogy();
1082  TMultiGraph* mg4B=new TMultiGraph();
1083  mg4B->SetTitle("cc=0.6;ANN;#Events");
1084  for (int h=0;h<nRHO;h++){
1085  gB[3*nRHO+h]->SetLineColor(4);
1086  if(gB[3*nRHO+h]->GetN()!=0) mg4B->Add(gB[3*nRHO+h]);
1087  }
1088  mg4B->Draw("al");
1089 
1090 
1091  //cout<<"dopo Draw()"<<endl;
1092  TString CnameS(name);
1093  TString CnameB(name);
1094  TString CnameSroot(name);
1095  TString CnameBroot(name);
1096  char CnameS2[1024];
1097  char CnameB2[1024];
1098  char CnameS2root[1024];
1099  char CnameB2root[1024];
1100  CnameS.ReplaceAll(".root",".png");
1101  CnameB.ReplaceAll(".root",".png");
1102  sprintf(CnameS2,"ANNthres_av/N_ANN_S_%s",CnameS.Data());
1103  sprintf(CnameB2,"ANNthres_av/N_ANN_B_%s",CnameB.Data());
1104  sprintf(CnameS2root,"ANNthres_av/N_ANN_S_%s",CnameSroot.Data());
1105  sprintf(CnameB2root,"ANNthres_av/N_ANN_B_%s",CnameBroot.Data());
1106  cS->SaveAs(CnameS2);
1107  cB->SaveAs(CnameB2);
1108  cS->Print(CnameS2root);
1109  cB->Print(CnameB2root);
1110  //cout<<"fine"<<endl;
1111 
1112 //CARTELLA Annth
1113 
1114 
1115 
1116 }
1117 
1118 
1119 void PlotsAv_cc(TString ifile){
1120  TString name(ifile);
1121  name.ReplaceAll("outfile/","");
1122  name.ReplaceAll("average_file/","");
1123  TFile* fTEST =TFile::Open(ifile.Data());
1124  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1125  //double outbis;
1126  double av;
1127  double out;
1128  double cc;
1129  double rho;
1130  int nSi;
1131 // NNTree2->SetBranchAddress("ANNout",&out);
1132  NNTree2->SetBranchAddress("Average",&av);
1133  NNTree2->SetBranchAddress("cc",&cc);
1134  NNTree2->SetBranchAddress("rho",&rho);
1135  NNTree2->SetBranchAddress("#TestSig",&nSi);
1136  NNTree2->GetEntry(0);
1137  int const nSig=nSi;
1138  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1139  int const nBg=NNTree2->GetEntries()-nSig;
1140  cout<<"nBg: "<<nBg<<" Entries: "<<NNTree2->GetEntries()<<endl;
1141 
1142  TGraph * gS[3];
1143  TGraph * gB[3];
1144  gB[0]=new TGraph();
1145  gS[0]=new TGraph();
1146 // cout<<"nSig: "<<nSig<<endl;
1147  for (int n = 0; n <NNTree2->GetEntries(); n++){
1148  NNTree2->GetEntry(n);
1149  if(n<nSig) {
1150  gS[0]->SetPoint(n,cc,av);
1151  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1152  }
1153  else {
1154  gB[0]->SetPoint(n-nSig,cc,av);
1155  cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1156  }
1157  }
1158 
1159 
1160  /*for (int a=0;a<nBg;a++){
1161  if(aRHOB[a]>=RHOth){
1162  for (int b=0;b<nth;b++){
1163  if(aCCB[a]>=cc1th[b]){
1164  for(int c=0;c<nth+1;c++){
1165  if(aANNB[a]>=ANN1th[c]) Z[b*nth+c]=Z[b*nth+c]+1;
1166  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1167  }
1168  }
1169  }
1170  }
1171  }*/
1172  //gS[0]->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1173  //gB[0]->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1174  //gS[0]->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1175  //gB[0]->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1176  gS[0]->SetMarkerColor(2);
1177  gB[0]->SetMarkerColor(4);
1178  gS[0]->SetMarkerStyle(6);
1179  gB[0]->SetMarkerStyle(7);
1180 
1181  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1182  c->Divide(1,2);
1183  c->cd(1);
1184  TMultiGraph* mg1=new TMultiGraph();
1185  mg1->SetTitle("Av_cc");
1186  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1187  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1188  mg1->Draw("ap");
1189  mg1->GetHistogram()->GetXaxis()->SetTitle("cc");
1190  mg1->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1191 // mg1->Draw("ap");
1192  c->cd(2);
1193  TMultiGraph* mg2=new TMultiGraph();
1194  mg2->SetTitle("Av_cc");
1195  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1196  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1197  mg2->Draw("ap");
1198  mg2->GetHistogram()->GetXaxis()->SetTitle("cc");
1199  mg2->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1200 
1201  cout<<" name "<<name<<endl;
1202  TString Cname(name);
1203  TString Cnameroot(name);
1204  char Cname2[1024];
1205  char Cname2root[1024];
1206  Cname.ReplaceAll(".root",".png");
1207  sprintf(Cname2,"average_png/out_Plots_%s",Cname.Data());
1208  sprintf(Cname2root,"average_png/out_Plots_%s",Cnameroot.Data());
1209  c->SaveAs(Cname2);
1210  c->Print(Cname2root);
1211 /*
1212  TCanvas* c2=new TCanvas("Plots_Bkg_on_Sig","Plots_Bkg_on_Sig",0,0,1200,700);
1213  c2->Divide(2,2);
1214  c2->cd(1);
1215  gS[0]->Draw("ap");
1216  gB[0]->Draw("p,same");
1217 
1218  c2->cd(2);
1219  gS[1]->Draw("ap");
1220  gB[1]->Draw("p,same");
1221  c2->cd(3);
1222  gS[2]->Draw("ap");
1223  gB[2]->Draw("p,same");
1224  c2->cd(4);
1225  TText* text2=new TText(0.37,0.0,"no cuts on ANN");
1226  hglitch->SetStats(0);
1227  hglitch->GetXaxis()->SetTitle("cc");
1228  hglitch->GetYaxis()->SetTitle("ANNoutput");
1229  hglitch->GetZaxis()->SetTitle("count");
1230  hglitch->Draw("colz");
1231  text2->Draw();
1232  gPad->SetLogz();
1233 
1234  TString Cname_2(name);
1235  TString Cname_2root(name);
1236  char Cname2_2[1024];
1237  char Cname2_2root[1024];
1238  Cname_2.ReplaceAll(".root",".png");
1239  sprintf(Cname2_2,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2.Data());
1240  sprintf(Cname2_2root,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2root.Data());
1241  c2->SaveAs(Cname2_2);
1242  c2->Print(Cname2_2root);
1243 
1244 //CARTELLA CCi_RHO_ANN_Plots */
1245 }
1246 void PlotsAv_Mc(TString ifile){
1247  TString name(ifile);
1248  name.ReplaceAll("outfile/","");
1249  name.ReplaceAll("average_file/","");
1250  TFile* fTEST =TFile::Open(ifile.Data());
1251  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1252  //double outbis;
1253  double av;
1254  double out;
1255  double cc;
1256  double Mc;
1257  double rho;
1258  int nSi;
1259 // NNTree2->SetBranchAddress("ANNout",&out);
1260  NNTree2->SetBranchAddress("Mchirp",&Mc);
1261  NNTree2->SetBranchAddress("Average",&av);
1262  NNTree2->SetBranchAddress("cc",&cc);
1263  NNTree2->SetBranchAddress("rho",&rho);
1264  NNTree2->SetBranchAddress("#TestSig",&nSi);
1265  NNTree2->GetEntry(0);
1266  int const nSig=nSi;
1267  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1268  int const nBg=NNTree2->GetEntries()-nSig;
1269 
1270  TGraph * gS[3];
1271  TGraph * gB[3];
1272  gB[0]=new TGraph();
1273  gS[0]=new TGraph();
1274 // cout<<"nSig: "<<nSig<<endl;
1275  for (int n = 0; n <NNTree2->GetEntries(); n++){
1276  NNTree2->GetEntry(n);
1277  if(n<nSig) {
1278  gS[0]->SetPoint(n,Mc,av);
1279  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1280  }
1281  else {
1282  gB[0]->SetPoint(n-nSig,Mc,av);
1283  cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1284  }
1285  }
1286 
1287 
1288  /*for (int a=0;a<nBg;a++){
1289  if(aRHOB[a]>=RHOth){
1290  for (int b=0;b<nth;b++){
1291  if(aCCB[a]>=cc1th[b]){
1292  for(int c=0;c<nth+1;c++){
1293  if(aANNB[a]>=ANN1th[c]) Z[b*nth+c]=Z[b*nth+c]+1;
1294  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1295  }
1296  }
1297  }
1298  }
1299  }*/
1300  gS[0]->SetMarkerColor(2);
1301  gB[0]->SetMarkerColor(4);
1302  gS[0]->SetMarkerStyle(6);
1303  gB[0]->SetMarkerStyle(7);
1304 
1305  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1306  c->Divide(1,2);
1307  c->cd(1);
1308  TMultiGraph* mg1=new TMultiGraph();
1309  mg1->SetTitle("Av_Mc");
1310  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1311  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1312  mg1->Draw("ap");
1313  mg1->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1314  mg1->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1315  c->cd(2);
1316  TMultiGraph* mg2=new TMultiGraph();
1317  mg2->SetTitle("Av_Mc");
1318  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1319  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1320  mg2->Draw("ap");
1321  mg2->GetHistogram()->GetXaxis()->SetTitle("Mchirp");
1322  mg2->GetHistogram()->GetYaxis()->SetTitle("Average on ANN ouputs");
1323 
1324  cout<<" name "<<name<<endl;
1325  TString Cname(name);
1326  TString Cnameroot(name);
1327  char Cname2[1024];
1328  char Cname2root[1024];
1329  Cname.ReplaceAll(".root",".png");
1330  sprintf(Cname2,"average_png/Mc_Plots_%s",Cname.Data());
1331  sprintf(Cname2root,"average_png/Mc_Plots_%s",Cnameroot.Data());
1332  c->SaveAs(Cname2);
1333  c->Print(Cname2root);
1334 }
rho
double rho
Definition: cwb_report_prod_2.C:149
graph
void graph(TString ifile)
Definition: Average.C:543
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:93
cwb_online.f
tuple f
Definition: cwb_online.py:91
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
TString
TString("c")
Definition: cwb_report_skymap.C:111
nRHO
#define nRHO
Definition: Average.C:25
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
std
STL namespace.
Average
void Average(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=1)
Definition: Average.C:47
ofile
char ofile[512]
Definition: cwb_condor_rescue.C:40
j
int j
Definition: cwb_net.C:10
deltaANN
#define deltaANN
Definition: Average.C:27
i
i drho i
Definition: cwb_epparameters.C:85
PlotsAv_Mc
void PlotsAv_Mc(TString ifile)
Definition: Average.C:1246
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
nCC
#define nCC
Definition: Average.C:28
h
wavearray< double > h
Definition: Regression_H1.C:25
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
output
char output[256]
Definition: cwb_report_prod_2.C:876
k
int k
Definition: cwb_report_prod_2.C:150
netevent::netcc
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:94
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
PlotsAv_cc
void PlotsAv_cc(TString ifile)
Definition: Average.C:1119
yy
wavearray< double > yy
Definition: TestFrame5.C:12
nANN
#define nANN
Definition: Average.C:26
s
s s
Definition: cwb_net.C:137
Annth
void Annth(TString ifile)
Definition: Average.C:866
deltarho
#define deltarho
Definition: Average_FAout_moreGraphs_3.C:32
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
iMc
#define iMc
Definition: Average.C:34
netevent::GetEntries
Int_t GetEntries()
Definition: netevent.cc:381
for
for(int i=0;i< 101;++i) Cos2[2][i]=0
Definition: CWB_Plugin_MDC_OTF_Config_BRST.C:36
nIFO
#define nIFO
Definition: Average.C:24
y
wavearray< double > y
Definition: Test10.C:31
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
exit
exit(0)