Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cwbANN_purity_BoS.C
Go to the documentation of this file.
1 //
2 // Neural Network Analysis
3 // Author : Serena Vinciguerra
4 // Note : run with ACLiC -> root -l macro/cwb_mlp.C+
5 //
6 
7 #define XIFO 4
8 
9 #pragma GCC system_header
10 
11 #include "cwb.hh"
12 #include "config.hh"
13 #include "network.hh"
14 #include "wavearray.hh"
15 #include "TString.h"
16 #include "TObjArray.h"
17 #include "TObjString.h"
18 #include "TPaletteAxis.h"
19 #include "TMultiLayerPerceptron.h"
20 #include "TMLPAnalyzer.h"
21 #include "TRandom.h"
22 #include "TComplex.h"
23 #include "TMath.h"
24 #include "mdc.hh"
25 #include "watplot.hh"
26 #include <vector>
27 
28 
29 //#define SIG_FILE "/home/serena.vinciguerra/test0/S6D_R11_SIM_EOBNRv2_advL1H1V1_2G_run1_2/merge/NN_wave_S6D_R11_SIM_EOBNRv2_advL1H1V1_2G_run1_2.M1.root"
30 //#define BG_FILE "/home/serena.vinciguerra/test0/S6D_R11_BKG_L1H1V1_2G_MP_run1/merge/NN_wave_S6D_R11_BKG_L1H1V1_2G_MP_run1.M1.root"
31 #define NRHO 50
32 #define NCC 50
33 #define nIFO 3 // number of detectors
34 //#define differentW
35 #define RHO_CC
36 #define CREATE_NETWORK
37 #define trainTEST
38 TString select7Train_non_es(TString ORIGINAL_FILE,std::vector<int>* choice,TString IDname);
39 
40 //NOTA_BENE:OGNI VOLTA CHE VIENE RICHIAMATA LA MACRO LA NN PRECEDENTE VIENE SOVRASCRITTA
41 //NOTA_BENE: OGNI VOLTA CHE FACCIO UN GRAFICO ROC SI SOVRASCRIVE SE CAMBIO SOLO COSE DIVERSE DA NTRAINBG E NTRAINS E Y
42 //NOTA_BENE: DOVE SI USA QST MACRO DEVONO ESISTERE LA CARRTELLA png(se def Graph), tree, NN
43 
44 void cwbANN_purity_BoS(int nTrain,int nTrainBg, TString option,TString TREE_FILE,int q,int p,double rm=0., double cm=0., int lm=5, Int_t nepoch=700) {
45 
46  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
47  gSystem->Load("libMLP");
48  }
49 
50  //gSystem->Load("select7_C.so");
51 // gROOT->LoadMacro("select7Train_non_es.C");
52 // cout<<soglia<<endl;
53  double error;
54  //int y=0;
55 // int q=0;//sig
56 // int p=0;//bg
57  double maxTrain=0;
58  // cout<<"y:"<<y<<endl;
59  if(nTrain>nTrainBg) maxTrain=nTrain;
60  else maxTrain=nTrainBg;
61 
62  TFile* ifile=TFile::Open(TREE_FILE.Data());
63  TTree* NNTree=(TTree*)ifile->Get("nnTree");
64  //NNTree->SetDirectory(0);
65  //ifile->Close();
66 
67  int entries=NNTree->GetEntries();
68  cout<<"entries: "<<entries<<endl;
69 
70 //estraggo le info che servono anche al tree di def dell'mlp
71  int ndim;
72  int ninp;
73  int y;
74  int entriesTot;
75  int bg_entries;
76  int sig_entries;
77  NNTree->SetBranchAddress("#Entries_type",&entriesTot);
78  NNTree->SetBranchAddress("Matrix_dim",&ndim);
79  NNTree->SetBranchAddress("#inputs",&ninp);
80  NNTree->SetBranchAddress("amplitude_mode",&y);
81 
82  NNTree->GetEntry(0);
83  sig_entries=entriesTot;
84  NNTree->GetEntry(entries-1);
85  bg_entries=entriesTot;
86 
87  int const NDIM=ndim;
88  int const nINP=ninp;
89  cout<<"NDIM "<<NDIM<<endl;
90  cout<<"nINP"<<nINP<<endl;
91 
92  cout<<"sig e "<<sig_entries<<endl;
93  cout<<"bg e "<<bg_entries<<endl;
94 
95 //contolli
96  if(sig_entries==bg_entries) cout<<"Error: training with a single type of event"<<endl;
97 
98  if(bg_entries<nTrainBg){cout<<"Error: Bg: num train>num events"<<endl;
99  exit(0);}
100  if(sig_entries<nTrain) {cout<<"Error: Sig: num train>num events"<<endl;
101  exit(0);}
102 
103  if(!nTrainBg) nTrainBg=bg_entries;
104  if(!nTrain) nTrain=sig_entries;
105 //fine controlli
106 
107 //definizione start
108  // q=gRandom->Integer(sig_entries-nTrain);
109  // p=gRandom->Integer(bg_entries-nTrainBg)+sig_entries;
110  //q=1;
111  //p=1;
112 
113  TString TreeF(TREE_FILE);
114  TreeF.ReplaceAll("/home/serena.vinciguerra/","");
115  TreeF.ReplaceAll("/","_");
116  TreeF.ReplaceAll(".root","");
117  cout<<"p "<<p<<endl;
118  cout<<"q "<<q<<endl;
119 //fine defnizione start
120  TString structure(option);
121  structure.ReplaceAll(":","_");
122  char nomeNN[512];
123  sprintf(nomeNN,"NN/mlpNetwork_nTS%i_nTB%i_structure%s_lm%i_epochs%i_s%i_b%i_ROC_TREEfile%s.root",nTrain,nTrainBg,structure.Data(),lm,nepoch,q,p,TreeF.Data());
124 cout<<nomeNN<<endl;
125  char nameTrain[512];
126  sprintf(nameTrain,"PropertyGraphs/mlpNetwork_nTS%i_nTB%i_s%i_b%i_TREEfile_%s.root",nTrain,nTrainBg,q,p,TreeF.Data());
127 #ifdef RHO_CC
128 TString SIG_FILE;
129 TString BG_FILE;
130 char FILE_NAME[516];
131 NNTree->SetBranchAddress("Files_name",&FILE_NAME);
132 NNTree->GetEntry(0);
133 SIG_FILE=FILE_NAME;
134 NNTree->GetEntry(entries-1);
135 BG_FILE=FILE_NAME;
136 cout<<"fine ifdef RHO_CC"<<endl;
137 #endif
138 
139 TString FileSelectName(nameTrain);
140 TString nomefile;
141 FileSelectName.ReplaceAll("PropertyGraphs/mlpNetwork","selectTREE/TreeTrain");
142 TFile *f0=TFile::Open(FileSelectName.Data(),"read");
143 
144 if (!f0) {
145 
146 //RICOSTRUZIONE TREE per Train
147 std::vector<int>* choice=new vector<int>;
148 cout<<"vector definition"<<endl;
149 //choice->push_back(q);
150 for (int z=0;z<nTrain;z++) choice->push_back(q+z);
151 //(*choice)[nTrain]=p;
152 for (int z=0;z<nTrainBg;z++) choice->push_back(z+p);
153 // gSystem->Load("select7_C.so");
154 
155 nomefile=select7Train_non_es(TREE_FILE.Data(),choice,nameTrain);
156 cout<<"new file created"<<endl;
157 }
158 else {
159 nomefile=FileSelectName.Data();
160 cout<<"open file-Train pre-existing: "<<nomefile.Data()<<endl;
161 f0->Close();}
162 //exit(0);
163 //char nomefile[516];
164 //sprintf(nomefile,"selectTREE/TreeTrain_size%i_s%i_b%i.root",choice->size(),q,p);
165 TFile* fsel=TFile::Open(nomefile.Data());
166 TTree* tTrain1=(TTree*)fsel->Get("nnTree");
167 cout<<"def Tree per Train"<<endl;
168  // add leaf
169  Int_t type;
170  Float_t x[nINP];
171  char ilabel[nINP][16];
172  // char llabel[nINP][16];
173  for(int i=0;i<nINP;i++) {
174  sprintf(ilabel[i],"x%d",i+1);
175  // cout<<"ilabel "<<ilabel[i]<<endl;
176  // sprintf(llabel[i],"%s/F",ilabel[i]);
177  // tTrain->Branch(ilabel[i], &x[i],llabel[i]);
178  tTrain1->SetBranchAddress(ilabel[i], &x[i]);
179  }
180  tTrain1->SetBranchAddress("type",&type);
181  #ifdef differentW
182  Float_t w[nINP];
183  char wlabel[nINP][16];
184  // char llabel2[nINP][16];
185 
186  for(int i=0;i<nINP;i++) {
187  sprintf(wlabel[i],"w%d",i+1);
188  // sprintf(llabel2[i],"%s/F",wlabel[i]);
189  tTrain1->SetBranchAddress(wlabel[i], &w[i]);
190  }
191  #else
192  float w;
193  tTrain1->SetBranchAddress("w",&w);
194  #endif
195  ifile->Close();
196 
197 #ifdef RHO_CC
198  TChain sigTree("waveburst");//cerca il Tree "waveburst nei file
199  sigTree.Add(SIG_FILE.Data());//determina i file
200  netevent signal(&sigTree,nIFO);
201  int sig_entries2 = signal.GetEntries();
202  cout << "sig entries2 : " << sig_entries2 << endl;
203 
204  TChain bgTree("waveburst");
205  bgTree.Add(BG_FILE.Data());
206  netevent background(&bgTree,nIFO);
207  int bg_entries2 = background.GetEntries();
208  cout << "bg entries2 : " << bg_entries2 << endl;
209 
210  TH1F *ibg = new TH1F("bgh", "rho[0]", 50, 0, 50);
211  TH1F *isig = new TH1F("sigh", "rho[0]", 50, 0, 18);
212  TH1F *ibg2 = new TH1F("bgh2", "netcc[1]", 50, 0, 1.2);
213  TH1F *isig2 = new TH1F("sigh2", "netcc[1]", 50, 0, 1.1);
214 // TH1F *ibg3 = new TH1F("bgh3", "test secondo giro", 50, -1, 2);
215 // TH1F *isig3 = new TH1F("sigh3", "test secondo giro", 50, -1, 2);
216 
217  ibg->SetDirectory(0);
218  isig->SetDirectory(0);
219  ibg2->SetDirectory(0);
220  isig2->SetDirectory(0);
221 // ibg3->SetDirectory(0);
222 // isig3->SetDirectory(0);
223  Int_t i;
224  double minRHO=0.;
225  double minCC=0.;
226  double maxRHO=0.;
227  double maxCC=0.;
228 
229 TTree* t_NNout=new TTree("rho_cc_NNout","rho_cc_NNout");
230 double rho_NN=0.;
231 double cc_NN=0.;
232 double NNout=0.;
233 int NNtype;
234 t_NNout->Branch("rho",&rho_NN,"rho/D");
235 t_NNout->Branch("cc",&cc_NN,"cc/D");
236 t_NNout->Branch("type",&NNtype,"type/I");
237 t_NNout->Branch("NNout",&NNout,"NNout/D");
238 
239 for(int h=q;h<q+nTrain;h++){
240  int resto0=0;
241  if(q%2==0) resto0=(h)%2;
242  else resto0=(h+1)%2;
243  sigTree.GetEntry(h);
244  if (resto0==0&&h<q+2) minCC=(double)signal.netcc[1];
245  if (resto0==0&&h<q+2) minRHO=(double)signal.rho[0];
246  if (signal.rho[0]<minRHO&&resto0==0)minRHO=(double)signal.rho[0];
247  if (signal.netcc[1]<minCC&&resto0==0)minCC=(double)signal.netcc[1];
248  if (signal.rho[0]>maxRHO&&resto0==0)maxRHO=(double)signal.rho[0];
249  if (signal.netcc[1]>maxCC&&resto0==0)maxCC=(double)signal.netcc[1];
250  isig->Fill(signal.rho[0]);
251  isig2->Fill(signal.netcc[1]);
252  }
253 
254 for(int h=(p-sig_entries);h<(p-sig_entries+nTrainBg);h++){
255  int resto0=0;
256  if((p-sig_entries+nTrainBg)%2==0) resto0=(h)%2;
257  else resto0=(h+1)%2;
258  bgTree.GetEntry(h);
259  if (resto0==0&&h<q+2) minCC=(double)signal.netcc[1];
260  if (resto0==0&&h<q+2) minRHO=(double)signal.rho[0];
261  if (background.rho[0]<minRHO&&resto0==0)minRHO=(double)background.rho[0];
262  if (background.netcc[1]<minCC&&resto0==0)minCC=(double)background.netcc[1];
263  if (background.rho[0]>maxRHO&&resto0==0)maxRHO=(double)background.rho[0];
264  if (background.netcc[1]>maxCC&&resto0==0)maxCC=(double)background.netcc[1];
265  ibg->Fill(background.rho[0]);
266  ibg2->Fill(background.netcc[1]);
267  }
268 cout<<"max RHO"<<maxRHO<<" maxCC "<<maxCC<<" minRHO "<<minRHO<<" minCC "<<minCC<<endl;
269 #endif
270 if (cm!=0) maxCC=cm;
271 if (rm!=0) maxRHO=rm;
272 
273 
274 TString layout;
275 for (int i=0;i<nINP-1;i++) {
276 char add[512];
277 sprintf(add,"@%s,",ilabel[i]);
278 layout+=add;
279 
280 }
281 char add1[512];
282 sprintf(add1,"@%s:%s:type",ilabel[nINP-1],option.Data());
283 layout+=add1;
284 cout<<layout.Data()<<endl;
285 
286 // char layout[5096]="";
287 // for (int i=0;i<nINP-1;i++) sprintf(layout,"%s@%s,",layout,ilabel[i]);
288 // sprintf(layout,"%s@%s:%s:type",layout,ilabel[nINP-1],option.Data());
289  // cout << "MultiLayer Perceptron layout -> " << layout << endl;
290 
291 #ifdef CREATE_NETWORK
292  TCanvas *errors=new TCanvas("errors","errors",0,0,600,600);
293  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(layout,"w",tTrain1,"Entry$%2","(Entry$+1)%2");
294  if (lm==1)mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
295  if(lm==2)mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
296  if(lm==3)mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
297  if(lm==4)mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
298  if(lm==5)mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
299  if(lm==6)mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
300  if(lm>6||lm<1){
301  cout<<"Invalid Learning Method:1<=lm<=6"<<endl;
302  exit(0);
303  }
304 
305  mlp->Train(nepoch, "text,current,graph,update=10");
306  // char nomePNG[128];
307  TString PNGname(nomeNN);
308  TString PNGnameroot(nomeNN);
309  PNGname.ReplaceAll("NN/","ErrorGraphs/");
310  PNGnameroot.ReplaceAll("NN/","ErrorGraphs/");
311  PNGname.ReplaceAll(".root",".png");
312 // sprintf(nomePNG,"ErrorGraphs/mlpNetwork_nTS%i_nTB%i_lm%i_epochs%i_s%i_b%i.png",nTrain,nTrainBg,q,p);
313  errors->SaveAs(PNGname.Data());
314  errors->Print(PNGnameroot.Data());
315  // char nomeNN[128];
316  // sprintf(nomeNN,"NN/mlpNetwork_nTS%i_nTB%i__lm%i_epochs%i_s%i_b%i.root",nTrain,nTrainBg,q,p);
317  // TFile fnet(nomeNN,"RECREATE");
318  TFile fnet(nomeNN,"RECREATE");
319 
320  mlp->Write();
321 // fnet.Close();
322 
323 #else
324 cout<<"ci sono"<<endl;
325  TFile fnet(nomeNN);
326 cout<<"non ci sono più"<<endl;
327 // TFile fnet = TFile::Open(nomeNN);
328  TMultiLayerPerceptron *mlp = (TMultiLayerPerceptron*)fnet.Get("TMultiLayerPerceptron");
329  if(mlp==NULL) {cout << "Error getting mlp" << endl;exit(1);}
330  //fnet->Close();
331 #endif
332 
333  // Use TMLPAnalyzer to see what it looks for
334  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas","Network analysis");
335 #ifdef RHO_CC
336  mlpa_canvas->Divide(2,2);
337 
338  //DISEGNA L'ISTROGRMMA DI RHO
339  mlpa_canvas->cd(1);
340  ibg->SetLineColor(kBlue);
341  ibg->SetFillStyle(3008); ibg->SetFillColor(kBlue);
342  isig->SetLineColor(kRed);
343  isig->SetFillStyle(3003); isig->SetFillColor(kRed);
344  ibg->SetStats(0);
345  isig->SetStats(0);
346 // isig->GetYaxis()->SetRangeUser(0,maxTrain/2);
347 // ibg->GetYaxis()->SetRangeUser(0,maxTrain/2);
348  ibg->Draw();
349  isig->Draw("same");
350  TLegend *ilegend = new TLegend(.75, .80, .95, .95);
351  ilegend->AddEntry(ibg, "Background");
352  ilegend->AddEntry(isig, "Signal");
353  ilegend->Draw();
354  cout << "Integral ibg : " << ibg->GetEntries() << endl;
355  cout << "Integral isig : " << isig->GetEntries() << endl;
356 
357 
358  //DISEGNA L'ISTOGRMMA CC
359  mlpa_canvas->cd(2);
360  // ibg2->GetYaxis()->SetRangeUser(0,maxTrain/2);
361  // isig2->GetYaxis()->SetRangeUser(0,maxTrain/2);
362  ibg2->SetLineColor(kBlue);
363  ibg2->SetFillStyle(3008); ibg2->SetFillColor(kBlue);
364  isig2->SetLineColor(kRed);
365  isig2->SetFillStyle(3003); isig2->SetFillColor(kRed);
366  ibg2->SetStats(0);
367  isig2->SetStats(0);
368  ibg2->Draw();
369  isig2->Draw("same");
370  TLegend *ilegend2 = new TLegend(.75, .80, .95, .95);
371  ilegend2->AddEntry(ibg2, "Background");
372  ilegend2->AddEntry(isig2, "Signal");
373  ilegend2->Draw();
374  cout << "Integral ibg : " << ibg2->GetEntries() << endl;
375  cout << "Integral isig : " << isig2->GetEntries() << endl;
376 #endif
377 
378  // Build and train the NN ptsumf is used as a weight since we are primarly
379 
380  TMLPAnalyzer ana(mlp);
381  // Initialisation
382  ana.GatherInformations();
383  // output to the console
384  ana.CheckNetwork();
385 #ifdef CREATE_NETWORK
386 #ifdef RHO_CC
387  mlpa_canvas->cd(4);
388  // shows how each variable influences the network
389  ana.DrawDInputs();
390  // mlpa_canvas->cd(2);
391  // shows the network structure
392  // mlp->Draw();
393  //ana.DrawDInputs();
394  mlpa_canvas->cd(3);
395  // draws the resulting network
396  gPad->SetLogy();
397  ana.DrawNetwork(0,"type==1","type==0"); //disegna solo il test sample
398 #else
399 mlpa_canvas->Divide(1,2);
400 mlpa_canvas->cd(1);
401  ana.DrawDInputs();
402 mlpa_canvas->cd(2);
403  gPad->SetLogy();
404  ana.DrawNetwork(0,"type==1","type==0"); //disegna solo il test sample
405 #endif
406 
407 //DEFINISCO IL TREE DOVE SALVO TUTTE LE INFO SUL SET TEST
408 
409 TTree* t= new TTree("info","info");
410 t->Branch("amplitude_mode",&y,"amplitude_mode/I");
411 // int nDim=NDIM;
412 t->Branch("Matrix_dim",&ndim,"Matrix_dim/I");
413 // int nInp=nINP;
414 t->Branch("#inputs",&ninp,"#inputs/I");
415  char NNname[516];
416  sprintf(NNname,"%s",option.Data());
417 t->Branch("NNname",&NNname,"NNname/C");
418 t->Branch("#trainBg",&nTrainBg,"#trainBg/I");
419 t->Branch("#trainSig",&nTrain,"#trainSig/I");
420 //double threshold;
421 //t->Branch("threshold",&soglia,"threshold/D");
422 t->Branch("Rand_start_Bg",&p,"Rand_start_Bg/I");
423 t->Branch("Rand_start_Sig",&q,"Rand_start_Sig/I");
424 t->Branch("epochs",&nepoch,"epocs/I");
425 t->Branch("learning_method",&lm,"learning_method/I");
426 t->Fill();
427 
428 t->Write();
429 //t->SetDirectory(fnet);
430 
431 double errTot=mlp->GetError(TMultiLayerPerceptron::kTest);
432 
433 
434 mlpa_canvas->cd(0);
435 mlpa_canvas->Write();
436 TString CanvNN(nomeNN);
437 TString CanvNNroot(nomeNN);
438 CanvNN.ReplaceAll("NN/","mlpa_canv/");
439 CanvNNroot.ReplaceAll("NN/","mlpa_canv/");
440 CanvNN.ReplaceAll(".root",".png");
441 mlpa_canvas->SaveAs(CanvNN.Data());
442 mlpa_canvas->Print(CanvNNroot.Data());
443 #endif
444 fnet.Close();
445 //fnet.Close();
446 
447 //Testin the Train sample
448  int const npunti=40;
449  double params[nINP];
450  float TA[npunti+1];
451  float FA[npunti+1];
452  float thres0[npunti+1];
453  float TA1[npunti+1];
454  float FA1[npunti+1];
455 // float thres01[npunti+1];
456 
457  for (int i=0;i<=npunti;i++){
458  TA[i]=0;
459  TA1[i]=0;
460  FA[i]=0;
461  FA1[i]=0;
462  // thres01[i]=0;
463  thres0[i]=0;
464  }
465  float passo=0;
466  passo=(log(1.1)-log(0.5))/(npunti/2);
467  cout<<"passo "<<passo<<endl;
468  thres0[npunti]=0.5;
469  cout<<"soglia0 "<<thres0[0]<<endl;
470  for(int i=npunti;i>npunti/2;i--) {
471  thres0[i-1]=thres0[i]*exp(passo);}
472  for(int i=0;i<(npunti/2);i++) {
473  thres0[npunti/2+i+1]=1.6 -thres0[npunti/2+i+1];
474  thres0[npunti/2-i-1]=1.-thres0[npunti/2+i+1];
475 // cout<<"i: "<<i+11<<" "<<thres0[11+i]<<endl;
476 // cout<<"i: "<<11-i<<" "<<thres0[11-i]<<endl;
477 // thres0[i+1]=thres0[i]+passo;
478  }
479 thres0[npunti/2]=0.5;
480  for(int i=0;i<(npunti);i++) {
481  cout<<thres0[i]<<endl;
482  }
483 
484 /*for(int i=0;i<(npunti/2);i++) {
485  thres0[npunti/2+i+1]=thres0[npunti/2+i]*exp(passo);}
486  for(int i=0;i<(npunti/2);i++) {
487  thres0[npunti/2+i+1]=1.5 -thres0[npunti/2+i+1];
488  thres0[npunti/2-i-1]=1.-thres0[npunti/2+i+1];
489 // thres0[i+1]=thres0[i]+passo;
490  }
491 */
492 //TH2D* PUREZZA_S=new TH2D("Purezza","Purezza",NCC,0.,1.2,NRHO,0.,50);
493 TH2D* PUREZZA_S=new TH2D("Purezza","Purezza",NCC,minCC,maxCC,NRHO,minRHO,maxRHO);
494 TH2D* h2BgO_R=new TH2D("rho_cc_out","rho_cc_out",50,0.,50,50,-0.5,1.5);
495 TH2D* h2SigO_R=new TH2D("rho_cc_out","rho_cc_out",50,0.,50,50,-0.5,1.5);
496 TH2D* h2BgR_C=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,0.,50);
497 TH2D* h2SigR_C=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,0.,50);
498 TH2D* h2Bg=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,-0.5,1.5);
499 TH2D* h2Sig=new TH2D("rho_cc_out","rho_cc_out",50,0.,1.2,50,-0.5,1.5);
500 h2Bg->SetDirectory(0);
501 h2Sig->SetDirectory(0);
502 PUREZZA_S->SetDirectory(0);
503  double tt1=0;
504  double yy1=0;
505  double ty1=0;
506  double tmedio1=0;
507  double ymedio1=0;
508  double tt=0;
509  double yy=0;
510  double ty=0;
511  double tmedio=0;
512  double ymedio=0;
513  double rho1[NRHO];
514  double cc1[NCC];
515 int const CCRHO=NRHO*NCC;
516  double pur[CCRHO];
517  double tot[CCRHO];
518  for (int ii=0;ii<NRHO; ii++) {
519  rho1[ii]=0.;
520  for (int ji=0;ji<NCC; ji++) {
521  if (ii==0) cc1[ji]=0.;
522  pur[NRHO*ii+ji]=0.;
523  tot[NRHO*ii+ji]=0.;
524  }
525  }
526  double deltacc=0.;
527  //deltacc=1.2/NRHO;
528  deltacc=(maxCC-minCC)/NCC;
529  double deltarho=0.;
530  //deltarho=50/NRHO;
531  deltarho=(maxRHO-minRHO)/NRHO;
532  //for (int ii=0;ii<NRHO;ii++) { rho1[ii]=50.0-ii*deltarho;
533  for (int ii=0;ii<NRHO;ii++) {
534  rho1[ii]=maxRHO-(ii+1)*deltarho;
535 /* cout<<" rho_i "<<rho1[ii];
536  rho1[ii]=100*rho1[ii];
537  int rho0=0;
538  rho0=(int)(rho1[ii]+0.5);
539  rho1[ii]=0.;
540  rho1[ii]=(double)rho0/100;
541  cout<<" rho_f "<<rho1[ii]<<endl;
542 */ }
543  for (int ji=0;ji<NCC;ji++) {
544  cc1[ji]=(ji)*deltacc+minCC;
545  /*cout<<" cc_i "<<cc1[ji]<<endl;
546  cc1[ji]=100*cc1[ji];
547  int cc0=0;
548  cc0=(int)(cc1[ji]+0.5);
549  cc1[ji]=0.;
550  cc1[ji]=(double)cc0/100;
551  cout<<" cc_f "<<cc1[ji]<<endl;
552 */ }
553  for (int i = 0; i <(double)tTrain1->GetEntries(); i++) //o va +1 o +0,5?
554  {
555  tTrain1->GetEntry(i);
556  for (int j=0;j<nINP;j++)
557  {
558  params[j]=0;
559  params[j]=(double)x[j];
560  }
561  int j=npunti;
562  double output=mlp->Evaluate(0,params);
563  cout<<"output "<<output;
564  //<<endl;
565  int resto=i%2;
566  if (resto==0){
567  NNout=output;
568  if(type==1){
569  NNtype=1;
570  rho_NN=(double)signal.rho[0];
571  cc_NN=(double)signal.netcc[1];
572  }
573  if(type==0){
574  NNtype=0;
575  rho_NN=(double)background.rho[0];
576  cc_NN=(double)background.netcc[1];
577  }
578  t_NNout->Fill();
579  }
580 #ifdef RHO_CC
581  if (type==1&&resto==0){
582  sigTree.GetEntry(i+q);
583  h2SigO_R->Fill(signal.rho[0],output);
584  for (int ii=0; ii<NRHO;ii++){
585  if (signal.rho[0]>=rho1[ii]) {
586  for (int ji=0;ji<NCC;ji++){
587  if (signal.netcc[1]>=cc1[ji]) {
588  if (output>0.5){ pur[ii*NRHO+ji]=pur[ii*NRHO+ji]+1;
589  tot[ii*NRHO+ji]=tot[ii*NRHO+ji]+1;
590  }
591  }
592  }
593  }
594  }
595  if (output<0.1) cout<<"output<0.1, while type=1 event_index: "<<i+q<<" of file: "<<SIG_FILE.Data()<<endl;
596  h2SigR_C->Fill(signal.netcc[1],signal.rho[0]);
597  h2Sig->Fill(signal.netcc[1],output);
598  }
599  if (type==0&&resto==0){
600  bgTree.GetEntry(p-sig_entries+i-nTrain);
601  h2BgO_R->Fill(background.rho[0],output);
602  for (int ii=0; ii<NRHO;ii++){
603  for (int ji=0;ji<NCC;ji++){
604  if (background.rho[0]>rho1[ii]) {
605  if (background.netcc[1]>cc1[ji]) {
606  if (output>0.5) { tot[ii*NRHO+ji]=tot[ii*NRHO+ji]+1;
607  }
608  }
609  }
610  }
611  }
612  if (output>0.9) cout<<"output>0.9, while type=0 event_index: "<<p-sig_entries+i-nTrain<<" of file: "<<BG_FILE.Data()<<endl;
613  h2BgR_C->Fill(background.netcc[1],background.rho[0]);
614  h2Bg->Fill(background.netcc[1],output);
615  }
616 #endif
617  // cout<<"output "<<output<<" type "<<type;
618  // cout<<" #entry "<<i<<endl;
619  int entries=nTrain;
620  if(resto==0){
621  while(output<thres0[j]&&(j>0)) j=j-1;
622  for (int k=npunti;k>=0;k--){
623  if(type==0) {
624  if(k<=j) FA[k]=FA[k]+1;}
625  // cout<<"FA di "<<k<<" : "<<FA[k]<<" relativa a soglia "<<thres0[k]<<endl;}
626  if(type==1) {
627  if(k<=j) TA[k]=TA[k]+1; }
628  }
629  tmedio=tmedio+type;
630  ymedio=ymedio+output;
631  yy=output*output+yy;
632  tt=type*type+tt;
633  ty=type*output+ty;
634  }
635  if(resto!=0){
636  while(output<thres0[j]&&(j>0)) j=j-1;
637  for (int k=npunti;k>=0;k--){
638  if(type==0) {
639  if(k<=j) FA1[k]=FA1[k]+1;}
640  // cout<<"FA di "<<k<<" : "<<FA[k]<<" relativa a soglia "<<thres0[k]<<endl;}
641  if(type==1) {
642  if(k<=j) TA1[k]=TA1[k]+1; }
643  }
644  tmedio1=tmedio1+type;
645  ymedio1=ymedio1+output;
646  yy1=output*output+yy1;
647  tt1=type*type+tt1;
648  ty1=type*output+ty1;
649  }
650 
651  }
652 /* h2Bg->SetLineColor(kBlue);
653  h2Bg->SetFillStyle(3008); ibg->SetFillColor(kBlue);
654  h2Sig->SetLineColor(kRed);
655  h2Sig->SetFillStyle(3003); isig->SetFillColor(kRed);
656 */
657 
658  TString NNoutName(nomeNN);
659  NNoutName.ReplaceAll("NN/","NNout/NNout_");
660  TFile* f_NNout=new TFile(NNoutName,"RECREATE");
661  t_NNout->Write();
662  f_NNout->Close();
663 
664 cout<<""<<endl;
665  for (int ii=0;ii<NRHO; ii++) for (int ji=0;ji<NCC; ji++) {
666  cout<<" pur "<<pur[NRHO*ii+ji]<<" tot "<<tot[NRHO*ii+ji]<<" cc "<<cc1[ji]<< " rho "<<rho1[ii];
667 //<<endl;
668  if(tot[NRHO*ii+ji]>0) pur[NRHO*ii+ji]=pur[NRHO*ii+ji]/tot[NRHO*ii+ji];
669  cout<<" pur_fin "<<pur[NRHO*ii+ji];
670  //cout<<" index "<<NRHO*ii+ji;
671  cout<<""<<endl;
672  //if(tot[NRHO*ii+ji]>0) PUREZZA_S->Fill(cc1[ji],rho1[ii],pur[NRHO*ii+ji]);
673  if(tot[NRHO*ii+ji]>0) PUREZZA_S->Fill(cc1[ji]+deltacc/2,rho1[ii]+deltarho/2,pur[NRHO*ii+ji]);
674  }
675  h2Bg->SetMarkerStyle(7);
676  h2Sig->SetMarkerStyle(6);
677  h2Bg->SetMarkerColor(4);
678  h2Sig->SetMarkerColor(2);
679  h2BgO_R->SetMarkerStyle(7);
680  h2SigO_R->SetMarkerStyle(6);
681  h2BgO_R->SetMarkerColor(4);
682  h2SigO_R->SetMarkerColor(2);
683  h2BgR_C->SetMarkerStyle(7);
684  h2SigR_C->SetMarkerStyle(6);
685  h2BgR_C->SetMarkerColor(4);
686  h2SigR_C->SetMarkerColor(2);
687  // h2Bg->SetStats(0);
688  // h2Bg->SetStats(0);
689 
690  TCanvas *cRHO_CC_OUT=new TCanvas("RHO_CC_NNOUT","RHO_CC_NNOUT",1200,700);
691  //cRHO_CC_OUT->Divide(3,1);
692  cRHO_CC_OUT->Divide(2,2);
693  cRHO_CC_OUT->cd(1)->SetTitle("OUT_CC");
694  h2Bg->GetXaxis()->SetTitle("cc");
695  h2Bg->GetYaxis()->SetTitle("NNoutput");
696  h2Bg->GetZaxis()->SetTitle("count");
697  h2Sig->Draw();
698  h2Bg->Draw("same");
699  cRHO_CC_OUT->cd(2)->SetTitle("OUT_RHO");
700  h2BgO_R->GetXaxis()->SetTitle("rho");
701  h2BgO_R->GetYaxis()->SetTitle("NNoutput");
702  h2BgO_R->GetZaxis()->SetTitle("count");
703  h2SigO_R->Draw();
704  h2BgO_R->Draw("same");
705  cRHO_CC_OUT->cd(3)->SetTitle("RHO_CC");
706  h2BgR_C->GetXaxis()->SetTitle("cc");
707  h2BgR_C->GetYaxis()->SetTitle("rho");
708  h2BgR_C->GetZaxis()->SetTitle("count");
709  h2SigR_C->Draw();
710  h2BgR_C->Draw("same");
711  cRHO_CC_OUT->cd(4)->SetTitle("PUREZZA");
712  PUREZZA_S->SetStats(0);
713  PUREZZA_S->GetXaxis()->SetTitle("cc");
714  PUREZZA_S->GetYaxis()->SetTitle("rho");
715  PUREZZA_S->GetZaxis()->SetTitle("count");
716  PUREZZA_S->Draw("colz");
717 // PUREZZA_B->Draw("same");
718  TString cRHO_CC_OUTname(nomeNN);
719  TString cRHO_CC_OUTnameroot(nomeNN);
720  cRHO_CC_OUTname.ReplaceAll(".root",".png");
721  cRHO_CC_OUTname.ReplaceAll("NN/mlpNetwork","CC_OUT/RHO_CC_OUTgraph");
722  cRHO_CC_OUTnameroot.ReplaceAll("NN/mlpNetwork","CC_OUT/RHO_CC_OUTgraph");
723  cRHO_CC_OUT->SaveAs(cRHO_CC_OUTname.Data());
724  cRHO_CC_OUT->Print(cRHO_CC_OUTnameroot.Data());
725 // }
726  tmedio1=tmedio1/nTrain;
727  ymedio1=ymedio1/nTrain;
728 // tt1=tt1/nTrain;
729 // yy1=yy1/nTrain;
730 // ty1=ty1/nTrain;
731  tmedio=tmedio/nTrain;
732  ymedio=ymedio/nTrain;
733 // tt=tt/nTrain;
734 // yy=yy/nTrain;
735 // ty=ty/nTrain;
736  cout<<"tmedio1: "<<tmedio1<<" ymedio1 "<<ymedio1<<" tt1 "<<tt1<<" yy1 "<<yy1<<" ty1 "<<ty1<<" entries "<<entries<<endl;
737  double corr=0.;
738  double corr1=0.;
739  corr1=(ty1-nTrain*ymedio1*tmedio1)/sqrt((tt1-nTrain*tmedio1*tmedio1)*(yy1-nTrain*ymedio1*ymedio1));
740  corr=(ty-nTrain*ymedio*tmedio)/sqrt((tt-nTrain*tmedio*tmedio)*(yy-nTrain*ymedio*ymedio));
741  cout<<"corr "<<corr<<endl;
742  cout<<"corr1 "<<corr1<<endl;
743 
744 
745  for (int k=0;k<=npunti;k++) {
746  // cout<<"su "<<nTrain/2<<" TA: "<<TA[k]<<endl;
747  // cout<<"su "<<nTrainBg/2<<" FA: "<<FA[k]<<endl;
748  TA1[k]=(TA1[k]/nTrain)*2;
749  TA[k]=(TA[k]/nTrain)*2;
750  //cout<<"Eff: "<<TA[k]<<endl;
751  FA1[k]=(FA1[k]/nTrainBg)*2;
752  FA[k]=(FA[k]/nTrainBg)*2;
753  }
754  TString nf(nomeNN);
755  cout<<"ng "<<nf.Data()<<endl;
756  char nomefileTest[516];
757  int npunti2=npunti;
758  //sprintf(nomefile,"thresFiles/thres_np%i_nTS%i_nTB%i.root",npunti,nTrainS,nTrainB);
759  sprintf(nomefileTest,"thresFiles/thres_np%i_Test_onTrainSet_",npunti2);
760  nf.ReplaceAll("NN/mlpNetwork",nomefileTest);
761  cout<<"nf "<<nf.Data()<<endl;
762  TFile *gfile=new TFile(nf.Data(),"RECREATE");
763 
764 //#ifdef ROC_GRAPH
765  TGraph* gROC1=new TGraph(npunti+1,FA1,TA1);
766  gROC1->SetTitle("ROC_train(threshold)");
767  gROC1->SetMarkerStyle(7);
768  TGraph* gFA1=new TGraph(npunti+1,thres0,FA1);
769  gFA1->SetTitle("False Alarm _train(threshold)");
770  gFA1->SetMarkerStyle(7);
771  TGraph* gTA1=new TGraph(npunti+1,thres0,TA1);
772  gTA1->SetTitle("Efficiency _train(threshold)");
773  gTA1->SetMarkerStyle(7);
774 
775  TGraph* gROC=new TGraph(npunti+1,FA,TA);
776  gROC->SetTitle("ROC_test(threshold)");
777  gROC->SetMarkerStyle(7);
778  TGraph* gFA=new TGraph(npunti+1,thres0,FA);
779  gFA->SetTitle("False Alarm _test(threshold)");
780  gFA->SetMarkerStyle(7);
781  TGraph* gTA=new TGraph(npunti+1,thres0,TA);
782  gTA->SetTitle("Efficiency _test(threshold)");
783  gTA->SetMarkerStyle(7);
784 
785  TCanvas *cROC=new TCanvas("ROC_graph","ROC_graph");
786  cROC->Divide(2,3);
787  cROC->cd(2)->SetLogx();
788  // cROC->cd(2)->SetLogy();
789  gROC->Draw("apl");
790  cROC->cd(6);
791  cROC->cd(6)->SetLogy();
792  gFA->Draw("apl");
793  cROC->cd(4);
794  cROC->cd(4)->SetLogy();
795  gTA->Draw("apl");
796  cROC->cd(1)->SetLogx();
797  // cROC->cd(1)->SetLogy();
798  gROC1->Draw("apl");
799  cROC->cd(5);
800  cROC->cd(5)->SetLogy();
801  gFA1->Draw("apl");
802  cROC->cd(3);
803  cROC->cd(3)->SetLogy();
804  gTA1->Draw("apl");
805 
806  cROC->Write();
807  TString CANV_NAME(nomeNN);
808  TString CANV_NAMEroot(nomeNN);
809  CANV_NAME.ReplaceAll("NN/mlpNetwork","TRAIN_GRAPHS/ROC");
810  CANV_NAMEroot.ReplaceAll("NN/mlpNetwork","TRAIN_GRAPHS/ROC");
811  CANV_NAME.ReplaceAll(".root",".png");
812  cROC->SaveAs(CANV_NAME.Data());
813  cROC->Print(CANV_NAMEroot.Data());
814 //#endif
815 //#ifdef TH_TREE
816  TTree *ThTree=new TTree("ROC","ROC");
817 char FAlabel1[128];
818 char TAlabel1[128];
819 char FAlabel[128];
820 char TAlabel[128];
821 char thlabel[128];
822 sprintf(FAlabel1,"FalseAlarm_train[%i]/F",npunti);
823 sprintf(TAlabel1,"TrueAlarm_train[%i]/F",npunti);
824 sprintf(FAlabel,"FalseAlarm_test[%i]/F",npunti);
825 sprintf(TAlabel,"TrueAlarm_test[%i]/F",npunti);
826 sprintf(thlabel,"thresholds[%i]/F",npunti);
827 // float FAl=0;
828  ThTree->Branch("FalseAlarm_test",&FA,FAlabel);
829  ThTree->Branch("FalseAlarm_train",&FA1,FAlabel1);
830 // float TAl=0;
831  ThTree->Branch("Efficiency_train",&TA1,TAlabel1);
832  ThTree->Branch("Efficiency_test",&TA,TAlabel);
833 // double thresholds=0;
834  ThTree->Branch("thresholds",&thres0,thlabel);
835  ThTree->Branch("correlation_train",&corr1,"correlation_train/D");
836  ThTree->Branch("correlation_test",&corr,"correlation_test/D");
837 // ThTree->Branch("#TrainNN_Sig",&nTrainSNN,"#TrainNN_Sig/I");
838  int u=npunti;
839  ThTree->Branch("#points",&u,"#points/I");
840 // ThTree->Branch("#TrainNN_Bg",&nTrainBNN,"#TrainNN_Bg/I");
841  ThTree->Fill();
842  ThTree->Write();
843 //#endif
844 cout<<nomeNN<<endl;
845  gfile->Close();
846 }
847 
848 
849 
850 
851 TString select7Train_non_es(TString ORIGINAL_FILE,std::vector<int>* choice,TString nomeNN) {
852 TString IDname(nomeNN);
853 
854 //IDname.ReplaceAll("NN/","PropertyGraphs/");
855 
856 TFile* fOriginal=TFile::Open(ORIGINAL_FILE.Data());
857 TTree* tOriginal=(TTree*)fOriginal->Get("nnTree");
858 
859 if((choice->size())>(tOriginal->GetEntries())){cout<<"Error:startB+nTrainB>Entries"<<endl;
860  // exit(0);
861  }
862 
863 //estrazione nINP
864 int ninp;
865 tOriginal->SetBranchAddress("#inputs",&ninp);
866 tOriginal->GetEntry(0);
867 
868 int const nINP=ninp;
869 int const NDIM=sqrt(nINP);
870 
871 //definizione nuovo tree
872 TTree* t=new TTree("nnTree","nnTree");
873 
874 int type;
875 tOriginal->SetBranchAddress("type",&type);
876 t->Branch("type",&type,"type/I");
877 
878 float x[nINP];
879 char xlabel[nINP][16];
880 char xllabel[nINP][16];
881 #ifdef defferentW
882 float w[nINP];
883 char wlabel[nINP][16];
884 char wllabel[nINP][16];
885 #else
886 float w;
887 tOriginal->SetBranchAddress("w",&w);
888 t->Branch("w",&w,"w/F");
889 #endif
890 
891 #ifdef trainTEST
892 float SDs[nINP];
893 char SDlabels[nINP][16];
894 char SDllabels[nINP][16];
895 float RMSs[nINP];
896 char RMSlabels[nINP][16];
897 char RMSllabels[nINP][16];
898 float averages[nINP];
899 char averagelabels[nINP][16];
900 char averagellabels[nINP][16];
901 float SDb[nINP];
902 char SDlabelb[nINP][16];
903 char SDllabelb[nINP][16];
904 float RMSb[nINP];
905 char RMSlabelb[nINP][16];
906 char RMSllabelb[nINP][16];
907 float averageb[nINP];
908 char averagelabelb[nINP][16];
909 char averagellabelb[nINP][16];
910 float SD[nINP];
911 char SDlabel[nINP][16];
912 char SDllabel[nINP][16];
913 float RMS[nINP];
914 char RMSlabel[nINP][16];
915 char RMSllabel[nINP][16];
916 float average[nINP];
917 char averagelabel[nINP][16];
918 char averagellabel[nINP][16];
919 
920 TTree *tTrainTest=new TTree("trainTEST","trainTEST");
921 #endif
922 
923 for (int i=0; i<nINP;i++){
924  #ifdef trainTEST
925  RMS[i]=0;
926  average[i]=0;
927  SD[i]=0;
928  RMSs[i]=0;
929  averages[i]=0;
930  SDs[i]=0;
931  RMSb[i]=0;
932  averageb[i]=0;
933  SDb[i]=0;
934  #endif
935  sprintf(xlabel[i],"x%i",i+1);
936  sprintf(xllabel[i],"%s/F",xlabel[i]);
937  tOriginal->SetBranchAddress(xlabel[i],&x[i]);
938  t->Branch(xlabel[i],&x[i],xllabel[i]);
939  #ifdef differentW
940  sprintf(wlabel[i],"w%i",i+1);
941  sprintf(wllabel[i],"%s/F",wlabel[i]);
942  tOriginal->SetBranchAddress(wlabel[i],&w[i]);
943  t->Branch(wlabel[i],w[i],wllabel[i]);
944  #endif
945  }
946 int size=0;
947 size=choice->size();
948 int nTSB=0;
949 nTSB=size/2;
950 for (int k=0;k<choice->size();k++){
951  //cout<<"choice di "<<k<<" : "<<(*choice)[k]<<endl;
952  tOriginal->GetEntry((*choice)[k]);
953  #ifdef trainTEST
954  for (int i=0;i<nINP;i++){
955  RMS[i]=RMS[i]+x[i]*x[i];
956  average[i]=average[i]+x[i];
957  if(k<nTSB){
958  RMSs[i]=RMSs[i]+x[i]*x[i];
959  averages[i]=averages[i]+x[i];
960  }
961  if(k>=nTSB){
962  RMSb[i]=RMSb[i]+x[i]*x[i];
963  averageb[i]=averageb[i]+x[i];
964  }
965  }
966  #endif
967  for (int u=0;u<k;u++) if( (*choice)[u]==(*choice)[k]) {cout<<"Error:same event taken 2 times"<<"u "<<u<<" choise "<<(*choice)[u]<<endl; exit(0);}
968 // h[k]=(*choice)[k];
969  t->Fill();
970  }
971 /*
972 for (int k=b;k<b+nTrainB;k++){
973  tOriginal->GetEntry(k);
974  t->Fill();
975 
976 */
977 //cout<<t->GetEntries();
978 int Ntest=choice->size();
979 #ifdef trainTEST
980 /*TGraph *RMSg=new TGraph();
981 TGraph *AVg=new TGraph();
982 TGraph *SDg=new TGraph();
983 */
984 TH2D *RMShs=new TH2D("RMShs","RMShs",NDIM,0,NDIM,NDIM,0,NDIM);
985 TH2D *AVhs=new TH2D("AVhs","Avhs",NDIM,0,NDIM,NDIM,0,NDIM);
986 TH2D *SDhs=new TH2D("SDhs","SDhs",NDIM,0,NDIM,NDIM,0,NDIM);
987 TH2D *RMShb=new TH2D("RMShb","RMShb",NDIM,0,NDIM,NDIM,0,NDIM);
988 TH2D *AVhb=new TH2D("AVhb","Avhb",NDIM,0,NDIM,NDIM,0,NDIM);
989 TH2D *SDhb=new TH2D("SDhb","SDhb",NDIM,0,NDIM,NDIM,0,NDIM);
990 TH2D *RMSh=new TH2D("RMSh","RMSh",NDIM,0,NDIM,NDIM,0,NDIM);
991 TH2D *AVh=new TH2D("AVh","Avh",NDIM,0,NDIM,NDIM,0,NDIM);
992 TH2D *SDh=new TH2D("SDh","SDh",NDIM,0,NDIM,NDIM,0,NDIM);
993 //TH1D *RMSh=new TH1D("RMSh","RMSh",nINP,0,nINP-1);
994 //TH1D *AVh=new TH1D("AVh","Avh",nINP,0,nINP-1);
995 //TH1D *SDh=new TH1D("SDh","SDh",nINP,0,nINP-1);
996 RMSh->SetStats(0);
997 SDh->SetStats(0);
998 AVh->SetStats(0);
999 SDhb->SetStats(0);
1000 AVhb->SetStats(0);
1001 RMShb->SetStats(0);
1002 SDhs->SetStats(0);
1003 AVhs->SetStats(0);
1004 RMShs->SetStats(0);
1005 
1006 
1007 
1008 for (int i=0;i<nINP;i++){
1009  RMSs[i]=RMSs[i]/Ntest*2;
1010  averages[i]=averages[i]/Ntest*2;
1011  SDs[i]=RMSs[i]-pow(averages[i],2);
1012  SDs[i]=pow(SDs[i],0.5);
1013  RMSs[i]=pow(RMSs[i],0.5);
1014 
1015  cout << "SIGNAL: " << averages[i] << " " << RMSs[i] << " " << SDs[i] << endl;
1016 
1017  RMSb[i]=RMSb[i]/Ntest*2;
1018  averageb[i]=averageb[i]/Ntest*2;
1019  SDb[i]=RMSb[i]-pow(averageb[i],2);
1020  SDb[i]=pow(SDb[i],0.5);
1021  RMSb[i]=pow(RMSb[i],0.5);
1022 
1023  cout << "BKG: " << averageb[i] << " " << RMSb[i] << " " << SDb[i] << endl;
1024 
1025  RMS[i]=RMS[i]/Ntest;
1026  average[i]=average[i]/Ntest;
1027  SD[i]=RMS[i]-pow(average[i],2);
1028  SD[i]=pow(SD[i],0.5);
1029  RMS[i]=pow(RMS[i],0.5);
1030 
1031  cout << "ALL: " << average[i] << " " << RMS[i] << " " << SD[i] << endl;
1032 
1033  sprintf(SDlabels[i],"SDs%i",i+1);
1034  sprintf(SDllabels[i],"%s/F",SDlabels[i]);
1035  tTrainTest->Branch(SDlabels[i],&SDs[i],SDllabels[i]);
1036  sprintf(RMSlabels[i],"RMSs%i",i+1);
1037  sprintf(RMSllabels[i],"%s/F",RMSlabels[i]);
1038  tTrainTest->Branch(RMSlabels[i],&RMSs[i],RMSllabels[i]);
1039  sprintf(averagelabels[i],"averages%i",i+1);
1040  sprintf(averagellabels[i],"%s/F",averagelabels[i]);
1041  tTrainTest->Branch(averagelabels[i],&averages[i],averagellabels[i]);
1042  sprintf(SDlabelb[i],"SDb%i",i+1);
1043  sprintf(SDllabelb[i],"%s/F",SDlabelb[i]);
1044  tTrainTest->Branch(SDlabelb[i],&SDb[i],SDllabelb[i]);
1045  sprintf(RMSlabelb[i],"RMSb%i",i+1);
1046  sprintf(RMSllabelb[i],"%s/F",RMSlabelb[i]);
1047  tTrainTest->Branch(RMSlabelb[i],&RMSb[i],RMSllabelb[i]);
1048  sprintf(averagelabelb[i],"averageb%i",i+1);
1049  sprintf(averagellabelb[i],"%s/F",averagelabelb[i]);
1050  tTrainTest->Branch(averagelabelb[i],&averageb[i],averagellabelb[i]);
1051  sprintf(SDlabel[i],"SD%i",i+1);
1052  sprintf(SDllabel[i],"%s/F",SDlabel[i]);
1053  tTrainTest->Branch(SDlabel[i],&SD[i],SDllabel[i]);
1054  sprintf(RMSlabel[i],"RMS%i",i+1);
1055  sprintf(RMSllabel[i],"%s/F",RMSlabel[i]);
1056  tTrainTest->Branch(RMSlabel[i],&RMS[i],RMSllabel[i]);
1057  sprintf(averagelabel[i],"average%i",i+1);
1058  sprintf(averagellabel[i],"%s/F",averagelabel[i]);
1059  tTrainTest->Branch(averagelabel[i],&average[i],averagellabel[i]);
1060  /* RMSg->SetPoint(i,i,RMS[i]);
1061  AVg->SetPoint(i,i,average[i]);
1062  SDg->SetPoint(i,i,SD[i]);
1063  */
1064  int fi=0;
1065  int ti=0;
1066  ti=i/NDIM;
1067  fi=i-ti*NDIM;
1068 // cout<<"ti: "<<ti<<" fi: "<<fi<<endl;
1069  RMShs->Fill(ti,fi,RMSs[i]);
1070  AVhs->Fill(ti,fi,averages[i]);
1071  SDhs->Fill(ti,fi,SDs[i]);
1072  RMShb->Fill(ti,fi,RMSb[i]);
1073  AVhb->Fill(ti,fi,averageb[i]);
1074  SDhb->Fill(ti,fi,SDb[i]);
1075  RMSh->Fill(ti,fi,RMS[i]);
1076  AVh->Fill(ti,fi,average[i]);
1077  SDh->Fill(ti,fi,SD[i]);
1078 // RMSh->Fill(i,RMS[i]);
1079 // AVh->Fill(i,average[i]);
1080 // SDh->Fill(i,SD[i]);
1081  }
1082 TCanvas *cProp=new TCanvas ("TrainSet_Properties","TrainSet_Properties",0,0,900,600);
1083 cProp->Divide(3,3);
1084 cProp->cd(1);
1085 RMSh->SetTitle("RMS");
1086 RMSh->Draw("COLZ");
1087 cProp->cd(2);
1088 AVh->SetTitle("Average");
1089 AVh->Draw("COLZ");
1090 cProp->cd(3);
1091 SDh->SetTitle("SD");
1092 SDh->Draw("COLZ");
1093 cProp->cd(4);
1094 RMShs->SetTitle("RMS_s");
1095 RMShs->Draw("COLZ");
1096 cProp->cd(5);
1097 AVhs->SetTitle("Average_s");
1098 AVhs->Draw("COLZ");
1099 cProp->cd(6);
1100 SDhs->SetTitle("SD_s");
1101 SDhs->Draw("COLZ");
1102 cProp->cd(7);
1103 RMShb->SetTitle("RMS_b");
1104 RMShb->Draw("COLZ");
1105 cProp->cd(8);
1106 AVhb->SetTitle("Average_b");
1107 AVhb->Draw("COLZ");
1108 cProp->cd(9);
1109 SDhb->SetTitle("SD_b");
1110 SDhb->Draw("COLZ");
1111 
1112 TString Prop_name(IDname);
1113 TString Prop_nameroot(IDname);
1114 Prop_name.ReplaceAll(".root","Prop.png");
1115 Prop_nameroot.ReplaceAll(".root","Prop.root");
1116 cProp->SaveAs(Prop_name.Data());
1117 cProp->Print(Prop_nameroot.Data());
1118 cProp->Close();
1119 #endif
1120 int ent=t->GetEntries();
1121 if (choice->size()!=ent) cout<<"Error in Tree dimension: Entries new tree="<<t->GetEntries()<<" vector->size"<<choice->size()<<endl;
1122 
1123 //char nomefile[516];
1124 //sprintf(nomefile,"selectTREE/TreeTrain_size%i_s%i_b%i.root",choice->size(),s,b);
1125 TString NOME_FILE(IDname);
1126 NOME_FILE.ReplaceAll("PropertyGraphs/mlpNetwork","selectTREE/TreeTrain");
1127 TFile* ofile=new TFile(NOME_FILE.Data(),"RECREATE");
1128 t->SetDirectory(ofile);
1129 t->Write();
1130 ofile->Close();
1131 fOriginal->Close();
1132 
1133 #ifdef trainTEST
1134 /*TCanvas *matrixIndex_RMS=new TCanvas("matrixIndex_RMS","matrixIndex_RMS",0,0,600,600);
1135 tTrainTest->Draw("RMS");
1136 matrixIndex*/
1137 TFile* TrainTest= new TFile(IDname.Data(),"RECREATE");
1138 tTrainTest->Write();
1139 TrainTest->Close();
1140 #endif
1141 cout<<NOME_FILE.Data()<<endl;
1142 return NOME_FILE;
1143 }
1144 
1145 
t
wavearray< double > t(hp.size())
cwbANN_purity_BoS
void cwbANN_purity_BoS(int nTrain, int nTrainBg, TString option, TString TREE_FILE, int q, int p, double rm=0., double cm=0., int lm=5, Int_t nepoch=700)
Definition: cwbANN_purity_BoS.C:44
minCC
#define minCC
Definition: Average_FAout_moreGraphs_3.C:26
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:93
select7Train_non_es
TString select7Train_non_es(TString ORIGINAL_FILE, std::vector< int > *choice, TString IDname)
Definition: cwbANN_purity_BoS.C:851
z
wavearray< double > z
Definition: Test10.C:32
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
NCC
#define NCC
Definition: cwbANN_purity_BoS.C:32
size
Long_t size
Definition: cwb_condor_check.C:48
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
w
wavearray< double > w
Definition: Test1.C:27
NDIM
#define NDIM
Definition: ClusterToFrameNew01_dtmin_approx_cent.C:31
h
wavearray< double > h
Definition: Regression_H1.C:25
log
bool log
Definition: WaveMDC.C:41
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
NRHO
#define NRHO
Definition: cwbANN_purity_BoS.C:31
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
add
r add(tfmap,"hchannel")
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
yy
wavearray< double > yy
Definition: TestFrame5.C:12
if
if[[-z $1]]
Definition: cbc_plots.sh:1
deltarho
#define deltarho
Definition: Average_FAout_moreGraphs_3.C:32
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
nIFO
#define nIFO
Definition: cwbANN_purity_BoS.C:33
gfile
TString gfile
Definition: DrawPhaseShift.C:44
netevent::GetEntries
Int_t GetEntries()
Definition: netevent.cc:381
y
wavearray< double > y
Definition: Test10.C:31
deltacc
#define deltacc
Definition: Average_FAout_moreGraphs_3.C:31
exit
exit(0)
maxCC
#define maxCC
Definition: Average_FAout_moreGraphs_3.C:25