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