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