Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Fisher_1.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 3
35 using namespace std;
36 void Fisher_ex(TString ifile, int nFs, int nFb, int Fs, int Fb);
37 void graph_Fish(TString ifile);
38 //void Annth(TString ifile);
39 void PlotsAv_cc(TString ifile);
40 void PlotsAv_Mc(TString ifile);
41 
42 
43 //uf !=0 means the Test file of events is the same used for train the network, ni==0 select the event used for the training, in this way we test the events alway for NNn-1 network, if ni!=0 also the events not used for the training are considered and there the average in made on all the nNN network
44 void Fisher_1(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 nFs, int nFb, int Fs=0, int Fb=0,int TS=0, int TB=0, int s=0, int b=0, int uf=1, int ni=0){
45 if(uf==0) ni=1;
46 int n_NN=0;
47 if (NN_FILE.CompareTo("")){
48  n_NN=n_NN+1;
49  }
50 if (NN_FILE2.CompareTo("")){
51  n_NN=n_NN+1;
52  }
53 if (NN_FILE3.CompareTo("")){
54  n_NN=n_NN+1;
55  }
56 if (NN_FILE4.CompareTo("")){
57  n_NN=n_NN+1;
58  }
59 if (NN_FILE5.CompareTo("")){
60  n_NN=n_NN+1;
61  }
62 if (NN_FILE6.CompareTo("")){
63  n_NN=n_NN+1;
64  }
65 if (NN_FILE7.CompareTo("")){
66  n_NN=n_NN+1;
67  }
68 if (NN_FILE8.CompareTo("")){
69  n_NN=n_NN+1;
70  }
71 int const nNN=n_NN;
72 
73 // int p=0;
74  char NNi[nNN][1024];
75  char NNi2[nNN][1024];
76  sprintf(NNi2[0],"%s",NN_FILE.Data());
77  if(nNN>1) sprintf(NNi2[1],"%s",NN_FILE2.Data());
78  if(nNN>2) sprintf(NNi2[2],"%s",NN_FILE3.Data());
79  if(nNN>3) sprintf(NNi2[3],"%s",NN_FILE4.Data());
80  if(nNN>4) sprintf(NNi2[4],"%s",NN_FILE5.Data());
81  if(nNN>5) sprintf(NNi2[5],"%s",NN_FILE6.Data());
82  if(nNN>6) sprintf(NNi2[6],"%s",NN_FILE7.Data());
83  if(nNN>7) sprintf(NNi2[7],"%s",NN_FILE8.Data());
84  for (int u=0;u<nNN;u++){
85  int p=0;
86  while (NNi2[u][p]){
87  //cout<<NNi2[p]<<endl;
88  if (NNi2[u][p]=='N'){
89 // cout<<p<<endl;
90  if((NNi2[u][p+1]=='N')&&(NNi2[u][p+2]=='\/')) {
91 // cout<<" dentro if"<<endl;
92  int hh=p+3;
93  while (NNi2[u][hh]!='\0'){NNi[u][hh-p-3]=NNi2[u][hh];hh=hh+1;}
94  break;
95  }
96  }
97  p=p+1;
98  }
99 }
100  int q=0;
101  char Filei[1024];
102  char Filei2[1024];
103  sprintf(Filei2,"%s",TEST_FILE.Data());
104  while (Filei2[q]){
105  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]=='\/')) {
106  int hh=q+7;
107  while (Filei2[hh]!='\0') {Filei[hh-q-7]=Filei2[hh];hh=hh+1;
108  }
109  for (int h0=hh-q-7;h0<1024;h0++) Filei[h0]='\0';
110  break;
111  }
112  q=q+1;
113  }
114 
115  cout<<Filei<<" original String: "<<Filei2<<endl;
116 
117 
118  TString NNi0[nNN];
119  NNi0[0]=NNi[0];
120  if(nNN>1) NNi0[1]=NNi[1];
121  if(nNN>2) NNi0[2]=NNi[2];
122  if(nNN>3) NNi0[3]=NNi[3];
123  if(nNN>4) NNi0[4]=NNi[4];
124  if(nNN>5) NNi0[5]=NNi[5];
125  if(nNN>6) NNi0[6]=NNi[6];
126  if(nNN>7) NNi0[7]=NNi[7];
127  TString Filei0(Filei);
128 cout<<"NN"<<endl;
129 int min_NNb=5000000;
130  for (int u=0;u<nNN;u++){
131  NNi0[u].ReplaceAll(".root","");
132  }
133  Filei0.ReplaceAll(".root","");
134 
135  TFile *fnet[nNN];
136  TMultiLayerPerceptron* mlp[nNN];
137  TTree* infot[nNN];
138  int NNs[nNN];
139  int NNb[nNN];
140  int NNnTS[nNN];
141  int NNnTB[nNN];
142  for (int u=0;u<nNN;u++){
143  fnet[u]=TFile::Open(NNi2[u]);
144  cout<<NNi2[u]<<endl;
145  cout<<"dopo file"<<endl;
146  mlp[u]=(TMultiLayerPerceptron*)fnet[u]->Get("TMultiLayerPerceptron");
147 
148  cout<<"dopo TMLP"<<endl;
149  if(mlp[u]==NULL) {cout << "Error getting mlp" << endl;exit(1);}
150  infot[u]=(TTree*)fnet[u]->Get("info");
151  infot[u]->SetBranchAddress("Rand_start_Sig",&NNs[u]);
152  infot[u]->SetBranchAddress("Rand_start_Bg",&NNb[u]);
153  infot[u]->SetBranchAddress("#trainSig",&NNnTS[u]);
154  infot[u]->SetBranchAddress("#trainBg",&NNnTB[u]);
155  infot[u]->GetEntry(0);
156  if(NNb[u]<min_NNb) min_NNb=NNb[u];
157  cout<<"n "<<u<<" b "<<NNb[u]<<" min_NNb "<<min_NNb<<endl;
158 
159  cout<<"da Tree"<<endl;
160  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;
161  }
162 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;}
163 //exit(0);
164 cout<<"Def. tree e mlp"<<endl;
165  TFile* fTEST =TFile::Open(TEST_FILE.Data());
166  TTree* NNTree=(TTree*)fTEST->Get("nnTree");
167  int entries=NNTree->GetEntries();
168  cout<<"entries: "<<entries<<endl;
169 
170 //estraggo le info che servono anche al tree di def dell'mlp
171  int ndim;
172  int ninp;
173  int y;
174  int entriesTot;
175  int bg_entries;
176  int sig_entries;
177 
178  NNTree->SetBranchAddress("#Entries_type",&entriesTot);
179  NNTree->SetBranchAddress("Matrix_dim",&ndim);
180  NNTree->SetBranchAddress("#inputs",&ninp);
181  NNTree->SetBranchAddress("amplitude_mode",&y);
182 
183  NNTree->GetEntry(0);
184  sig_entries=entriesTot;
185  NNTree->GetEntry(entries-1);
186  bg_entries=entriesTot;
187 
188  int const NDIM=ndim;
189  int const nINP=ninp;
190 
191  cout<<"NDIM "<<NDIM<<endl;
192  cout<<"nINP"<<nINP<<endl;
193  cout<<"sig e "<<sig_entries<<endl;
194  cout<<"bg e "<<bg_entries<<endl;
195  int minevents=0;
196  if (sig_entries>bg_entries) minevents=bg_entries;
197  else minevents=sig_entries;
198 
199  if(b==0) b=sig_entries;
200  if(TB==0 && TS==0){
201  TS=minevents;
202  TB=minevents;
203  }
204  if (b<sig_entries){
205  cout<<"Error: Bg index<sig_entries"<<endl;
206  exit(0);
207  }
208  if((TS>sig_entries||TB>bg_entries)&&(TS==TB)) {TS=minevents-s;TB=minevents-b+sig_entries;}
209  if((TS>sig_entries||TB>bg_entries)&&(TS!=TB)) {TS=sig_entries-s;TB=bg_entries-b+sig_entries;}
210 
211 
212  char NOMEtot[1024];
213  sprintf(NOMEtot,"%s",NOMEtot_S.Data());
214  cout<<"nome: "<<NOMEtot<<endl;
215 
216  TString SIG_FILE;
217  TString BG_FILE;
218  char FILE_NAME[516];
219  NNTree->SetBranchAddress("Files_name",&FILE_NAME);
220  NNTree->GetEntry(0);
221  SIG_FILE=FILE_NAME;
222  NNTree->GetEntry(entries-1);
223  BG_FILE=FILE_NAME;
224  cout<<"fine ifdef RHO_CC"<<endl;
225 
226  TChain sigTree("waveburst");//cerca il Tree "waveburst nei file
227  sigTree.Add(SIG_FILE.Data());//determina i file
228  netevent signal(&sigTree,nIFO);
229  int sig_entries2 = signal.GetEntries();
230  cout << "sig entries2 : " << sig_entries2 << endl;
231 
232  TChain bgTree("waveburst");
233  bgTree.Add(BG_FILE.Data());
234  netevent background(&bgTree,nIFO);
235  int bg_entries2 = background.GetEntries();
236  cout << "bg entries2 : " << bg_entries2 << endl;
237 
238  cout<<"b: "<<b<<endl;
239  cout<<"s: "<<s<<endl;
240 
241  // add leaf
242  Float_t x[nINP];
243  for (int jj=0; jj<nINP;jj++) x[jj]=0.;
244  char ilabel[nINP][16];
245 
246  //costruzione una foglia per ogni input
247  for(int i=0;i<nINP;i++) {
248  sprintf(ilabel[i],"x%i",i+1);
249  NNTree->SetBranchAddress(ilabel[i], &x[i]);
250  }
251 
252 char ofile[1024];
253 sprintf(ofile,"average_file/%s.root",NOMEtot);
254 TFile*f =new TFile(ofile,"RECREATE");
255  TTree* NNTree2=new TTree("Parameters","Parameters");
256  NNTree2->SetDirectory(f);
257  double average=0.;
258  double out[nNN];
259  for(int y=0; y<nNN; y++) out[y]=0.;
260  double Mc=0.;
261  double NNcc=0.;
262  double NNrho=0.;
263  double std=0.;
264  NNTree2->Branch("Average",&average,"Average/D");
265  NNTree2->Branch("StandardDevaition",&std,"StandardDeviation/D");
266  NNTree2->Branch("cc",&NNcc,"cc/D");
267  NNTree2->Branch("Mchirp",&Mc,"Mchirp/D");
268  NNTree2->Branch("rho",&NNrho,"rho/D");
269  for(int u=0;u<nNN;u++){
270  char NNoutl[512];
271  char NNoutl2[512];
272  char NNnamel[512];
273  char NNnamel2[512];
274  sprintf(NNoutl,"NNout%i",u);
275  sprintf(NNoutl2,"NNout%i/D",u);
276  sprintf(NNnamel,"NNname%i",u);
277  sprintf(NNnamel2,"NNname%i/C",u);
278  NNTree2->Branch(NNoutl,&out[u],NNoutl2);
279  NNTree2->Branch(NNnamel,&NNi[u],NNnamel2);}
280  char Testf[1024];
281  NNTree2->Branch("TestFile",&Testf,"TestFile/C");
282  sprintf(Testf,"%s",TEST_FILE.Data());
283  int nTestS=0;
284  NNTree2->Branch("#TestSig",&nTestS,"#TestSig/I");
285  cout<<"nTestS: "<<nTestS<<" TS: "<<TS<<endl;
286  cout<<"dopo def tree"<<endl;
287 int scount=0;
288  if(uf==0) nTestS=TS;
289  else {
290 for(int n=s;n<s+TS;n++) {
291  int countNN=0;
292  for(int u=0;u<nNN;u++){
293  if (uf!=0&&n>=NNs[u]&&n<(NNs[u]+NNnTS[u])) countNN=countNN+1;
294  }
295  //if(nNN==1){ni=1;cout<<"out==Averege:only 1NN is considered"<<endl;}
296  if(countNN==1&&ni==0) scount=scount+1;
297  if(countNN<2&&ni!=0) scount=scount+1;
298  if(countNN>1){cout<<"Error: training non independent"; exit(0);}
299  }
300 nTestS=scount;
301 }
302 
303 cout<<"s test "<<s<< " s+TS "<<s+TS<<" b "<<b<<" b+ TB "<<b+TB<<" nTestS "<<nTestS<<endl;
304 //exit(0);
305  double params[nINP];
306  int sig_05=0;
307  for (int i=0; i<nINP; i++) params[i]=0.;
308  //for(int n=0;n<entryS.GetEntries();n++) {
309 
310  int S_eff=0;
311  int TD=0;
312  for(int n=s;n<s+TS;n++) {
313  average=0.;
314  int indNN=-1;
315  int countNN=0;
316  for(int u=0;u<nNN;u++){
317  if (uf!=0&&n>=NNs[u]&&n<(NNs[u]+NNnTS[u])) { indNN=u;countNN=countNN+1;}
318  }
319  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
320  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
321  if(ni==0&&countNN==0)continue;
322  NNTree->GetEntry(n);
323  signal.GetEntry(n);
324  NNcc=(double)signal.netcc[1];
325  NNrho=(double)signal.rho[0];
326  Mc=(double)signal.chirp[iMc];
327  for (int i=0; i<nINP; i++){
328  params[i]=x[i];
329  }
330  for (int u=0;u<nNN;u++) {
331  double output=mlp[u]->Evaluate(0,params);
332  cout<<output<<endl;
333  out[u]=output;
334  }
335  for (int u=0;u<nNN;u++){
336  if((u!=indNN&&ni==0)||ni!=0) {average=out[u]+average;std+=out[u]*out[u];}
337  }
338  //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);}
339  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);}
340  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);}
341  cout<<"average: "<<average<<endl;
342  // cout<<"Mc "<<Mc<<endl;
343  if (average<0.6) TD=TD+1;
344  //if (outbis<0.5) sig_05=sig_05+1;
345  //if (outbis>0.5&&outbis<0.6) if(out<0.1) sig_05=sig_05+1;
346  // if (outbis<0.5 || out<0.5) sig_05=sig_05+1;
347  //if (outbis<0.6) sig_05=sig_05+1;
348  // cout<<"rho: "<<signal.rho[0]<<" cc "<<signal.netcc[1]<<" out: "<<out<<" outbis "<<outbis<<endl;
349  NNTree2->Fill();
350  S_eff=S_eff+1;
351  }
352 int B_eff=0;
353 int FA=0;
354 cout<<"riempito sig"<<endl;
355 int bg_05=0;
356 cout<<"b "<<b<<" TB "<<TB<<endl;
357 //for(int n=b;n<b+TB;n++) {
358 //for(int u=0;u<nNN;u++){
359  // cout<<" u "<<u<<" NN b[u] "<<NNb[u]<<" NNnTB[u] "<<NNnTB[u]<<" NNb[u]+NNnTB[u]"<<NNb[u]+NNnTB[u]<<endl;
360  // }
361 //}
362 //exit(0);
363 cout<<" prima ciclo fro bg"<<endl;
364  //for(int n=sig_entries+b;n<sig_entries+b+TB;n++) {
365  for(int n=b;n<b+TB;n++) {
366  cout<<n<<endl;
367  average=0.;
368  int indNN=-1;
369  int countNN=0;
370  for(int u=0;u<nNN;u++){
371  cout<<" uf "<<uf<<" n "<<n<<" u "<<u<<" NN b[u] "<<NNb[u]<<" (primo non preso) NNb[u]+NNnTB[u]"<<NNb[u]+NNnTB[u]<<endl;
372  if (uf!=0&&n>=NNb[u]&&n<(NNb[u]+NNnTB[u])) { indNN=u;countNN=countNN+1;
373  // cout<<"n "<<n<<" u "<<u<<" NN b[u] "<<NNb[u]<<" (primo non preso)NNb[u]+NNnTB[u]"<<NNb[u]+NNnTB[u]<<endl;
374  }
375  }
376  if (countNN>1){cout<<"Error: training non independent"; exit(0);}
377  if(nNN==1){cout<<"out==Averege:only 1NN is considered"<<endl;}
378  if(ni==0&&countNN==0){
379  cout<<"countNN"<<countNN<<endl;
380  continue;}
381  //if (uf!=0&&(n>=NNb&&n<=(NNb+NNnTB) || n>=NNbbis&&n<=(NNbbis+NNnTBbis))) continue;
382  NNTree->GetEntry(n);
383  cout<<"n: "<<n<<"Bg index"<<(n-sig_entries)<<endl;
384  background.GetEntry(n-sig_entries);
385  NNcc=(double)background.netcc[1];
386  NNrho=(double)background.rho[0];
387  Mc=(double)background.chirp[iMc];
388  for (int i=0; i<nINP; i++){
389  //x[i]=(*frameS)[i];
390  params[i]=x[i];
391  }
392 
393  for(int u=0;u<nNN;u++) {
394  double output=mlp[u]->Evaluate(0,params);
395  cout<<output<<endl;
396  out[u]=output;
397  }
398  for (int u=0;u<nNN;u++){
399  if((u!=indNN&&ni==0)||ni!=0) { average=out[u]+average;std=out[u]*out[u];}
400  }
401  if(nNN!=1&&ni==0) {cout<<average<<endl;average=average/(nNN-countNN);
402 //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);}
403 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);}
404 
405  //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;}
406  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;}
407  cout<<"average: "<<average<<" nNN-1-countNN "<<nNN-1-countNN<<endl;
408  if(average>0.6) FA=FA+1;
409 /* d uble output=mlp->Evaluate(0,params);
410  out=output;
411  double outputbis=mlpbis->Evaluate(0,params);
412  outbis=outputbis;*/
413 // if(outbis>0.6) bg_05=bg_05+1;
414 // if(outbis<=0.6&&outbis>0.5) if(out>0.1)bg_05=bg_05+1;
415 
416 // if(outbis>0.5 && out>0.5) bg_05=bg_05+1;
417  //if(outbis>0.6) bg_05=bg_05+1;
418 
419 // cout<<"rho: "<<background.rho[0]<<" cc "<<background.netcc[1]<<" out: "<<out<<endl;
420  NNTree2->Fill();
421  B_eff=B_eff+1;
422  }
423 cout<<"riempito bg"<<endl;
424 //TFile*f0 =new TFile(ofile,"RECREATE");
425 //NNTree2->SetDirectory(f0);
426 NNTree2->Write();
427 f->Close();
428 //f0->Close();
429 cout<<"chiuso file"<<endl;
430 
431 
432 //graph(ofile.Data());
433 //graph(ofile);
434 cout<<ofile<<endl;
435 cout<<"dopo richiamo funzione"<<endl;
436 Fisher_ex(ofile,nFs,nFb,Fs,Fb);
437 //Annth(ofile);
438 //exit(0);
439 graph_Fish(ofile);
440 //PlotsAv_cc(ofile);
441 //PlotsAv_Mc(ofile);
442 cout<<" S_eff "<<S_eff<<" B_eff "<<B_eff<<endl;
443 cout<<" FA>0.6 "<<FA<<" TD<0.6 "<<TD<<" FA/tot_BG "<<FA/B_eff<<" TD/tot_S "<<TD/S_eff<<endl;
444 //cout<<" Sig with out<06: "<<sig_05<<endl;
445 //cout<<" Bg with double threshold "<<bg_05<<endl;
446 }
447 
448 
449 void Fisher_ex(TString ifile, int nFs, int nFb, int Fs, int Fb){
450 
451 
452 TString name(ifile);
453  name.ReplaceAll("average_file/","");
454  TFile* fTEST =TFile::Open(ifile.Data());
455  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
456  double av;
457  double cc;
458  int nSi;
459  NNTree2->SetBranchAddress("Average",&av);
460  NNTree2->SetBranchAddress("cc",&cc);
461  NNTree2->SetBranchAddress("#TestSig",&nSi);
462  NNTree2->GetEntry(0);
463 
464 //INIZIO DISCRIMINAZIONE FISHER
465 
466 double qS_o=0.;
467 double qS_c=0.;
468 double sS_co=0.;
469 double mS_o=0.;
470 double mS_c=0.;
471 double mB_o=0.;
472 double mB_c=0.;
473 double qB_c=0.;
474 double qB_o=0.;
475 double sB_co=0.;
476 
477 if (nFs+Fs>nSi){ Fs=0; if(nFs>nSi) nFs=nSi; cout<<" change: Fs "<<Fs<<" nFs "<<nFs<<endl;}
478 if(Fb==0)Fb=nSi;
479 if (Fb<nSi) {Fb=nSi+Fb; cout<<" change: Fb "<<Fb<<endl;}
480 
481 for(int n=Fs; n<Fs+nFs;n++){
482  NNTree2->GetEntry(n);
483  mS_o+=av;
484  mS_c+=cc;
485  qS_o+=av*av;
486  qS_c+=cc*cc;
487  sS_co=av*cc;
488 
489 }
490 mS_o=mS_o/nFs;
491 mS_c=mS_c/nFs;
492 
493 int const TOTent=NNTree2->GetEntries();
494 int const nBg=NNTree2->GetEntries()-nSi;
495 cout<<"nSi "<<nSi<<" nBg "<<nBg<<" entries "<<TOTent<<endl;
496 
497 if (nFb+Fb>TOTent){ Fb=nSi; if(nFs>nBg) nFb=nBg; cout<<" change: Fb "<<Fb<<" nFb "<<nFb<<endl;}
498 
499 for(int n=Fb; n<Fb+nFb;n++){
500  NNTree2->GetEntry(n);
501  mB_o+=av;
502  mB_c+=cc;
503  qB_o+=av*av;
504  qB_c+=cc*cc;
505  sB_co=av*cc;
506 
507 }
508 mB_o=mB_o/nFb;
509 mB_c=mB_c/nFb;
510 
511 //ELEMENTS OF WITHIN MATRIX
512 double a_s=0.;
513 double b_s=0.;
514 double d_s=0.;
515 double a_b=0.;
516 double b_b=0.;
517 double d_b=0.;
518 
519 a_b=qB_c-nFb*(mB_c*mB_c);
520 d_b=qB_o-nFb*(mB_o*mB_o);
521 b_b=sB_co-nFb*(mB_c*mB_o);
522 
523 cout<<"a_b "<<a_b<<endl;
524 cout<<"d_b "<<d_b<<endl;
525 cout<<"a_s "<<a_s<<endl;
526 cout<<"d_s "<<d_s<<endl;
527 a_s=qS_c-nFs*(mS_c*mS_c);
528 d_s=qS_o-nFs*(mS_o*mS_o);
529 b_s=sS_co-nFs*(mS_c*mS_o);
530 
531 
532 double a=0.;
533 double b=0.;
534 double d=0.;
535 
536 a=a_b+a_s;
537 b=b_b+b_s;
538 d=d_b+d_s;
539 
540 //ELEMENT OF THE INVERSE WITHIN MATRIX
541 
542 double det=0.;
543 double a_i=0.;
544 double b_i=0.;
545 double d_i=0.;
546 
547 det=a*d-b*b;
548 cout<<"ca standard deviation cc bg "<<pow(a_b/nFb,0.5)<<" ca standard deviation out bg "<<pow(d_b/nFb,0.5)<<" det "<<det<<" media cc bg "<<mB_c<<" media out bg "<<mB_o<<endl;
549 cout<<"ca standard deviation cc sig "<<pow(a_s/nFs,0.5)<<" ca standard deviation out sig "<<pow(d_s/nFs,0.5)<<" det "<<det<<" media cc sig "<<mS_c<<" media out sig "<<mS_o<<endl;
550 if(det==0){cout<<"Problem: det==0"<<endl;exit(0);}
551 a_i=d/det;
552 b_i=-b/det;
553 d_i=a/det;
554 
555 //DIFFERENCES OF THE AVERAGES
556 double dm_o=0.;
557 double dm_c=0.;
558 dm_c=mS_c-mB_c;
559 dm_o=mS_o-mB_o;
560 cout<<" a_i "<<a_i<<" d_i "<<d_i<<" b_i "<<b_i<<" dm_c "<<dm_c<<" dm_o "<<dm_o<<endl;
561 //w1 and w2
562 double w_c=0.;
563 double w_o=0.;
564 w_c=a_i*dm_c+b_i*dm_o;
565 w_o=b_i*dm_c+d_i*dm_o;
566 
567 //FINDING THE THRESHOLD
568 TGraph *ms=new TGraph();
569 TGraph *mb=new TGraph();
570 TGraph *pg=new TGraph();
571 TGraph *g=new TGraph();
572 TGraph *g0=new TGraph();
573 TGraph *g1=new TGraph();
574 int count=0;
575 double w0=0.;
576 for(int n=Fs; n<Fs+nFs;n++){
577  NNTree2->GetEntry(n);
578 // w0*=-(w_c*cc+w_o*av);
579  cout<<" sig"<<n<<" "<<w_c*cc+w_o*av<<endl;
580  g0->SetPoint(count,cc,av);
581  count=count+1;
582 
583 }
584 cout<<"count sig "<<count<<endl;
585 count=0;
586 /*for(int n=Fb; n<Fb+nFb;n++){
587  NNTree2->GetEntry(n);
588  // w0+=-(w_c*cc+w_o*av);
589 count=0;*/
590 for(int n=Fb; n<Fb+nFb;n++){
591  NNTree2->GetEntry(n);
592  // w0+=-(w_c*cc+w_o*av);
593  cout<<" bg"<<n<<" "<<w_c*cc+w_o*av<<endl;
594  g1->SetPoint(count,cc,av);
595  count=count+1;
596 }
597 double mb0=w_c*(mB_c)+w_o*(mB_o);
598 mb0=mb0/2;
599 
600 double ms0=w_c*(mS_c)+w_o*(mS_o);
601 ms0=ms0/2;
602 
603 w0=w_c*(mS_c+mB_c)+w_o*(mS_o+mB_o);
604 cout<<" mS_c+mB_c "<<mS_c+mB_c<<" mS_o+mB_o "<<mS_o+mB_o<<endl;
605 w0=w0/2.;
606 cout<<" count bg "<<count<<endl;
607 double x_ms=0.;
608 double y_ms=0.;
609 double x_mb=0.;
610 double y_mb=0.;
611 double x_thres=0.;
612 double y_thres=0.;
613 /*x_mb=pow(mb0*mb0/(1+(w_c*w_c/(w_o*w_o))),0.5);
614 y_mb=-w_c/w_o*x_mb;
615 x_ms=pow(ms0*ms0/(1+(w_c*w_c/(w_o*w_o))),0.5);
616 y_ms=-w_c/w_o*x_ms;*/
617 x_thres=pow(w0*w0/(1+(w_c*w_c/(w_o*w_o))),0.5);
618 y_thres=-w_c/w_o*x_thres;
619 pg->SetPoint(0,x_thres,y_thres);
620 TF1 retta_ms("retta_p","[1]*x+[0]",-0.5,1.);
621 TF1 retta_mb("retta","[1]*x+[0]",-0.5,1.);
622 TF1 retta_p("retta_p","[1]*x+[0]",-0.5,1.);
623 TF1 retta("retta","[1]*x+[0]",-0.5,1.);
624 retta.SetParameter(0,0.);
625 double c_ang=0.;
626 c_ang=-w_c/w_o;
627 retta.SetParameter(1,c_ang);
628 double c_ang_p=-1./c_ang;
629 retta_p.SetParameter(1,c_ang_p);
630 double q_p=0.;
631 q_p=(mS_o+mB_o)/2.-c_ang_p*(mS_c+mB_c)/2.;
632 //q_p=y_thres-c_ang_p*x_thres;
633 retta_p.SetParameter(0,q_p);
634 retta_ms.SetParameter(1,c_ang_p);
635 double q_ms=0.;
636 //q_ms=y_ms-c_ang_p*x_ms;
637 q_ms=mS_o-c_ang_p*mS_c;
638 retta_ms.SetParameter(0,q_ms);
639 retta_mb.SetParameter(1,c_ang_p);
640 double q_mb=0.;
641 //q_mb=y_mb-c_ang_p*x_mb;
642 q_mb=mB_o-c_ang_p*mB_c;
643 retta_mb.SetParameter(0,q_mb);
644 cout<<"q_mb "<<q_mb<<" q_ms "<<q_ms<<" mS_o "<<mS_o<<" mS_c"<<mS_c<<" mB_o "<<mB_o<<" mB_c"<<mB_c<<endl;
645 cout<<" wc "<<w_c<<" w_o "<<w_o<<" coef ang "<<c_ang<<" x_thres "<<x_thres<<" y_thres "<<y_thres<<" w0 "<<w0<<" q_p "<<q_p<<" c_ang_p"<<c_ang_p<<endl;
646 TCanvas *canv=new TCanvas("reg_wErr_woErr_media","reg_wErr_woErr_media",0,0,600,400);
647 g1->SetMarkerStyle(6);
648 g1->SetMarkerColor(4);
649 g0->SetMarkerStyle(7);
650 g0->SetMarkerColor(2);
651 g->SetPoint(0,1,1.3);
652 g->SetPoint(1,-0.3,-0.3);
653 g->SetMarkerColor(10);
654 g->Draw("ap");
655 g0->Draw("p,same");
656 g1->Draw("p,same");
657 pg->SetMarkerStyle(29);
658 pg->SetMarkerColor(1);
659 retta_ms.SetLineColor(3);
660 retta_mb.SetLineColor(3);
661 retta_p.SetLineColor(6);
662 retta_ms.Draw("same");
663 retta_mb.Draw("same");
664 retta_p.Draw("same");
665 retta.SetLineColor(1);
666 retta.Draw("same");
667 //retta.Draw("al");
668 char nomec[1024];
669 char nomecroot[1024];
670 sprintf(nomec,"Fisher_png/%s_.png",name.Data());
671 sprintf(nomecroot,"Fisher_png/%s_.root",name.Data());
672 canv->SaveAs(nomec);
673 canv->SaveAs(nomecroot);
674 canv->Close();
675 fTEST->Close();
676 
677 TString NOMEtot(ifile);
678 NOMEtot.ReplaceAll("average_file/","");
679 NOMEtot.ReplaceAll(".root","");
680 char namefF[1024];
681 sprintf(namefF,"Fisher_file/%s_Fisher.root",NOMEtot.Data());
682 cout<<" namefF "<<namefF<<" NOMEtot "<<NOMEtot.Data()<<endl;
683 TFile* Fishf =new TFile(namefF,"RECREATE");
684  TTree* TreeFis=new TTree("Fisher_Discriminant","Fisher_Discriminant");
685  TreeFis->SetDirectory(Fishf);
686  double wc=0.;
687  double wo=0.;
688  TreeFis->Branch("Fish_coef_cc",&wc,"Fish_coef_cc/D");
689  TreeFis->Branch("Fish_coef_ANNout",&wo,"Fish_coef_ANNout/D");
690  TreeFis->Branch("nFs",&nFs,"nFs/I");
691  TreeFis->Branch("nFb",&nFb,"nFb/I");
692  TreeFis->Branch("Fb_start",&Fb,"Fb_start/I");
693  TreeFis->Branch("Fs_start",&Fs,"Fs_start/I");
694  char nomeif[2048];
695  sprintf(nomeif,"%s",ifile.Data());
696  TreeFis->Branch("File_name",&nomeif,"File_name/C");
697 
698  wo=w_o;
699  wc=w_c;
700  TreeFis->Fill();
701  cout<<" wc "<<wc<<" wo "<<wo<<endl;
702  cout<<" entries "<<TreeFis->GetEntries()<<endl;
703  TreeFis->Write();
704 
705 Fishf->Close();
706 
707 }
708 
709 
710 void graph_Fish(TString ifile){
711  TString name(ifile);
712  name.ReplaceAll("average_file/","Fisher_file/");
713  TFile* fTEST =TFile::Open(ifile.Data());
714  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
715 
716  double w_o;
717  double w_c;
718  int nFs;
719  int nFb;
720  int Fs;
721  int Fb;
722 
723  TString namefF(name);
724  namefF.ReplaceAll(".root","_Fisher.root");
725  TFile* fF =TFile::Open(namefF.Data());
726  TTree* FTree= (TTree*)fF->Get("Fisher_Discriminant");
727  FTree->SetBranchAddress("Fs_start",&Fs);
728  FTree->SetBranchAddress("Fb_start",&Fb);
729  FTree->SetBranchAddress("nFs",&nFs);
730  FTree->SetBranchAddress("nFb",&nFb);
731  FTree->SetBranchAddress("Fish_coef_cc",&w_c);
732  FTree->SetBranchAddress("Fish_coef_ANNout",&w_o);
733  FTree->GetEntry(0);
734  double w_i=0.;
735  int ind=0;
736  FTree->Branch("w_i",&w_i,"w_i/D");
737  FTree->Branch("index",&ind,"index/I");
738 
739  //fF->Close();
740  name.ReplaceAll("Fisher_file/","");
741 
742  double av;
743  double cc;
744  double rho;
745  int nSi;
746  NNTree2->SetBranchAddress("Average",&av);
747  NNTree2->SetBranchAddress("cc",&cc);
748  NNTree2->SetBranchAddress("rho",&rho);
749  NNTree2->SetBranchAddress("#TestSig",&nSi);
750  NNTree2->GetEntry(0);
751 
752  int const nSig=nSi;
753  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
754  cout<<" BG: "<<NNTree2->GetEntries()-nSi<<endl;
755  int nBgi=NNTree2->GetEntries()-nSi;
756  double minw_i=10000000;
757  double maxw_i=0.;
758 
759  for(int u=0;u<nSi;u++ ){
760  if (u>=Fs&&u<(Fs+nFs))continue;
761  ind=u;
762  NNTree2->GetEntry(u);
763  cout<<" u "<<u<<endl;
764  w_i=w_c*cc+w_o*av;
765  FTree->Fill();
766  if(w_i<minw_i)minw_i=w_i;
767  if(w_i>maxw_i)maxw_i=w_i;
768  }
769  int const nSigF=FTree->GetEntries();
770  cout<<" nSigF "<<nSigF<<endl;
771  for(int u=nSi;u<nSi+nBgi;u++ ){
772  if (u>=Fb&&u<(Fb+nFb))continue;
773  ind=u;
774  NNTree2->GetEntry(u);
775  cout<<" u "<<u<<endl;
776  w_i=w_c*cc+w_o*av;
777  FTree->Fill();
778  if(w_i<minw_i)minw_i=w_i;
779  if(w_i>maxw_i)maxw_i=w_i;
780  }
781 
782  cout<<"F Entries "<<FTree->GetEntries()<<endl;
783 
784  double deltaw=0.;
785  deltaw=(maxw_i-minw_i)/50.;
786  TH1D *hbg=new TH1D("Hbg", "projection[0]", 50, minw_i-deltaw, maxw_i+ deltaw);
787  TH1D *hsig=new TH1D("Hsig", "projection[0]", 50, minw_i-deltaw, maxw_i+deltaw);
788  hsig->SetDirectory(0);
789  hbg->SetDirectory(0);
790 
791 
792 
793  for( int u=0; u<FTree->GetEntries()-1;u++){ // //-1 perchèxk hai tot c'è in npiù l'uno da Fisher_ex: i nuovi branches ripartono da zero
794  FTree->GetEntry(u);
795  cout<<" u "<<u<<" ind "<<ind<<" w_i "<<w_i<<" nSigF-1 "<<nSigF-1<<endl;
796  if(u<nSigF-1) hsig->Fill(w_i);
797  else hbg->Fill(w_i);
798  }
799 
800 //exit(0);
801 //Efficiency_vs_rho","Efficiency_vs_rho",0,0,1200,700
802  TCanvas *c_histo= new TCanvas("Histogram_Fisher", "Histogram_Fisher", 0,0,900,600);
803  hbg->SetLineColor(kBlue);
804  hbg->SetFillStyle(3008); hbg->SetFillColor(kBlue);
805  hsig->SetLineColor(kRed);
806  hsig->SetFillStyle(3003); hsig->SetFillColor(kRed);
807  hsig->SetStats(0);
808  hbg->SetStats(0);
809 // isig->GetYaxis()->SetRangeUser(0,maxTrain/2);
810 // ibg->GetYaxis()->SetRangeUser(0,maxTrain/2);
811  hsig->Draw();
812  hbg->Draw("same");
813  TLegend *hlegend = new TLegend(.75, .80, .95, .95);
814  hlegend->AddEntry(hbg, "Background");
815  hlegend->AddEntry(hsig, "Signal");
816  hlegend->Draw();
817  char namech[2048];
818  char namechroot[2048];
819  TString namch(name);
820  namch.ReplaceAll(".root",".png");
821  sprintf(namech,"Fisher_png/Histog_%s",namch.Data());
822  sprintf(namechroot,"Fisher_png/Histog_%s",name.Data());
823  c_histo->SaveAs(namech);
824  c_histo->SaveAs(namechroot);
825 
826 /* int const ncurve=nANN*nCC;
827  cout<<"dentro funzione dopodef"<<endl;
828 //LOG(NBg)vs RHO------------------------------------------
829  double* rhoSig[ncurve];
830  for (int i=0;i<ncurve;i++) rhoSig[i]=new double[nSig];
831  cout<<"dopo def rhoSig"<<endl;
832  //double rhoSig[20][10];
833  int NSig[ncurve];
834 // int nBg=0;
835  int const nBg=NNTree2->GetEntries()-nSig;
836  for (int i=0;i<ncurve;i++) {
837  NSig[i]=0;
838  for (int j=0;j<nSig;j++) rhoSig[i][j]=0.;
839  }
840  cout<<"dopo def rhoSig"<<endl;
841  double* rhoBg[ncurve];
842  for (int i=0;i<ncurve;i++) rhoBg[i]=new double[nBg];
843  //double rhoBg[20][10];
844  int NBg[ncurve];
845  for (int i=0;i<ncurve;i++) {
846  NBg[i]=0;
847  for (int j=0;j<nBg;j++) rhoBg[i][j]=0.;
848  }
849  cout<<"dopo def rhoBg"<<endl;
850  double ccTh[nCC];
851  for (int i=0;i<nCC;i++) ccTh[i]=0.;
852  double NNTh[nANN];
853  for (int i=0;i<nANN;i++) NNTh[i]=0.;
854  double deltacc=0.;
855  double deltaANN=0.;
856  deltacc=0.2/nCC;
857  //deltaANN=0.6/(nANN-1);
858  deltaANN=(maxw_i-minw_i)/nANN;
859 
860  cout<<NNTree2->GetEntries()<<endl;
861  for(int n=0;n<NNTree2->GetEntries();n++){
862  cout<<n<<endl;
863  NNTree2->GetEntry(n);
864  w_i=w_c*cc+w_o*av;
865  cout<<"rho "<<rho<<" cc "<<cc<<" w_i "<<w_i<<endl;
866  for(int i=0; i<nCC;i++){
867  ccTh[i]=0.5+i*deltacc;
868  if(cc<ccTh[i]) continue;
869  for(int j=0; j<nANN;j++){
870  if(j==0) NNTh[j]=-1000.;
871  //else NNTh[j]=0.5+(j-1)*deltaANN;
872  //else NNTh[j]=0.1+(j-1)*deltaANN;
873  //else NNTh[j]=0.+(j)*deltaANN;
874  else NNTh[j]=minw_i+(j)*deltaANN;
875  //else NNTh[j]=0.+(j-1)*deltaANN/1000.;
876  //else NNTh[j]=1.;
877  if(w_i<NNTh[j]) continue;
878  int ni=0;
879  if(n>nSig) {
880 
881  NBg[i*nANN+j]= NBg[i*nANN+j]+1;
882  while(rhoBg[i*nANN+j][ni]!=0)ni=ni+1;
883  rhoBg[i*nANN+j][ni]=rho;
884  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
885 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
886  }
887  else {
888  NSig[i*nANN+j]= NSig[i*nANN+j]+1;
889  while(rhoSig[i*nANN+j][ni]!=0)ni=ni+1;
890  rhoSig[i*nANN+j][ni]=rho;
891  // cout<<"rho: "<<rho<<" colonna "<<i*nANN+j<<" riga "<<ni<<endl;
892 // cout<<" soglia_cc "<<ccTh[i]<<" soglia_ANN "<<NNTh[j]<<endl;
893  }
894  // }
895  }
896  }
897  }
898  cout<<"dopo riempimento variabili"<<endl;
899  int* indexS[ncurve];
900  for (int i=0;i<ncurve;i++) indexS[i]=new int[nSig];
901 
902  TGraph * gS[ncurve];
903  for (int y=0;y<ncurve;y++) {
904  int igS=0;
905  int igS_p=0;
906 // int indexS[nSig];
907 // double rhoS[nSig];
908 // for(int i=0;i<nSig;i++) {
909 // rhoS[i]=0.;
910 // indexS[i]=0;
911 // }
912  // ig=1;
913  // for(int i=0;i<nSig;i++) rhoS[i]=rhoSig[y][i];
914  gS[y]=new TGraph();
915 // gS[y]->SetMarkerStyle(7);
916  TMath::Sort(nSig,rhoSig[y],indexS[y],false);
917 // cout<<"dopo Sort "<<y<<endl;
918  for (int k=0;k<nSig;k++) {
919  int ii=indexS[y][k];
920  int yy=0;
921  if (k>0){
922 // cout<<"k "<<k<<endl;
923  int ij=indexS[y][k-1];
924  //if(rhoSig[y][ii]!=0&&rhoSig[y][ii]!=rhoSig[y][ij]) {
925  if(rhoSig[y][ii]!=0) {
926  yy=NSig[y]-igS;
927  //gS[y]->SetPoint(igS,rhoSig[y][ii],yy);
928  if(rhoSig[y][ii]!=rhoSig[y][ij]) gS[y]->SetPoint(igS_p++,rhoSig[y][ii],yy);
929  cout<<"igS"<<igS<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
930  igS=igS+1;
931  }
932  }
933  else {
934  if(rhoSig[y][ii]!=0){
935  yy=NSig[y]-igS;
936  gS[y]->SetPoint(0,rhoSig[y][ii],yy);
937  igS=igS+1;
938 // cout<<" x "<<rhoSig[y][ii]<<" y: "<<yy<<endl;
939  }
940 
941  }
942  }
943  }
944  // cout<<"dopo inserimento puntiiS"<<endl;
945 // cout<<ncurve<<endl;
946  int* indexB[ncurve];
947  for (int i=0;i<ncurve;i++) indexB[i]=new int[nBg];
948 
949  TGraph * gB[ncurve];
950  for (int y=0;y<ncurve;y++) {
951  int igB=0;
952  int igB_p=0;
953  gB[y]=new TGraph();
954  TMath::Sort(nBg,rhoBg[y],indexB[y],false);
955 // cout<<"dopo Sort "<<y<<endl;
956  for (int k=0;k<nBg;k++) {
957  int ii=indexB[y][k];
958  int yy=0;
959  if (k>0){
960  int ij=indexB[y][k-1];
961  //if(rhoBg[y][ii]!=0&&rhoBg[y][ii]!=rhoBg[y][ij]) {
962  if(rhoBg[y][ii]!=0) {
963  yy=NBg[y]-igB;
964  //gB[y]->SetPoint(igB,rhoBg[y][ii],yy);
965  if(rhoBg[y][ii]!=rhoBg[y][ij]) gB[y]->SetPoint(igB_p++,rhoBg[y][ii],yy);
966  igB=igB+1;
967 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
968  }
969  }
970  else {
971  if(rhoBg[y][ii]!=0) {
972  yy=NBg[y]-igB;
973  gB[y]->SetPoint(0,rhoBg[y][ii],yy);
974  igB=igB+1;
975 // cout<<"igB"<<igB<<" x "<<rhoBg[y][ii]<<" y: "<<yy<<endl;
976  }
977  }
978  }
979  }
980  cout<<"dopo inserimento puntiB"<<endl;
981  fTEST->Close();
982  fF->Close();
983  //gB[1]->SetMarkerStyle(7);
984  //gB[1]->SetPoint(1,1,2);
985  TCanvas* cS=new TCanvas("Efficiency_vs_rho","Efficiency_vs_rho",0,0,1200,700);
986  cS->Divide(2,2);
987  cS->cd(1)->SetLogy();
988  TMultiGraph* mg1=new TMultiGraph();
989 // gS[0]->SetMarkerStyle(7);
990  gS[0]->SetMarkerColor(2);
991  gS[0]->SetLineColor(2);
992  mg1->SetTitle("cc=0.5;rho;#Events");
993  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
994 // gS[0]->Draw("apl");
995  for (int h=1;h<nANN;h++){
996  gS[h]->SetMarkerColor(3);
997  gS[h]->SetLineColor(3);
998 // gS[h]->SetMarkerStyle(h+1);
999  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
1000  // gS[h]->Draw("apl,same");
1001  }
1002  mg1->Draw("apl");
1003  cS->cd(2)->SetLogy();
1004  TMultiGraph* mg2=new TMultiGraph();
1005 // gS[nANN]->SetMarkerStyle(7);
1006  gS[nANN]->SetMarkerColor(2);
1007  gS[nANN]->SetLineColor(2);
1008  mg2->SetTitle("cc=0.55;rho;#Events");
1009  if(gS[nANN]->GetN()!=0) mg2->Add(gS[nANN]);
1010  for (int h=1;h<nANN;h++){
1011  gS[nANN+h]->SetMarkerColor(3);
1012  gS[nANN+h]->SetLineColor(3);
1013 // gS[nANN+h]->SetMarkerStyle(h+1);
1014  if(gS[nANN+h]->GetN()!=0) mg2->Add(gS[nANN+h]);
1015  // gS[nANN+h]->Draw("apl,same");
1016  }
1017  mg2->Draw("apl");
1018  cS->cd(3)->SetLogy();
1019  TMultiGraph* mg3=new TMultiGraph();
1020 // gS[nANN*2]->SetMarkerStyle(7);
1021  gS[nANN*2]->SetMarkerColor(2);
1022  gS[nANN*2]->SetLineColor(2);
1023  mg3->SetTitle("cc=0.6;rho;#Events");
1024  if(gS[nANN*2]->GetN()!=0) mg3->Add(gS[nANN*2]);
1025  for (int h=1;h<nANN;h++){
1026  gS[2*nANN+h]->SetMarkerColor(3);
1027  gS[2*nANN+h]->SetLineColor(3);
1028 // gS[2*nANN+h]->SetMarkerStyle(h+1);
1029  if(gS[2*nANN+h]->GetN()!=0) mg3->Add(gS[2*nANN+h]);
1030  // gS[2*nANN+h]->Draw("apl,same");
1031  }
1032  mg3->Draw("apl");
1033  cS->cd(4)->SetLogy();
1034  TMultiGraph* mg4=new TMultiGraph();
1035 // gS[nANN*3]->SetMarkerStyle(7);
1036  gS[nANN*3]->SetMarkerColor(2);
1037  gS[nANN*3]->SetLineColor(2);
1038  mg4->SetTitle("cc=0.65;rho;#Events");
1039  if(gS[nANN*3]->GetN()!=0) mg4->Add(gS[nANN*3]);
1040  // gS[nANN*3]->Draw("apl");
1041  for (int h=1;h<nANN;h++){
1042  gS[3*nANN+h]->SetMarkerColor(3);
1043  gS[3*nANN+h]->SetLineColor(3);
1044 // gS[3*nANN+h]->SetMarkerStyle(h+1);
1045  if(gS[3*nANN+h]->GetN()!=0) mg4->Add(gS[3*nANN+h]);
1046  // gS[3*nANN+h]->Draw("apl,same");
1047  }
1048  mg4->Draw("apl");
1049  cout<<"nuovo canv"<<endl;
1050  TCanvas* cB=new TCanvas("Number_vs_rho","Number_vs_rho",0,0,1200,700);
1051  cB->Divide(2,2);
1052  cB->cd(1)->SetLogy();
1053  TMultiGraph* mg1B=new TMultiGraph();
1054 // gB[0]->SetMarkerStyle(7);
1055  gB[0]->SetMarkerColor(2);
1056  gB[0]->SetLineColor(2);
1057  mg1B->SetTitle("cc=0.5;rho;#Events");
1058  if(gB[0]->GetN()!=0) mg1B->Add(gB[0]);
1059  for (int h=1;h<nANN;h++){
1060  gB[h]->SetMarkerColor(3);
1061  gB[h]->SetLineColor(3);
1062  // gB[h]>SetMarkerStyle(h+1);
1063  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
1064  //gB[h]->Draw("apl,same");
1065  }
1066  mg1B->Draw("apl");
1067  cB->cd(2)->SetLogy();
1068  TMultiGraph* mg2B=new TMultiGraph();
1069  //gB[nANN]->SetMarkerStyle(7);
1070  gB[nANN]->SetMarkerColor(2);
1071  gB[nANN]->SetLineColor(2);
1072  mg2B->SetTitle("cc=0.55;rho;#Events");
1073  if(gB[nANN]->GetN()!=0) mg2B->Add(gB[nANN]);
1074  for (int h=1;h<nANN;h++){
1075  gB[nANN+h]->SetMarkerColor(3);
1076  gB[nANN+h]->SetLineColor(3);
1077  //gB[nANN+h]>SetMarkerStyle(h+1);
1078  if(gB[nANN+h]->GetN()!=0) mg2B->Add(gB[nANN+h]);
1079  //gB[h]->Draw("apl,same");
1080  }
1081  mg2B->Draw("apl");
1082  cB->cd(3)->SetLogy();
1083  TMultiGraph* mg3B=new TMultiGraph();
1084 // gB[2*nANN]->SetMarkerStyle(7);
1085  gB[2*nANN]->SetMarkerColor(2);
1086  gB[2*nANN]->SetLineColor(2);
1087  mg3B->SetTitle("cc=0.6;rho;#Events");
1088  if(gB[2*nANN]->GetN()!=0) mg3B->Add(gB[2*nANN]);
1089  for (int h=1;h<nANN;h++){
1090  gB[2*nANN+h]->SetMarkerColor(3);
1091  gB[2*nANN+h]->SetLineColor(3);
1092 // gB[2*nANN+h]>SetMarkerStyle(h+1);
1093  if(gB[2*nANN+h]->GetN()!=0) mg3B->Add(gB[2*nANN+h]);
1094  //gB[h]->Draw("apl,same");
1095  }
1096  mg3B->Draw("apl");
1097  cB->cd(4)->SetLogy();
1098  TMultiGraph* mg4B=new TMultiGraph();
1099 // gB[3*nANN]->SetMarkerStyle(7);
1100  gB[3*nANN]->SetMarkerColor(2);
1101  gB[3*nANN]->SetLineColor(2);
1102  mg4B->SetTitle("cc=0.65;rho;#Events");
1103  if(gB[3*nANN]->GetN()!=0) mg4B->Add(gB[3*nANN]);
1104  for (int h=1;h<nANN;h++){
1105  gB[3*nANN+h]->SetMarkerColor(3);
1106  gB[3*nANN+h]->SetLineColor(3);
1107 // gB[3*nANN+h]>SetMarkerStyle(h+1);
1108  if(gB[3*nANN+h]->GetN()!=0) mg4B->Add(gB[3*nANN+h]);
1109  //gB[h]->Draw("apl,same");
1110  }
1111  mg4B->Draw("apl");
1112 
1113 // cout<<"dopo Draw()"<<endl;
1114  TString CnameS(name);
1115  TString CnameB(name);
1116  TString CnameSroot(name);
1117  TString CnameBroot(name);
1118  char CnameS2[1024];
1119  char CnameB2[1024];
1120  char CnameS2root[1024];
1121  char CnameB2root[1024];
1122  CnameS.ReplaceAll(".root",".png");
1123  CnameB.ReplaceAll(".root",".png");
1124  sprintf(CnameS2,"logN_rho_Fish/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameS.Data());
1125  sprintf(CnameB2,"logN_rho_Fish/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameB.Data());
1126  sprintf(CnameS2root,"logN_rho_Fish/logN_rho_S_dANN%1.2f_%s",deltaANN,CnameSroot.Data());
1127  sprintf(CnameB2root,"logN_rho_Fish/logN_rho_B_dANN%1.2f_%s",deltaANN,CnameBroot.Data());
1128  cS->SaveAs(CnameS2);
1129  cB->SaveAs(CnameB2);
1130  cS->Print(CnameS2root);
1131  cB->Print(CnameB2root);
1132 // cout<<"fine"<<endl;
1133  */
1134 }
1135 
1136 /*
1137 void Annth(TString ifile){
1138  TString name(ifile);
1139  name.ReplaceAll("average_file/","Fisher_file/");
1140  TFile* fTEST =TFile::Open(ifile.Data());
1141  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1142 
1143  double w_o;
1144  double w_c;
1145  int nFs;
1146  int nFb;
1147  int Fs;
1148  int Fb;
1149 
1150  TString namefF(name);
1151  namefF.ReplaceAll(".root","_Fisher.root");
1152  TFile* fF =TFile::Open(namefF.Data());
1153  TTree* FTree= (TTree*)fF->Get("Fisher_Discriminant");
1154  FTree->SetBranchAddress("Fs_start",&Fs);
1155  FTree->SetBranchAddress("Fb_start",&Fb);
1156  FTree->SetBranchAddress("nFs",&nFs);
1157  FTree->SetBranchAddress("nFb",&nFb);
1158  FTree->SetBranchAddress("Fish_coef_cc",&w_c);
1159  FTree->SetBranchAddress("Fish_coef_ANNout",&w_o);
1160  FTree->GetEntry(0);
1161  double w_i;
1162  int ind;
1163  FTree->SetBranchAddress("w_i",&w_i,"w_i/D");
1164  FTree->SetBranchAddress("index",&ind,"index/I");
1165 
1166  //TString name(ifile);
1167  //name.ReplaceAll("average_file/","");
1168  //TFile* fTEST =TFile::Open(ifile.Data());
1169  //TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1170  double av;
1171  double cc;
1172  double rho;
1173  int nSi;
1174  NNTree2->SetBranchAddress("Average",&av);
1175  NNTree2->SetBranchAddress("cc",&cc);
1176  NNTree2->SetBranchAddress("rho",&rho);
1177  NNTree2->SetBranchAddress("#TestSig",&nSi);
1178  NNTree2->GetEntry(0);
1179  int const nSig=nSi;
1180 // cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1181  int const ncurve2=nRHO*nCC;
1182 
1183 
1184  double* ANNSig[ncurve2];
1185  for (int i=0;i<ncurve2;i++) ANNSig[i]=new double[nSig];
1186  //double rhoSig[20][10];
1187  int NSig[ncurve2];
1188 // int nBg=0;
1189  int const nBg=NNTree2->GetEntries()-nSig;
1190  for (int i=0;i<ncurve2;i++) {
1191  NSig[i]=0;
1192  for (int j=0;j<nSig;j++) ANNSig[i][j]=0.;
1193  }
1194  double* ANNBg[ncurve2];
1195  for (int i=0;i<ncurve2;i++) ANNBg[i]=new double[nBg];
1196 
1197  //double rhoBg[20][10];
1198  int NBg[ncurve2];
1199  for (int i=0;i<ncurve2;i++) {
1200  NBg[i]=0;
1201  for (int j=0;j<nBg;j++) ANNBg[i][j]=0.;
1202  }
1203  double ccTh[nCC];
1204  for (int i=0;i<nCC;i++) ccTh[i]=0.;
1205  double rhoTh[nRHO];
1206  for (int i=0;i<nRHO;i++) rhoTh[i]=0.;
1207  double deltacc=0.;
1208  double deltarho=0.;
1209  deltacc=0.2/nCC;
1210  deltarho=1./(nRHO);
1211 
1212 
1213  for(int n=0;n<FTree->GetEntries();n++){
1214  FTree->GetEntry(n)
1215  NNTree2->GetEntry(ind);
1216  cout<<"rho "<<rho<<" cc "<<cc<<" w_i "<<w_i<<endl;
1217  for(int i=0; i<nCC;i++){
1218  ccTh[i]=0.5+i*deltacc;
1219  if(cc<ccTh[i]) continue;
1220  for(int j=0; j<nRHO;j++){
1221  rhoTh[j]=5+j*deltarho;
1222  if(rho<rhoTh[j]) continue;
1223  int ni=0;
1224  if(n>nSig) {
1225  NBg[i*nRHO+j]= NBg[i*nRHO+j]+1;
1226  while(ANNBg[i*nRHO+j][ni]!=0)ni=ni+1;
1227  ANNBg[i*nRHO+j][ni]=w_i;
1228  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1229  }
1230  else {
1231  NSig[i*nRHO+j]= NSig[i*nRHO+j]+1;
1232  while(ANNSig[i*nRHO+j][ni]!=0)ni=ni+1;
1233  ANNSig[i*nRHO+j][ni]=w_i;
1234  // cout<<" soglia_cc "<<ccTh[i]<<" soglia_RHO "<<rhoTh[j]<<endl;
1235  }
1236  }
1237  }
1238  }
1239 
1240  int* indexS[ncurve2];
1241  for (int i=0;i<ncurve2;i++) indexS[i]=new int[nSig];
1242 
1243  TGraph * gS[ncurve2];
1244  for (int y=0;y<ncurve2;y++) {
1245  int igS=0;
1246  int igS_p=0;
1247  gS[y]=new TGraph();
1248  TMath::Sort(nSig,ANNSig[y],indexS[y],false);
1249  for (int k=0;k<nSig;k++) {
1250  int ii=indexS[y][k];
1251  int yy=0;
1252  if (k>0){
1253  //cout<<"k "<<k<<endl;
1254  int ij=indexS[y][k-1];
1255  //if(ANNSig[y][ii]!=0&&ANNSig[y][ii]!=ANNSig[y][ij]) {
1256  if(ANNSig[y][ii]!=0) {
1257  yy=NSig[y]-igS;
1258  //gS[y]->SetPoint(igS,ANNSig[y][ii],yy);
1259  if(ANNSig[y][ii]!=ANNSig[y][ij]) gS[y]->SetPoint(igS_p++,ANNSig[y][ii],yy);
1260  // cout<<"igS"<<igS<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1261  igS=igS+1;
1262 
1263  }
1264  }
1265  else {
1266  if(ANNSig[y][ii]!=0){
1267  yy=NSig[y]-igS;
1268  gS[y]->SetPoint(0,ANNSig[y][ii],yy);
1269  igS=igS+1;
1270  // cout<<" x "<<ANNSig[y][ii]<<" y: "<<yy<<endl;
1271  }
1272 
1273  }
1274  }
1275  }
1276 
1277 
1278 
1279  int* indexB[ncurve2];
1280  for (int i=0;i<ncurve2;i++) indexB[i]=new int[nBg];
1281 
1282  TGraph * gB[ncurve2];
1283  for (int y=0;y<ncurve2;y++) {
1284  int igB=0;
1285  int igB_p=0;
1286  gB[y]=new TGraph();
1287  TMath::Sort(nBg,ANNBg[y],indexB[y],false);
1288  // cout<<"dopo Sort "<<y<<endl;
1289  for (int k=0;k<nBg;k++) {
1290  int ii=indexB[y][k];
1291  int yy=0;
1292  if (k>0){
1293  int ij=indexB[y][k-1];
1294  //if(ANNBg[y][ii]!=0&&ANNBg[y][ii]!=ANNBg[y][ij]) {
1295  if(ANNBg[y][ii]!=0) {
1296  yy=NBg[y]-igB;
1297  //gB[y]->SetPoint(igB,ANNBg[y][ii],yy);
1298  if(ANNBg[y][ii]!=ANNBg[y][ij]) gB[y]->SetPoint(igB_p++,ANNBg[y][ii],yy);
1299  igB=igB+1;
1300  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1301  }
1302  }
1303  else {
1304  if(ANNBg[y][ii]!=0) {
1305  yy=NBg[y]-igB;
1306  gB[y]->SetPoint(0,ANNBg[y][ii],yy);
1307  igB=igB+1;
1308  // cout<<"igB"<<igB<<" x "<<ANNBg[y][ii]<<" y: "<<yy<<endl;
1309  }
1310  }
1311  }
1312  }
1313 
1314 
1315  TCanvas* cS=new TCanvas("Efficiency_vs_ANN","Efficiency_vs_ANN",0,0,1200,700);
1316  cS->Divide(2,2);
1317  //cS->cd(1)->SetLogy();
1318  cS->cd(1);
1319  TMultiGraph* mg1=new TMultiGraph();
1320  mg1->SetTitle("cc=0.5;ANN;#Events");
1321  for (int h=0;h<nRHO;h++){
1322  gS[h]->SetLineColor(4);
1323  if(gS[h]->GetN()!=0) mg1->Add(gS[h]);
1324  }
1325  mg1->Draw("al");
1326  cS->cd(2);
1327  // cS->cd(2)->SetLogy();
1328  TMultiGraph* mg2=new TMultiGraph();
1329  mg2->SetTitle("cc=0.55;ANN;#Events");
1330  for (int h=0;h<nRHO;h++){
1331  gS[nRHO+h]->SetLineColor(4);
1332  if(gS[nRHO+h]->GetN()!=0) mg2->Add(gS[nRHO+h]);
1333  }
1334  mg2->Draw("al");
1335 
1336  //cS->cd(3)->SetLogy();
1337  cS->cd(3);
1338  TMultiGraph* mg3=new TMultiGraph();
1339  mg3->SetTitle("cc=0.6;ANN;#Events");
1340  for (int h=0;h<nRHO;h++){
1341  gS[2*nRHO+h]->SetLineColor(4);
1342  if(gS[2*nRHO+h]->GetN()!=0) mg3->Add(gS[2*nRHO+h]);
1343  }
1344  mg3->Draw("al");
1345  //cS->cd(4)->SetLogy();
1346  cS->cd(4);
1347  TMultiGraph* mg4=new TMultiGraph();
1348  mg4->SetTitle("cc=0.65;ANN;#Events");
1349  for (int h=0;h<nRHO;h++){
1350  gS[3*nRHO+h]->SetLineColor(4);
1351  if(gS[3*nRHO+h]->GetN()!=0) mg4->Add(gS[3*nRHO+h]);
1352  }
1353  mg4->Draw("al");
1354 
1355  TCanvas* cB=new TCanvas("Number_vs_ANN","Number_vs_ANN",0,0,1200,700);
1356  cB->Divide(2,2);
1357  cB->cd(1)->SetLogy();
1358  TMultiGraph* mg1B=new TMultiGraph();
1359  mg1B->SetTitle("cc=0.5;ANN;#Events");
1360  for (int h=0;h<nRHO;h++){
1361  gB[h]->SetLineColor(4);
1362  if(gB[h]->GetN()!=0) mg1B->Add(gB[h]);
1363  }
1364  mg1B->Draw("al");
1365  cB->cd(2)->SetLogy();
1366  TMultiGraph* mg2B=new TMultiGraph();
1367  mg2B->SetTitle("cc=0.55;ANN;#Events");
1368  for (int h=0;h<nRHO;h++){
1369  gB[nRHO+h]->SetLineColor(4);
1370  if(gB[nRHO+h]->GetN()!=0) mg2B->Add(gB[nRHO+h]);
1371  }
1372  mg2B->Draw("al");
1373  cB->cd(3)->SetLogy();
1374  TMultiGraph* mg3B=new TMultiGraph();
1375  mg3B->SetTitle("cc=0.6;ANN;#Events");
1376  for (int h=0;h<nRHO;h++){
1377  gB[2*nRHO+h]->SetLineColor(4);
1378  if(gB[2*nRHO+h]->GetN()!=0) mg3B->Add(gB[2*nRHO+h]);
1379  }
1380  mg3B->Draw("al");
1381  cB->cd(4)->SetLogy();
1382  TMultiGraph* mg4B=new TMultiGraph();
1383  mg4B->SetTitle("cc=0.6;ANN;#Events");
1384  for (int h=0;h<nRHO;h++){
1385  gB[3*nRHO+h]->SetLineColor(4);
1386  if(gB[3*nRHO+h]->GetN()!=0) mg4B->Add(gB[3*nRHO+h]);
1387  }
1388  mg4B->Draw("al");
1389 
1390 
1391  //cout<<"dopo Draw()"<<endl;
1392  TString CnameS(name);
1393  TString CnameB(name);
1394  TString CnameSroot(name);
1395  TString CnameBroot(name);
1396  char CnameS2[1024];
1397  char CnameB2[1024];
1398  char CnameS2root[1024];
1399  char CnameB2root[1024];
1400  CnameS.ReplaceAll(".root",".png");
1401  CnameB.ReplaceAll(".root",".png");
1402  sprintf(CnameS2,"ANNthres_Fish/N_ANN_S_%s",CnameS.Data());
1403  sprintf(CnameB2,"ANNthres_Fish/N_ANN_B_%s",CnameB.Data());
1404  sprintf(CnameS2root,"ANNthres_Fish/N_ANN_S_%s",CnameSroot.Data());
1405  sprintf(CnameB2root,"ANNthres_Fish/N_ANN_B_%s",CnameBroot.Data());
1406  cS->SaveAs(CnameS2);
1407  cB->SaveAs(CnameB2);
1408  cS->Print(CnameS2root);
1409  cB->Print(CnameB2root);
1410  //cout<<"fine"<<endl;
1411 
1412 //CARTELLA Annth
1413 
1414 
1415 
1416 }
1417 
1418 */
1419 void PlotsAv_cc(TString ifile){
1420  TString name(ifile);
1421  name.ReplaceAll("outfile/","");
1422  name.ReplaceAll("average_file/","");
1423  TFile* fTEST =TFile::Open(ifile.Data());
1424  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1425  //double outbis;
1426  double av;
1427  double out;
1428  double cc;
1429  double rho;
1430  int nSi;
1431 // NNTree2->SetBranchAddress("ANNout",&out);
1432  NNTree2->SetBranchAddress("Average",&av);
1433  NNTree2->SetBranchAddress("cc",&cc);
1434  NNTree2->SetBranchAddress("rho",&rho);
1435  NNTree2->SetBranchAddress("#TestSig",&nSi);
1436  NNTree2->GetEntry(0);
1437  int const nSig=nSi;
1438  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1439  int const nBg=NNTree2->GetEntries()-nSig;
1440  cout<<"nBg: "<<nBg<<" Entries: "<<NNTree2->GetEntries()<<endl;
1441 
1442  TGraph * gS[3];
1443  TGraph * gB[3];
1444  gB[0]=new TGraph();
1445  gS[0]=new TGraph();
1446 // cout<<"nSig: "<<nSig<<endl;
1447  for (int n = 0; n <NNTree2->GetEntries(); n++){
1448  NNTree2->GetEntry(n);
1449  if(n<nSig) {
1450  gS[0]->SetPoint(n,cc,av);
1451  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1452  }
1453  else {
1454  gB[0]->SetPoint(n-nSig,cc,av);
1455  cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1456  }
1457  }
1458 
1459 
1460  /*for (int a=0;a<nBg;a++){
1461  if(aRHOB[a]>=RHOth){
1462  for (int b=0;b<nth;b++){
1463  if(aCCB[a]>=cc1th[b]){
1464  for(int c=0;c<nth+1;c++){
1465  if(aANNB[a]>=ANN1th[c]) Z[b*nth+c]=Z[b*nth+c]+1;
1466  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1467  }
1468  }
1469  }
1470  }
1471  }*/
1472  gS[0]->SetMarkerColor(2);
1473  gB[0]->SetMarkerColor(4);
1474  gS[0]->SetMarkerStyle(6);
1475  gB[0]->SetMarkerStyle(7);
1476 
1477  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1478  c->Divide(1,2);
1479  c->cd(1);
1480  TMultiGraph* mg1=new TMultiGraph();
1481  mg1->SetTitle("Av_cc");
1482  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1483  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1484  mg1->Draw("ap");
1485  c->cd(2);
1486  TMultiGraph* mg2=new TMultiGraph();
1487  mg2->SetTitle("Av_cc");
1488  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1489  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1490  mg2->Draw("ap");
1491 
1492  cout<<" name "<<name<<endl;
1493  TString Cname(name);
1494  TString Cnameroot(name);
1495  char Cname2[1024];
1496  char Cname2root[1024];
1497  Cname.ReplaceAll(".root",".png");
1498  sprintf(Cname2,"average_png/out_Plots_%s",Cname.Data());
1499  sprintf(Cname2root,"average_png/out_Plots_%s",Cnameroot.Data());
1500  c->SaveAs(Cname2);
1501  c->Print(Cname2root);
1502 /*
1503  TCanvas* c2=new TCanvas("Plots_Bkg_on_Sig","Plots_Bkg_on_Sig",0,0,1200,700);
1504  c2->Divide(2,2);
1505  c2->cd(1);
1506  gS[0]->Draw("ap");
1507  gB[0]->Draw("p,same");
1508 
1509  c2->cd(2);
1510  gS[1]->Draw("ap");
1511  gB[1]->Draw("p,same");
1512  c2->cd(3);
1513  gS[2]->Draw("ap");
1514  gB[2]->Draw("p,same");
1515  c2->cd(4);
1516  TText* text2=new TText(0.37,0.0,"no cuts on ANN");
1517  hglitch->SetStats(0);
1518  hglitch->GetXaxis()->SetTitle("cc");
1519  hglitch->GetYaxis()->SetTitle("ANNoutput");
1520  hglitch->GetZaxis()->SetTitle("count");
1521  hglitch->Draw("colz");
1522  text2->Draw();
1523  gPad->SetLogz();
1524 
1525  TString Cname_2(name);
1526  TString Cname_2root(name);
1527  char Cname2_2[1024];
1528  char Cname2_2root[1024];
1529  Cname_2.ReplaceAll(".root",".png");
1530  sprintf(Cname2_2,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2.Data());
1531  sprintf(Cname2_2root,"CC_RHO_ANN_Plots/CC_RHO_ANN_Plots_BoS_%s",Cname_2root.Data());
1532  c2->SaveAs(Cname2_2);
1533  c2->Print(Cname2_2root);
1534 
1535 //CARTELLA CCi_RHO_ANN_Plots */
1536 }
1537 void PlotsAv_Mc(TString ifile){
1538  TString name(ifile);
1539  name.ReplaceAll("outfile/","");
1540  name.ReplaceAll("average_file/","");
1541  TFile* fTEST =TFile::Open(ifile.Data());
1542  TTree* NNTree2=(TTree*)fTEST->Get("Parameters");
1543  //double outbis;
1544  double av;
1545  double out;
1546  double cc;
1547  double Mc;
1548  double rho;
1549  int nSi;
1550 // NNTree2->SetBranchAddress("ANNout",&out);
1551  NNTree2->SetBranchAddress("Mchirp",&Mc);
1552  NNTree2->SetBranchAddress("Average",&av);
1553  NNTree2->SetBranchAddress("cc",&cc);
1554  NNTree2->SetBranchAddress("rho",&rho);
1555  NNTree2->SetBranchAddress("#TestSig",&nSi);
1556  NNTree2->GetEntry(0);
1557  int const nSig=nSi;
1558  cout<<"nSig: "<<nSig<<" nSi: "<<nSi<<endl;
1559  int const nBg=NNTree2->GetEntries()-nSig;
1560 
1561  TGraph * gS[3];
1562  TGraph * gB[3];
1563  gB[0]=new TGraph();
1564  gS[0]=new TGraph();
1565 // cout<<"nSig: "<<nSig<<endl;
1566  for (int n = 0; n <NNTree2->GetEntries(); n++){
1567  NNTree2->GetEntry(n);
1568  if(n<nSig) {
1569  gS[0]->SetPoint(n,Mc,av);
1570  cout<<"Sig_graph1: x="<<cc<<" y: "<<av<<endl;
1571  }
1572  else {
1573  gB[0]->SetPoint(n-nSig,Mc,av);
1574  cout<<"Bg_graph1: x="<<cc<<" y: "<<av<<endl;
1575  }
1576  }
1577 
1578 
1579  /*for (int a=0;a<nBg;a++){
1580  if(aRHOB[a]>=RHOth){
1581  if(c=10) cout<<" ANN "<<ANN1th[c]<<" Zcount: "<<Z[b*nth+c]<<" aNNB "<<aANNB[a]<<" a "<<a<<" b "<<b<<endl;
1582  }
1583  }
1584  }
1585  }
1586  }*/
1587  gS[0]->SetMarkerColor(2);
1588  gB[0]->SetMarkerColor(4);
1589  gS[0]->SetMarkerStyle(6);
1590  gB[0]->SetMarkerStyle(7);
1591 
1592  TCanvas* c=new TCanvas("Plots","Plots",0,0,1200,700);
1593  c->Divide(1,2);
1594  c->cd(1);
1595  TMultiGraph* mg1=new TMultiGraph();
1596  mg1->SetTitle("Av_Mc");
1597  if(gB[0]->GetN()!=0) mg1->Add(gB[0]);
1598  if(gS[0]->GetN()!=0) mg1->Add(gS[0]);
1599  mg1->Draw("ap");
1600  c->cd(2);
1601  TMultiGraph* mg2=new TMultiGraph();
1602  mg2->SetTitle("Av_Mc");
1603  if(gS[0]->GetN()!=0) mg2->Add(gS[0]);
1604  if(gB[0]->GetN()!=0) mg2->Add(gB[0]);
1605  mg2->Draw("ap");
1606 
1607  cout<<" name "<<name<<endl;
1608  TString Cname(name);
1609  TString Cnameroot(name);
1610  char Cname2[1024];
1611  char Cname2root[1024];
1612  Cname.ReplaceAll(".root",".png");
1613  sprintf(Cname2,"average_png/Mc_Plots_%s",Cname.Data());
1614  sprintf(Cname2root,"average_png/Mc_Plots_%s",Cnameroot.Data());
1615  c->SaveAs(Cname2);
1616  c->Print(Cname2root);
1617 }
rho
double rho
Definition: cwb_report_prod_2.C:149
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:93
g
static double g(double e)
Definition: GNGen.cc:98
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
count
int count
Definition: compare_bkg.C:373
graph_Fish
void graph_Fish(TString ifile)
Definition: Fisher_1.C:710
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.
ofile
char ofile[512]
Definition: cwb_condor_rescue.C:40
i
i drho i
Definition: cwb_epparameters.C:85
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
g1
TF1 * g1
Definition: compare_bkg.C:505
PlotsAv_cc
void PlotsAv_cc(TString ifile)
Definition: Fisher_1.C:1419
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
netevent::netcc
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:94
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
ind
int * ind
Definition: cwb_fix_live_slag_missed.C:89
Fisher_ex
void Fisher_ex(TString ifile, int nFs, int nFb, int Fs, int Fb)
Definition: Fisher_1.C:449
s
s s
Definition: cwb_net.C:137
Fisher_1
void Fisher_1(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 nFs, int nFb, int Fs=0, int Fb=0, int TS=0, int TB=0, int s=0, int b=0, int uf=1, int ni=0)
Definition: Fisher_1.C:44
PlotsAv_Mc
void PlotsAv_Mc(TString ifile)
Definition: Fisher_1.C:1537
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
wc
netcluster wc
Definition: DrawClusterWithSkyplot.C:39
nIFO
#define nIFO
Definition: Fisher_1.C:24
iMc
#define iMc
Definition: Fisher_1.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
y
wavearray< double > y
Definition: Test10.C:31
netevent::chirp
Float_t * chirp
range to source: [0/1]-rec/inj
Definition: netevent.hh:109
exit
exit(0)