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