mlpPIDelecallEetot.C
Go to the documentation of this file.
1 void mlpPIDelecallEetot(Int_t ntrain= 600 ) {
2 //Author: Jianming Bian
3  double Elow = 1.0;
4  double Ehigh = 3.0;
5  double rsig = 1.0;
6  double rbg = 0.6;
7 
8  TChain signal("showerrec/fEvent");
9  TChain background("showerrec/fEvent");
10 
11  int nf = 0;
12  for(int i=0;i<600;i++){
13  char fname[200];
14  sprintf(fname,"/swap/caf_swap_fd_nhc_%d.root",i);
15  ifstream fin(fname);
16  if(!fin){cout<<fname<<" "<<"doesn't exist"<<endl;continue;}
17  signal.Add(fname);
18  nf++;cout<<"nf = "<<nf<<endl;
19  if(nf==501)break;
20  }
21  int nfb = 0;
22  for(int i=0;i<600;i++){
23  char fname[200];
24 // sprintf(fname,"/export/data/local/novadata/nuevt/mc/S12-10-07/FD/newhistnewfkcaf_fd_nhc_%d.root",i);
25  sprintf(fname,"/nonswap/ptcaf_fd_nhc_%d.root",i);
26  ifstream fin(fname);
27  if(!fin){cout<<fname<<" "<<"doesn't exist"<<endl;continue;}
28  background.Add(fname);
29  nfb++;cout<<"nfb = "<<nfb<<endl;
30  if(nfb==501)break;
31  }
32 
33 
34 
35  TCut MCNUECCQE = "evtTrueNuPdg==12&&(evtRwRand[0]>0&&evtRwRand[0]<evtRwProb[0])";
36  TCut MCBKG = "(evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])";
37  TH1F* htruenue=new TH1F("htruenue","htruenue", 100,0, 5);
38  TH1F* htruenupi0=new TH1F("htruenupi0","htruenupi0", 100,0, 5);
39  signal.Draw("evtTrueNuEnergy>>htruenue",MCNUECCQE);
40  background.Draw("evtTrueNuEnergy>>htruenupi0", MCBKG);
41 
42  double nsig = htruenue->GetEntries();
43  double nbkg = htruenupi0->GetEntries();
44 
45  TCanvas *c0 = new TCanvas("c0","c0",200,10, 1200,600);
46  c0->Divide(2,1);
47  c0->cd(1);
48  htruenue->Draw();
49  c0->cd(2);
50  htruenupi0->Draw();
51  c0->Update();
52 
53 
54  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
55  gSystem->Load("libMLP");
56  }
57  gSystem->Load("libRooFit");
58 
59  int evtTrueNuCCNC;
60  int evtTrueNuMode;
61  int evtTrueNuPdg;
62  double evtTrueNuEnergy;
63  int evtTrueNuIsFid;
64  double evtTrueNuP4[4];
65  double evtTrueNuVtx[3];
66 
67  int evtOsTag;
68  double evtOsProb[20];
69  double evtOsRand;
70 
71  double evtRwProb[20];
72  double evtRwRand[20];
73 
74  int evtTruePdg[20];
75  int evtTrueMother[20];
76  double evtTrueEnergy[20];
77 
78  int evtRun;
79  int evtSRun;
80  int evtEvent;
81  int evtNSlice;
82  int evtNTrk;
83  double evtPi0Mgg0;
84  double evtPi0Mgg;
85  double evtSh1Pi0Mgg;
86  double evtSh1VtxGeV;
87  double evtSh1CosTheta;
88  double evtSh1Length;
89  double evtPi0Eg1;
90  double evtPi0Eg2;
91  double evtPi0TrueDirMgg0;
92  double evtPi0TrueDirMgg;
93  double evtPi0StartDist;
94  double evtPi0LineDist;
95 
96  double evtEtot0;
97  double evtEtot;
98  double evtSh1SumPt;
99 
100  double evtPEtot;
101  double evtTrueEtot;
102  double evtSh1Enue;
103 
104  double evtSh1Energy;
105  double evtSh1NueEnergy;
106  double evtSh1Energy0;
107  double evtTotalRecoGeV;
108  double evtSh1SliceGeV;
109  double evtSheEnergy;
110  double evtSh1Gap;
111  double evtSh1IsFid;
112  double evtSh1P4[4];
113  double evtSh1Start[3];
114  double evtSh1DedxLLL[20];
115  double evtSh1DedxLLT[20];
116  double evtSh1NPlane;
117  double evtSh1NMIPPlane;
118  double evtSh1ContStartPlane;
119  double evtSh2Energy;
120 
121  int evtSh1TruePdg;
122  double evtSh1TrueDang;
123  double evtSh1TrueEDang;
124 
125  Int_t type;
126 
127 
128  signal.SetBranchAddress("evtTrueNuCCNC", &evtTrueNuCCNC);
129  signal.SetBranchAddress("evtTrueNuMode", &evtTrueNuMode);
130  signal.SetBranchAddress("evtTrueNuPdg", &evtTrueNuPdg);
131  signal.SetBranchAddress("evtTrueNuEnergy", &evtTrueNuEnergy);
132  signal.SetBranchAddress("evtTrueNuIsFid", &evtTrueNuIsFid);
133  signal.SetBranchAddress("evtTrueNuP4[4]", evtTrueNuP4);
134  signal.SetBranchAddress("evtTrueNuVtx[3]",evtTrueNuVtx);
135  signal.SetBranchAddress("evtOsTag", &evtOsTag);
136  signal.SetBranchAddress("evtOsProb[20]",evtOsProb);
137  signal.SetBranchAddress("evtOsRand", &evtOsRand);
138 
139  signal.SetBranchAddress("evtRwProb[20]",evtRwProb);
140  signal.SetBranchAddress("evtRwRand[20]",evtRwRand);
141 
142  signal.SetBranchAddress("evtTruePdg[20]",evtTruePdg);
143  signal.SetBranchAddress("evtTrueMother[20]",evtTrueMother);
144  signal.SetBranchAddress("evtTrueEnergy[20]",evtTrueEnergy);
145 
146  signal.SetBranchAddress("evtSh1TruePdg",&evtSh1TruePdg);
147  signal.SetBranchAddress("evtSh1TrueDang",&evtSh1TrueDang);
148  signal.SetBranchAddress("evtSh1TrueEDang",&evtSh1TrueEDang);
149  signal.SetBranchAddress("evtSh1Energy",&evtSh1Energy);
150  signal.SetBranchAddress("evtSh1NueEnergy",&evtSh1NueEnergy);
151  signal.SetBranchAddress("evtSh1Energy0",&evtSh1Energy0);
152  signal.SetBranchAddress("evtTotalRecoGeV",&evtTotalRecoGeV);
153  signal.SetBranchAddress("evtSh1SliceGeV",&evtSh1SliceGeV);
154  signal.SetBranchAddress("evtSheEnergy",&evtSheEnergy);
155  signal.SetBranchAddress("evtSh1Gap",&evtSh1Gap);
156  signal.SetBranchAddress("evtSh1IsFid",&evtSh1IsFid);
157  signal.SetBranchAddress("evtSh1P4[4]",evtSh1P4);
158  signal.SetBranchAddress("evtSh1Start[3]",evtSh1Start);
159  signal.SetBranchAddress("evtSh1DedxLLL[20]",evtSh1DedxLLL);
160  signal.SetBranchAddress("evtSh1DedxLLT[20]",evtSh1DedxLLT);
161  signal.SetBranchAddress("evtSh1NPlane",&evtSh1NPlane);
162  signal.SetBranchAddress("evtSh1NMIPPlane",&evtSh1NMIPPlane);
163  signal.SetBranchAddress("evtSh1ContStartPlane",&evtSh1ContStartPlane);
164  signal.SetBranchAddress("evtEtot",&evtEtot);
165  signal.SetBranchAddress("evtSh1SumPt",&evtSh1SumPt);
166  signal.SetBranchAddress("evtNTrk",&evtNTrk);
167  signal.SetBranchAddress("evtSh1Enue",&evtSh1Enue);
168  signal.SetBranchAddress("evtSh2Energy",&evtSh2Energy);
169  signal.SetBranchAddress("evtPi0Mgg",&evtPi0Mgg);
170  signal.SetBranchAddress("evtSh1Pi0Mgg",&evtSh1Pi0Mgg);
171  signal.SetBranchAddress("evtSh1VtxGeV",&evtSh1VtxGeV);
172  signal.SetBranchAddress("evtSh1CosTheta",&evtSh1CosTheta);
173  signal.SetBranchAddress("evtSh1Length",&evtSh1Length);
174  signal.SetBranchAddress("evtPi0StartDist",&evtPi0StartDist);
175  signal.SetBranchAddress("evtPi0LineDist",&evtPi0LineDist);
176 
177 
178 
179  background.SetBranchAddress("evtTrueNuCCNC", &evtTrueNuCCNC);
180  background.SetBranchAddress("evtTrueNuMode", &evtTrueNuMode);
181  background.SetBranchAddress("evtTrueNuPdg", &evtTrueNuPdg);
182  background.SetBranchAddress("evtTrueNuEnergy", &evtTrueNuEnergy);
183  background.SetBranchAddress("evtTrueNuIsFid", &evtTrueNuIsFid);
184  background.SetBranchAddress("evtTrueNuP4[4]", evtTrueNuP4);
185  background.SetBranchAddress("evtTrueNuVtx[3]",evtTrueNuVtx);
186  background.SetBranchAddress("evtOsTag", &evtOsTag);
187  background.SetBranchAddress("evtOsProb[20]",evtOsProb);
188  background.SetBranchAddress("evtOsRand", &evtOsRand);
189 
190  background.SetBranchAddress("evtRwProb[20]",evtRwProb);
191  background.SetBranchAddress("evtRwRand[20]",evtRwRand);
192 
193  background.SetBranchAddress("evtTruePdg[20]",evtTruePdg);
194  background.SetBranchAddress("evtTrueMother[20]",evtTrueMother);
195  background.SetBranchAddress("evtTrueEnergy[20]",evtTrueEnergy);
196 
197  background.SetBranchAddress("evtSh1TruePdg",&evtSh1TruePdg);
198  background.SetBranchAddress("evtSh1TrueDang",&evtSh1TrueDang);
199  background.SetBranchAddress("evtSh1TrueEDang",&evtSh1TrueEDang);
200  background.SetBranchAddress("evtSh1Energy",&evtSh1Energy);
201  background.SetBranchAddress("evtSh1NueEnergy",&evtSh1NueEnergy);
202  background.SetBranchAddress("evtSh1Energy0",&evtSh1Energy0);
203  background.SetBranchAddress("evtTotalRecoGeV",&evtTotalRecoGeV);
204  background.SetBranchAddress("evtSh1SliceGeV",&evtSh1SliceGeV);
205  background.SetBranchAddress("evtSheEnergy",&evtSheEnergy);
206  background.SetBranchAddress("evtSh1Gap",&evtSh1Gap);
207  background.SetBranchAddress("evtSh1IsFid",&evtSh1IsFid);
208  background.SetBranchAddress("evtSh1P4[4]",evtSh1P4);
209  background.SetBranchAddress("evtSh1Start[3]",evtSh1Start);
210  background.SetBranchAddress("evtSh1DedxLLL[20]",evtSh1DedxLLL);
211  background.SetBranchAddress("evtSh1DedxLLT[20]",evtSh1DedxLLT);
212  background.SetBranchAddress("evtSh1NPlane",&evtSh1NPlane);
213  background.SetBranchAddress("evtSh1NMIPPlane",&evtSh1NMIPPlane);
214  background.SetBranchAddress("evtSh1ContStartPlane",&evtSh1ContStartPlane);
215  background.SetBranchAddress("evtEtot",&evtEtot);
216  background.SetBranchAddress("evtSh1SumPt",&evtSh1SumPt);
217  background.SetBranchAddress("evtNTrk",&evtNTrk);
218  background.SetBranchAddress("evtSh1Enue",&evtSh1Enue);
219  background.SetBranchAddress("evtSh2Energy",&evtSh2Energy);
220  background.SetBranchAddress("evtPi0Mgg",&evtPi0Mgg);
221  background.SetBranchAddress("evtSh1Pi0Mgg",&evtSh1Pi0Mgg);
222  background.SetBranchAddress("evtSh1VtxGeV",&evtSh1VtxGeV);
223  background.SetBranchAddress("evtSh1CosTheta",&evtSh1CosTheta);
224  background.SetBranchAddress("evtSh1Length",&evtSh1Length);
225  background.SetBranchAddress("evtPi0StartDist",&evtPi0StartDist);
226  background.SetBranchAddress("evtPi0LineDist",&evtPi0LineDist);
227 
228 
229  TTree *simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events");
230  double egLLL, egLLT, emuLLL, emuLLT, epi0LLL, epi0LLT, epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT,gap, pi0mass;
231  double vtxgev;
232  double shE;
233  double nuE;
234  double costheta;
235  double length;
236  int ntrk;
237 
238  simu->Branch("egLLL", &egLLL, "egLLL/D");
239  simu->Branch("egLLT", &egLLT, "egLLT/D");
240  simu->Branch("emuLLL", &emuLLL, "emuLLL/D");
241  simu->Branch("emuLLT", &emuLLT, "emuLLT/D");
242  simu->Branch("epi0LLL", &epi0LLL, "epi0LLL/D");
243  simu->Branch("epi0LLT", &epi0LLT, "epi0LLT/D");
244  simu->Branch("epLLL", &epLLL, "epLLL/D");
245  simu->Branch("epLLT", &epLLT, "epLLT/D");
246  simu->Branch("enLLL", &enLLL, "enLLL/D");
247  simu->Branch("enLLT", &enLLT, "enLLT/D");
248  simu->Branch("epiLLL", &epiLLL, "epiLLL/D");
249  simu->Branch("epiLLT", &epiLLT, "epiLLT/D");
250  simu->Branch("gap", &gap, "gap/D");
251  simu->Branch("pi0mass", &pi0mass, "pi0mass/D");
252  simu->Branch("vtxgev", &vtxgev, "vtxgev/D");
253  simu->Branch("shE", &shE, "shE/D");
254  simu->Branch("nuE", &nuE, "nuE/D");
255 // simu->Branch("ntrk", &ntrk, "ntrk/I");
256  simu->Branch("costheta", &costheta, "costheta/D");
257  simu->Branch("length", &length, "length/D");
258  simu->Branch("type", &type, "type/I");
259 
260  TTree *annMC = new TTree("annMC", "Filtered Monte Carlo Events");
261  double trueNuE, etot, enue, sh1Energy, ANN, wSig, wBkg;
262  double annPi0StartDist, annPi0LineDist;
263  annMC->Branch("egLLL", &egLLL, "egLLL/D");
264  annMC->Branch("egLLT", &egLLT, "egLLT/D");
265  annMC->Branch("emuLLL", &emuLLL, "emuLLL/D");
266  annMC->Branch("emuLLT", &emuLLT, "emuLLT/D");
267  annMC->Branch("epi0LLL", &epi0LLL, "epi0LLL/D");
268  annMC->Branch("epi0LLT", &epi0LLT, "epi0LLT/D");
269  annMC->Branch("epLLL", &epLLL, "epLLL/D");
270  annMC->Branch("epLLT", &epLLT, "epLLT/D");
271  annMC->Branch("enLLL", &enLLL, "enLLL/D");
272  annMC->Branch("enLLT", &enLLT, "enLLT/D");
273  annMC->Branch("epiLLL", &epiLLL, "epiLLL/D");
274  annMC->Branch("epiLLT", &epiLLT, "epiLLT/D");
275  annMC->Branch("gap", &gap, "gap/D");
276  annMC->Branch("pi0mass", &pi0mass, "pi0mass/D");
277  annMC->Branch("vtxgev", &vtxgev, "vtxgev/D");
278  annMC->Branch("length", &length, "length/D");
279  annMC->Branch("annPi0StartDist", &annPi0StartDist, "annPi0StartDist/D");
280  annMC->Branch("annPi0LineDist", &annPi0LineDist, "annPi0LineDist/D");
281  annMC->Branch("evtTrueNuCCNC", &evtTrueNuCCNC, "evtTrueNuCCNC/I");
282  annMC->Branch("evtTrueNuPdg", &evtTrueNuPdg, "evtTrueNuPdg/I");
283  annMC->Branch("evtTrueNuMode", &evtTrueNuMode, "evtTrueNuMode/I");
284  annMC->Branch("trueNuE", &trueNuE, "trueNuE/D");
285  annMC->Branch("etot", &etot, "etot/D");
286  annMC->Branch("enue", &enue, "enue/D");
287  annMC->Branch("sh1Energy", &sh1Energy, "sh1Energy/D");
288  annMC->Branch("type", &type, "type/I");
289  annMC->Branch("wSig", &wSig, "wSig/D");
290  annMC->Branch("wBkg", &wBkg, "wBkg/D");
291  annMC->Branch("ANN", &ANN, "ANN/D");
292 
293  int nb0=0,nb1=0,nb2=0,nb3=0;
294  int nbNueCC0=0,nbNueCC1=0,nbNueCC2=0,nbNueCC2_a=0,nbNueCC3=0;
295  int nbNumuCC0=0,nbNumuCC1=0,nbNumuCC2=0,nbNumuCC3=0;
296  int nbNC0=0,nbNC1=0,nbNC2=0,nbNC3=0;
297 
298  int ns0=0,ns1=0,ns2=0,ns3=0;
299  int normSig = 0, normBkg =0, totBkg=0;
300  type = 1;
301 
302  for (int i = 0; i < signal.GetEntries(); i++) {
303  signal.GetEntry(i);
304  if(!(evtTrueNuCCNC==0&&evtTrueNuPdg==12&&(evtRwRand[0]>0&&evtRwRand[0]<evtRwProb[0])))continue;
305  if(i%2==1)ns0++;
306 // if(evtSh1Energy!=evtSheEnergy)continue;
307 // if(fabs(evtSh1TruePdg)!=11)continue;
308  if(evtSh1TrueDang>30||fabs(evtSh1TruePdg)!=11)continue;
309 // if(!(evtSh1TrueEDang<30&&evtSh1TrueEDang>0))continue;
310 // if(!(evtSh1IsFid==1&&evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
311  if(!(evtSh1Energy>0.0))continue;
312  if(i%2==1)ns1++;
313 // if(evtSh1DedxLLL[0]<0)continue;
314  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
315 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
316 // if(!(evtSh1Enue/1.40694>1.0&&evtSh1Enue/1.40694<Ehigh))continue;
317 // if(!((evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
318 // if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
319 
320  if(i%2==1)ns2++;
321 
322  if((double)i>(double)signal.GetEntries()*rsig)continue;
323 
324  egLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
325  egLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
326  emuLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
327  emuLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
328  epi0LLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
329  epi0LLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
330  epLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
331  epLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
332  enLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
333  enLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
334  epiLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
335  epiLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
336  gap = evtSh1Gap;
337  pi0mass = MAX(evtSh1Pi0Mgg,0);
338  shE = evtSh1Energy/evtEtot;
339  nuE = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
340  ntrk = evtNTrk;
341 // shE = evtSh1Energy*5.63616e-01/evtSh1SliceGeV;
342  vtxgev = evtSh1VtxGeV;
343 // costheta = evtSh1SumPt/evtSh1SliceGeV;
344  costheta = evtSh1CosTheta;
345  length = evtSh1Length/evtSh1Energy;
346  simu->Fill();
347  normSig++;
348  }
349 
350  type = 0;
351  for (int i = 0; i < background.GetEntries(); i++) {
352 // if(i>36383)continue;
353  background.GetEntry(i);
354  if(!((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))continue;
355  if(i%2==1)nb0++;
356 
357  if(!((fabs(evtTrueNuPdg)==14)||(fabs(evtTrueNuPdg)==12)))continue;
358  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])nbNumuCC0++;
359  if(i%2==1&&evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14)||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))nbNC0++;
360  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])nbNueCC0++;
361 
362  if(!(evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
363  if(i%2==1)nb1++;
364 
365  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])nbNumuCC1++;
366  if(i%2==1&&evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))nbNC1++;
367  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])nbNueCC1++;
368 
369 
370  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
371 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
372 // if(evtSh1DedxLLL[0]<0)continue;
373 // if(!(evtSh1Enue/1.40694>Elow&&evtSh1Enue/1.40694<Ehigh))continue;
374 
375  if(i%2==1)nb2++;
376  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])nbNumuCC2++;
377  if(i%2==1&&evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))nbNC2++;
378  if(i%2==1&&evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])nbNueCC2++;
379 // if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12)continue;
380 // nb3++;
381 
382  if((double)i>(double)background.GetEntries()*rbg)continue;
383 
384  egLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
385  egLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
386  emuLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
387  emuLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
388  epi0LLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
389  epi0LLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
390  epLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
391  epLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
392  enLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
393  enLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
394  epiLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
395  epiLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
396  gap = evtSh1Gap;
397  pi0mass = MAX(evtSh1Pi0Mgg,0);
398  vtxgev = evtSh1VtxGeV;
399 // costheta = evtSh1SumPt/evtSh1SliceGeV;
400  costheta = evtSh1CosTheta;
401  length = evtSh1Length/evtSh1Energy;
402  shE = evtSh1Energy/evtEtot;
403  nuE = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
404  ntrk = evtNTrk;
405  simu->Fill();
406  normBkg++;
407  }
408 
409  cout<<"Nsig, Nbig, totBkg, scaleBkg "<<normSig<<" "<<normBkg<<" "<<endl;
410 
411  double scalepotsig = 18.0e+20/((ns0/(2.87128765282905034e-02*6938.035))*2.99649003626026520e+21);//2.38399959425089536e+19);//NDOS
412  double scalepotNumuCC = 18.0e+20/((nbNumuCC0/(6938.035-199.21-5511.0658))*3.02657807308939830e+21);//2.38399959425089536e+19);//NDOS
413  double scalepotNueCC = 18.0e+20/((nbNueCC0/82.02)*3.02657807308939830e+21);//2.38399959425089536e+19);//NDOS
414  double scalepotNC = 18.0e+20/((nbNC0/2664.42968)*3.02657807308939830e+21);//2.38399959425089536e+19);//NDOS
415  cout<<"nbNumuCC0, scalepotNumuCC "<<nbNumuCC0<<" "<<scalepotNumuCC<<endl;
416 
417  // Build and train the NN ptsumf is used as a weight since we are primarly
418  // interested by high pt events.
419  // The datasets used here are the same as the default ones.
420 /*
421  TMultiLayerPerceptron *mlp =
422  new TMultiLayerPerceptron("epi0LLL,epi0LLT,epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT:15:3:type",
423  simu,"Entry$%2","(Entry$+1)%2");
424 */
425 /*
426  TMultiLayerPerceptron *mlp =
427  new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@vtxgev,@gap:18:12:6:type", simu,"Entry$%2","(Entry$+1)%2");
428 */
429 
430  TMultiLayerPerceptron *mlp =
431  new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@costheta,@nuE:22:12:6:type", simu,"Entry$%2","(Entry$+1)%2");
432 
433 
434 
435 /*
436  TMultiLayerPerceptron *mlp =
437  new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@vtxgev,@gap:22:14:8:type", simu,"Entry$%2","(Entry$+1)%2");
438 */
439 /*
440  TMultiLayerPerceptron *mlp =
441  new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@gap:18:12:6:type", simu,"Entry$%2","(Entry$+1)%2");
442 */
443  mlp->Train(ntrain, "text,graph,update=10");
444 // mlp->DumpWeights("weights_shape_fd_elec_Enewfk.txt");
445  mlp->DumpWeights("weights_shape_fd_elec_E_noecut.txt");
446  mlp->Export("test","python");
447  // Use TMLPAnalyzer to see what it looks for
448  TCanvas* mlpa_canvas1 = new TCanvas("mlpa_canvas1","Network analysis 1", 10,10, 1200,600);
449  mlpa_canvas1->Divide(2,1);
450  TMLPAnalyzer ana(mlp);
451  // Initialisation
452  ana.GatherInformations();
453  // output to the console
454  ana.CheckNetwork();
455  mlpa_canvas1->cd(1);
456  // shows how each variable influences the network
457  ana.DrawDInputs();
458  mlpa_canvas1->cd(2);
459  // shows the network structure
460  mlp->Draw();
461  mlpa_canvas1->Update();
462 
463  TCanvas* mlpa_canvas2 = new TCanvas("mlpa_canvas2","Network analysis 2", 10,10, 1200,600);
464  mlpa_canvas2->Divide(2,1);
465  mlpa_canvas2->cd(1);
466  // draws the resulting network
467  ana.DrawNetwork(0,"type==1","type==0");
468 
469  mlpa_canvas2->cd(2);
470  // Use the NN to plot the results for each sample
471  // This will give approx. the same result as DrawNetwork.
472  // All entries are used, while DrawNetwork focuses on
473  // the test sample. Also the xaxis range is manually set.
474  TH1F *bg = new TH1F("bgh", "NN output", 50, -.5, 1.5);
475  TH1F *bgNumuCC = new TH1F("bgNumuCC", "NN output", 50, -.5, 1.5);
476  TH1F *bgNC = new TH1F("bgNC", "NN output", 50, -.5, 1.5);
477  TH1F *bgNueCC= new TH1F("bgNueCC", "NN output", 50, -.5, 1.5);
478 
479  TH1F *sig = new TH1F("sigh", "NN output", 50, -.5, 1.5);
480  TH1F *sigbg = new TH1F("sigbgh", "NN output", 50, -.5, 1.5);
481 
482 
483 // bg->SetDirectory(0);
484 // bgNumuCC->SetDirectory(0);
485 // bgNC->SetDirectory(0);
486 // bgNueCC->SetDirectory(0);
487 // sig->SetDirectory(0);
488 
489  Double_t params[18];
490  for (int i = 0; i < background.GetEntries(); i++) {
491 // if((double)i>(double)background.GetEntries()*rbg)continue;
492  if(i%2==0)continue;
493 // if(i>36383)continue;
494  background.GetEntry(i);
495  if(!((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))continue;
496  if(!(evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
497  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
498 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
499 // if(evtSh1DedxLLL[0]<0)continue;
500 // if(!((evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
501  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
502 // if(!(evtSh1Energy>Elow&&evtSh1Energy<Ehigh))continue;
503 // if(!(evtSh1Enue/1.40694>1.0&&evtSh1Enue/1.40694<Ehigh))continue;
504 
505 
506  params[0] = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
507  params[1] = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
508  params[2] = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
509  params[3] = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
510  params[4] = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
511  params[5] = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
512  params[6] = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
513  params[7] = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
514  params[8] = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
515  params[9] = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
516  params[10] = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
517  params[11] = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
518  params[12] = MAX(evtSh1Pi0Mgg,0);
519  params[13] = evtSh1Energy/evtEtot;
520  params[14] = evtSh1VtxGeV;
521  params[15] = evtSh1Gap;
522  params[17] = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
523  params[16] = evtSh1CosTheta;
524 
525 // params[17] = evtSh1SumPt/evtSh1SliceGeV;
526 // params[18] = evtSh1Length/evtSh1Energy;
527 
528 
529 
530 // params[8] = sqrt((evtSh1Start[0]-evtTrueNuVtx[0])**2+(evtSh1Start[1]-evtTrueNuVtx[1])**2+(evtSh1Start[2]-evtTrueNuVtx[2])**2);
531 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[2])>10)continue;
532 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[2])>10)continue;
533 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[3])>10)continue;
534 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[3])>10)continue;
535 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[5])>10)continue;
536 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[5])>10)continue;
537 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[6])>10)continue;
538 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[6])>10)continue;
539 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[7])>10)continue;
540 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[7])>10)continue;
541 // bg->Fill(mlp->Evaluate(0, params));
542  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2]){
543  bgNumuCC->Fill(mlp->Evaluate(0, params), scalepotNumuCC);
544  bg->Fill(mlp->Evaluate(0, params), scalepotNumuCC);
545  }
546  if(evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1]))){
547  bgNC->Fill(mlp->Evaluate(0, params), scalepotNC);
548  bg->Fill(mlp->Evaluate(0, params), scalepotNC);
549  }
550  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1]){
551  bgNueCC->Fill(mlp->Evaluate(0, params), scalepotNueCC);
552  bg->Fill(mlp->Evaluate(0, params), scalepotNueCC);
553  }
554  }
555 
556  for (int i = 0; i < signal.GetEntries(); i++) {
557 // if((double)i>(double)signal.GetEntries()*rsig)continue;
558  if(i%2==0)continue;
559  signal.GetEntry(i);
560  if(!(evtTrueNuCCNC==0&&evtTrueNuPdg==12&&(evtRwRand[0]>0&&evtRwRand[0]<evtRwProb[0])))continue;
561 // if(evtSh1TrueDang>10||fabs(evtSh1TruePdg)!=11)continue;
562  if(!(evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
563 
564 
565  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
566 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
567 // if(evtSh1DedxLLL[0]<0)continue;
568 // if(!((evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
569  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
570 
571 // if(!(evtSh1Energy>Elow&&evtSh1Energy<Ehigh))continue;
572 // if(!(evtSh1Enue/1.40694>1.0&&evtSh1Enue/1.40694<Ehigh))continue;
573 
574 
575  params[0] = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
576  params[1] = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
577  params[2] = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
578  params[3] = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
579  params[4] = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
580  params[5] = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
581  params[6] = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
582  params[7] = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
583  params[8] = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
584  params[9] = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
585  params[10] = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
586  params[11] = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
587  params[12] = MAX(evtSh1Pi0Mgg,0);
588  params[13] = evtSh1Energy/evtEtot;
589  params[14] = evtSh1VtxGeV;
590  params[15] = evtSh1Gap;
591  params[17] = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
592  params[16] = evtSh1CosTheta;
593 
594 // params[17] = evtSh1SumPt/evtSh1SliceGeV;
595 // params[18] = evtSh1Length/evtSh1Energy;
596 
597 
598 // params[8] = sqrt((evtSh1Start[0]-evtTrueNuVtx[0])**2+(evtSh1Start[1]-evtTrueNuVtx[1])**2+(evtSh1Start[2]-evtTrueNuVtx[2])**2);
599 
600 // if(fabs( evtSh1DedxLLL[0] - evtSh1DedxLLL[2])>10)continue;
601 // if(fabs( evtSh1DedxLLT[0] - evtSh1DedxLLT[2])>10)continue;
602 // if(fabs( evtSh1DedxLLL[0] - evtSh1DedxLLL[3])>10)continue;
603 // if(fabs( evtSh1DedxLLT[0] - evtSh1DedxLLT[3])>10)continue;
604 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[5])>10)continue;
605 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[5])>10)continue;
606 // if(fabs(evtSh1DedxLLL[0] - evtSh1DedxLLL[6])>10)continue;
607 // if(fabs(evtSh1DedxLLT[0] - evtSh1DedxLLT[6])>10)continue;
608 // if(fabs( evtSh1DedxLLL[0] - evtSh1DedxLLL[7])>10)continue;
609 // if(fabs( evtSh1DedxLLT[0] - evtSh1DedxLLT[7])>10)continue;
610 // cout<<"sig eval "<<mlp->Evaluate(0,params)<<endl;
611  sig->Fill(mlp->Evaluate(0,params),scalepotsig);
612  if(fabs(evtSh1TruePdg)!=11)sigbg->Fill(mlp->Evaluate(0,params),scalepotsig);
613  }
614 
615 // TH1F *bg = new TH1F("bgh", "NN output", 50, -.5, 1.5);
616 // TH1F *sig = new TH1F("sigh", "NN output", 50, -.5, 1.5);
617  TH1F *hfom = new TH1F("hfom", "NN output", 50, -.5, 1.5);
618  TH1F *heff = new TH1F("heff", "NN output", 50, -.5, 1.5);
619  TH1F *hbgeff = new TH1F("hbgeff", "NN output", 50, -.5, 1.5);
620  TH1F *hbgNumuCCeff = new TH1F("hbgNumuCCeff", "NN output", 50, -.5, 1.5);
621  TH1F *hbgNCeff = new TH1F("hbgNCeff", "NN output", 50, -.5, 1.5);
622  TH1F *hbgNueCCeff = new TH1F("hbgNueCCeff", "NN output", 50, -.5, 1.5);
623  for(int i=0;i<50;i++){
624  double sumnbg = bg->Integral(i+1, 50);
625  double sumnbgNueCC = bgNueCC->Integral(i+1, 50);
626  double sumnsig = sig->Integral(i+1, 50);
627  double eff = double(sig->Integral(i+1, 50))/double(sig->Integral(1, 50));
628  double bgeff = double(bg->Integral(i+1, 50))/double(bg->Integral(1, 50));
629  double bgNumuCCeff = double(bgNumuCC->Integral(i+1, 50))/double(bgNumuCC->Integral(1, 50));
630  double bgNCeff = double(bgNC->Integral(i+1, 50))/double(bgNC->Integral(1, 50));
631  double bgNueCCeff = double(bgNueCC->Integral(i+1, 50))/double(bgNueCC->Integral(1, 50));
632  heff->SetBinContent(i+1,eff);
633  hbgeff->SetBinContent(i+1,bgeff);
634  hbgNumuCCeff->SetBinContent(i+1,bgNumuCCeff);
635  hbgNCeff->SetBinContent(i+1,bgNCeff);
636  hbgNueCCeff->SetBinContent(i+1,bgNueCCeff);
637  if(sumnsig+sumnbg-sumnbgNueCC<0.000001){
638  hfom->SetBinContent(i+1,0);
639  continue;
640  }
641  double significance = 0;
642  if(sumnsig+sumnbg-sumnbgNueCC>0) significance = double(sumnsig)/sqrt(sumnsig+sumnbg-sumnbgNueCC);
643 // if(sumnbg>0)significance = double(sumnsig)/sqrt(sumnbg);
644  cout<<"FOM, sig, bkg "<<significance<<" "<<sumnsig<<" "<<sumnbg<<" "<<endl;
645  hfom->SetBinContent(i+1,significance);
646  }
647 
648  bg->SetLineColor(kBlue);
649 // bg->SetFillStyle(3008);
650 // bg->SetFillColor(kBlue);
651  sig->SetLineColor(kRed);
652  sig->SetFillStyle(3003);
653  sig->SetFillColor(kRed);
654 // bg->SetStats(0);
655 // sig->SetStats(0);
656  bg->Draw();
657  sig->Draw("same");
658  sigbg->Draw("same");
659 // sig->DrawNormalized("same", sig->GetEntries()*scalesig);
660  TLegend *legend = new TLegend(.75, .80, .95, .95);
661  legend->AddEntry(bg, "Background (nhc)");
662  legend->AddEntry(sig, "Signal (Nue_CCQE)");
663  legend->Draw();
664  mlpa_canvas2->cd(0);
665  mlpa_canvas2->Update();
666 
667  cout<<"scalepotsig "<<scalepotsig<<endl;
668  cout<<"scalepotNumuCC "<<scalepotNumuCC<<endl;
669  cout<<"scalepotNC "<<scalepotNC<<endl;
670  cout<<"scalepotNueCC "<<scalepotNueCC<<endl;
671 
672  TCanvas* mlpa_canvas3 = new TCanvas("mlpa_canvas3","Network analysis 3", 10,10, 800,400);
673  mlpa_canvas3->Divide(2,1);
674  mlpa_canvas3->cd(1);
675  hfom->GetYaxis()->SetTitle("FOM");
676  hfom->GetXaxis()->SetTitle("ANN");
677  hfom->GetXaxis()->SetRangeUser(0.4,1.2);
678  hfom->Draw("l");
679  line1 = new TLine(hfom->GetBinCenter(hfom->GetMaximumBin()),0,hfom->GetBinCenter(hfom->GetMaximumBin()),hfom->GetMaximum());
680  line1->SetLineWidth(4);
681  line1->SetLineColor(kRed);
682  line1->Draw("same");
683 
684  mlpa_canvas3->Update();
685 
686  mlpa_canvas3->cd(2);
687  heff->GetYaxis()->SetTitle("PID Eff.");
688  heff->GetXaxis()->SetTitle("ANN");
689  heff->GetYaxis()->SetRangeUser(0,1);
690  heff->GetXaxis()->SetRangeUser(0.4,1.2);
691  heff->Draw("l");
692  line2 = new TLine(hfom->GetBinCenter(hfom->GetMaximumBin()),0,hfom->GetBinCenter(hfom->GetMaximumBin()),1.0);
693  line2->SetLineWidth(4);
694  line2->SetLineColor(kRed);
695  line2->Draw("same");
696 
697  mlpa_canvas3->Update();
698 
699  TCanvas* mlpa_canvas4 = new TCanvas("mlpa_canvas4","Network analysis 4", 10,10, 600,400);
700  mlpa_canvas4->Divide(1,1);
701  mlpa_canvas4->cd(1);
702  bg->GetYaxis()->SetTitle("Events");
703  bg->GetXaxis()->SetTitle("ANN");
704  bg->SetLineColor(kBlue);
705  bgNumuCC->SetLineColor(9);
706  bgNC->SetLineColor(7);
707  bgNueCC->SetLineColor(6);
708  sig->SetLineColor(kRed);
709  sig->SetFillStyle(3003);
710  sig->SetFillColor(kRed);
711 // bg->SetStats(0);
712 // sig->SetStats(0);
713 
714  bg->Draw();
715  bgNumuCC->Draw("same");
716  bgNC->Draw("same");
717  bgNueCC->Draw("same");
718  sig->Draw("same");
719 // sig->DrawNormalized("same", sig->GetEntries()*scalesig);
720  TLegend *legend = new TLegend(.75, .80, .95, .95);
721  legend->AddEntry(bg, "Background (nhc)");
722  legend->AddEntry(bgNumuCC, "Beam #nu_{#mu} CC");
723  legend->AddEntry(bgNC, "Beam NC");
724  legend->AddEntry(bgNueCC, "Beam #nu_{e} CC");
725  legend->AddEntry(sig, "Signal #nu_{e} CC");
726  legend->Draw();
727  mlpa_canvas4->Update();
728 
729 
730 
731  type = 1;
732  for (int i = 0; i < signal.GetEntries(); i++) {
733 // if((double)i>(double)signal.GetEntries()*rsig)continue;
734  if(i%2==0)continue;
735  signal.GetEntry(i);
736  if(!(evtTrueNuCCNC==0&&evtTrueNuPdg==12&&(evtRwRand[0]>0&&evtRwRand[0]<evtRwProb[0])))continue;
737 // if(evtSh1TrueDang>10||fabs(evtSh1TruePdg)!=11)continue;
738 // if(!(evtSh1IsFid==1&&evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
739  if(!(evtSh1Energy>0.0))continue;
740  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
741 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
742 // if(!((evtSh1Energy+0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
743  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
744 // if(evtSh1DedxLLL[0]<0)continue;
745  egLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
746  egLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
747  emuLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
748  emuLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
749  epi0LLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
750  epi0LLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
751  epLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
752  epLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
753  enLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
754  enLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
755  epiLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
756  epiLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
757  gap = evtSh1Gap;
758 // costheta = evtSh1SumPt/evtSh1SliceGeV;
759  costheta = evtSh1CosTheta;
760  length = evtSh1Length/evtSh1Energy;
761  shE = evtSh1Energy/evtSh1SliceGeV;
762  nuE = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
763  ntrk = evtNTrk;
764  pi0mass = MAX(evtSh1Pi0Mgg,0);
765  annPi0LineDist = evtPi0LineDist;
766  annPi0StartDist = evtPi0StartDist;
767 
768  params[0] = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
769  params[1] = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
770  params[2] = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
771  params[3] = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
772  params[4] = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
773  params[5] = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
774  params[6] = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
775  params[7] = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
776  params[8] = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
777  params[9] = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
778  params[10] = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
779  params[11] = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
780  params[12] = MAX(evtSh1Pi0Mgg,0);
781  params[13] = evtSh1Energy/evtEtot;
782  params[14] = evtSh1VtxGeV;
783  params[15] = evtSh1Gap;
784  params[17] = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
785  params[16] = evtSh1CosTheta;
786 
787 // params[17] = evtSh1SumPt/evtSh1SliceGeV;
788 // params[18] = evtSh1Length/evtSh1Energy;
789 
790 
791 // params[8] = sqrt((evtSh1Start[0]-evtTrueNuVtx[0])**2+(evtSh1Start[1]-evtTrueNuVtx[1])**2+(evtSh1Start[2]-evtTrueNuVtx[2])**2);
792  trueNuE = evtTrueNuEnergy;
793  etot = evtEtot;
794  enue = evtSh1Enue;
795  sh1Energy = evtSh1Energy;
796  ANN = mlp->Evaluate(0,params);
797  wSig = scalepotsig;
798  wBkg = 0;
799  annMC->Fill();
800  if(ANN<hfom->GetBinCenter(hfom->GetMaximumBin()))continue;
801 // if(!((evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
802  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
803 // if(!(evtSh1Energy>Elow&&evtSh1Energy<Ehigh))continue;
804 // if(!(evtSh1Enue/1.40694>1.0&&evtSh1Enue/1.40694<Ehigh))continue;
805 
806  ns3++;
807  }
808 
809  type = 0;
810  for (int i = 0; i < background.GetEntries(); i++) {
811 // if((double)i>(double)background.GetEntries()*rbg)continue;
812  if(i%2==0)continue;
813 // if(i>36383)continue;
814  background.GetEntry(i);
815  if(!((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))continue;
816 // if(!(evtSh1IsFid==1&&evtSh1NPlane>-15&&evtSh1Energy>0.0))continue;
817  if(!(evtSh1Energy>0.0))continue;
818  if(evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
819 // if(evtPi0Mgg>0.05&&evtPi0Mgg<0.22)continue;
820 // if(!((evtSh1Energy+0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
821  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
822 // if(evtSh1DedxLLL[0]<0)continue;
823  egLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
824  egLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
825  emuLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
826  emuLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
827  epi0LLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
828  epi0LLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
829  epLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
830  epLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
831  enLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
832  enLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
833  epiLLL = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
834  epiLLT = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
835  gap = evtSh1Gap;
836 // costheta = evtSh1SumPt/evtSh1SliceGeV;
837  costheta = evtSh1CosTheta;
838  length = evtSh1Length/evtSh1Energy;
839  shE = evtSh1Energy/evtEtot;
840  nuE = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
841  ntrk = evtNTrk;
842  pi0mass = MAX(evtSh1Pi0Mgg,0);
843  annPi0LineDist = evtPi0LineDist;
844  annPi0StartDist = evtPi0StartDist;
845 
846 
847 
848  params[0] = evtSh1DedxLLL[0] - evtSh1DedxLLL[1];
849  params[1] = evtSh1DedxLLT[0] - evtSh1DedxLLT[1];
850  params[2] = evtSh1DedxLLL[0] - evtSh1DedxLLL[2];
851  params[3] = evtSh1DedxLLT[0] - evtSh1DedxLLT[2];
852  params[4] = evtSh1DedxLLL[0] - evtSh1DedxLLL[3];
853  params[5] = evtSh1DedxLLT[0] - evtSh1DedxLLT[3];
854  params[6] = evtSh1DedxLLL[0] - evtSh1DedxLLL[5];
855  params[7] = evtSh1DedxLLT[0] - evtSh1DedxLLT[5];
856  params[8] = evtSh1DedxLLL[0] - evtSh1DedxLLL[6];
857  params[9] = evtSh1DedxLLT[0] - evtSh1DedxLLT[6];
858  params[10] = evtSh1DedxLLL[0] - evtSh1DedxLLL[7];
859  params[11] = evtSh1DedxLLT[0] - evtSh1DedxLLT[7];
860  params[12] = MAX(evtSh1Pi0Mgg,0);
861  params[13] = evtSh1Energy/evtEtot;
862  params[14] = evtSh1VtxGeV;
863  params[15] = evtSh1Gap;
864  params[17] = evtSh1NueEnergy;//(evtSh1Energy + 0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015;
865  params[16] = evtSh1CosTheta;
866 
867 // params[17] = evtSh1SumPt/evtSh1SliceGeV;
868 // params[18] = evtSh1Length/evtSh1Energy;
869 
870 
871 // params[8] = sqrt((evtSh1Start[0]-evtTrueNuVtx[0])**2+(evtSh1Start[1]-evtTrueNuVtx[1])**2+(evtSh1Start[2]-evtTrueNuVtx[2])**2);
872  trueNuE = evtTrueNuEnergy;
873  etot = evtEtot;
874  enue = evtSh1Enue;
875  sh1Energy = evtSh1Energy;
876  ANN = mlp->Evaluate(0,params);
877  wSig = 0;
878  wBkg = 0;
879  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])wBkg=scalepotNumuCC;
880  if(evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))wBkg=scalepotNC;
881  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])wBkg=scalepotNueCC;
882  annMC->Fill();
883 // if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12)nbNueCC2_a++;
884  if(ANN<hfom->GetBinCenter(hfom->GetMaximumBin()))continue;
885 // if(!((evtSh1Energy+0.282525 + 1.0766*(evtEtot-evtSh1Energy))/1.015>Elow&&(evtSh1Energy+0.282525+1.0766*(evtEtot-evtSh1Energy))/1.015<Ehigh))continue;
886  if(!(evtSh1NueEnergy>Elow&&evtSh1NueEnergy<Ehigh))continue;
887 // if(!(evtSh1Energy>Elow&&evtSh1Energy<Ehigh))continue;
888 // if(!(evtSh1Enue/1.40694>1.0&&evtSh1Enue/1.40694<Ehigh))continue;
889  nb3++;
890  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])nbNumuCC3++;
891  if(evtTrueNuCCNC==1&&((fabs(evtTrueNuPdg)==14&&evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])||(fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])))nbNC3++;
892  if(evtTrueNuCCNC==0&&fabs(evtTrueNuPdg)==12&&evtRwRand[1]>0&&evtRwRand[1]<evtRwProb[1])nbNueCC3++;
893  }
894  cout<<"nb3 "<<nb3<<endl;
895  cout<<"nbNumuCC3 "<<nbNumuCC3<<endl;
896  cout<<"nbNC3 "<<nbNC3<<endl;
897  cout<<"nbNueCC3 "<<nbNueCC3<<endl;
898 
899  TFile f("ann.root","recreate");
900  annMC->Write("annMC");
901  f.Close();
902 
903  std::cout<<"ns0,ns1,ns2,ns3,scalepotsig "<<ns0<<" "<<ns1<<" "<<ns2<<" "<<ns3<<" "<<scalepotsig<<endl;
904  double Dns0 = double(ns0*scalepotsig);
905  double Dns1 = double(ns1*scalepotsig);
906  double Dns2 = double(ns2*scalepotsig);
907  double Dns3 = double(ns3*scalepotsig);
908 
909  double DnbNumuCC0 = double(nbNumuCC0*scalepotNumuCC);
910  double DnbNumuCC1 = double(nbNumuCC1*scalepotNumuCC);
911  double DnbNumuCC2 = double(nbNumuCC2*scalepotNumuCC);
912  double DnbNumuCC3 = double(nbNumuCC3*scalepotNumuCC);
913 
914  double DnbNC0 = double(nbNC0*scalepotNC);
915  double DnbNC1 = double(nbNC1*scalepotNC);
916  double DnbNC2 = double(nbNC2*scalepotNC);
917  double DnbNC3 = double(nbNC3*scalepotNC);
918 
919  double DnbNueCC0 = double(nbNueCC0*scalepotNueCC);
920  double DnbNueCC1 = double(nbNueCC1*scalepotNueCC);
921  double DnbNueCC2 = double(nbNueCC2*scalepotNueCC);
922  double DnbNueCC3 = double(nbNueCC3*scalepotNueCC);
923 
924  double Dnb0 = double(DnbNumuCC0+DnbNC0+DnbNueCC0);
925  double Dnb1 = double(DnbNumuCC1+DnbNC1+DnbNueCC1);
926  double Dnb2 = double(DnbNumuCC2+DnbNC2+DnbNueCC2);
927  double Dnb3 = double(DnbNumuCC3+DnbNC3+DnbNueCC3);
928 
929 // double scalepot = 18.0e+20/((nbNumuCC0/7175.2)*2.8123e+21);//2.38399959425089536e+19);//NDOS
930 
931 
932  cout<<"Scaled sig "<<Dns0<<" "<<Dns1<<" "<<Dns2<<" "<<Dns3<<endl;
933  cout<<"Scaled bg "<<Dnb0<<" "<<Dnb1<<" "<<Dnb2<<" "<<Dnb3<<endl;
934  cout<<"Scaled bg NumuCC"<<DnbNumuCC0<<" "<<DnbNumuCC1<<" "<<DnbNumuCC2<<" "<<DnbNumuCC3<<endl;
935  cout<<"Scaled bg NC "<<DnbNC0<<" "<<DnbNC1<<" "<<DnbNC2<<" "<<DnbNC3<<endl;
936  cout<<"Scaled bg NueCC "<<DnbNueCC0<<" "<<DnbNueCC1<<" "<<DnbNueCC2<<" "<<DnbNueCC3<<endl;
937 
938  cout<<"ANN = "<<hfom->GetBinCenter(hfom->GetMaximumBin())<<endl;
939  cout<<"FOM hist = "<<hfom->GetBinContent(hfom->GetMaximumBin())<<endl;
940  cout<<"FOM = "<<double(Dns3)/sqrt(Dns3+Dnb3)<<endl;
941  cout<<"Eff. = "<<double(Dns3)/double(Dns0)<<endl;
942  cout<<"ns2 "<<ns2<<endl;
943  cout<<"nb2, nb3 "<<nb2<<" "<<nb3<<endl;
944 
945 /*
946  cout<<"BG Eff . = "<<hbgeff->GetBinContent(hfom->GetMaximumBin())<<endl;
947  cout<<"BG NumuCC Eff. = "<<hbgNumuCCeff->GetBinContent(hfom->GetMaximumBin())<<endl;
948  cout<<"BG NC Eff. = "<<hbgNCeff->GetBinContent(hfom->GetMaximumBin())<<endl;
949  cout<<"BG NueCC Eff. = "<<hbgNueCCeff->GetBinContent(hfom->GetMaximumBin())<<endl;
950 */
951  cout<<"BG Eff . = "<<double(Dnb3)/double(Dnb0)<<endl;
952  cout<<"BG NumuCC Eff. = "<<double(DnbNumuCC3)/double(DnbNumuCC0)<<endl;
953  cout<<"BG NC Eff. = "<<double(DnbNC3)/double(DnbNC0)<<endl;
954  cout<<"BG NueCC Eff. = "<<double(DnbNueCC3)/double(DnbNueCC0)<<endl;
955 }
956 double MAX(double a, double b){
957  if(a>b){
958  return a;
959  }
960  else{
961  return b;
962  }
963 }
void mlpPIDelecallEetot(Int_t ntrain=600)
TString fin
Definition: Style.C:24
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
length
Definition: demo0.py:21
const double a
double MAX(double a, double b)
OStream cout
Definition: OStream.cxx:6
const hit & b
Definition: hits.cxx:21
Float_t e
Definition: plot.C:35
enum BeamMode kBlue