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