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