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