00001 #ifndef madedanalysis_cxx
00002 #define madedanalysis_cxx
00003 #include <iostream>
00004 #include <fstream>
00005 #include <cmath>
00006 #include "TVector3.h"
00007 #include "TEventList.h"
00008 #include "TProfile.h"
00009 #include "TTree.h"
00010 #include "Mad/MadEdAnalysis.h"
00011
00012 using namespace std;
00013
00014
00015 MadEdAnalysis::MadEdAnalysis(TChain *chainSR,TChain *chainMC,
00016 TChain *chainTH,TChain *chainEM)
00017 {
00018
00019 if(!chainSR) {
00020 record = 0;
00021 strecord = 0;
00022 emrecord = 0;
00023 mcrecord = 0;
00024 threcord = 0;
00025 Clear();
00026 whichSource = -1;
00027 isMC = true;
00028 isTH = true;
00029 isEM = true;
00030 Nentries = -1;
00031 return;
00032 }
00033
00034 InitChain(chainSR,chainMC,chainTH,chainEM);
00035 whichSource = 0;
00036
00037 }
00038
00039
00040 MadEdAnalysis::MadEdAnalysis(JobC *j,string path,int entries)
00041 {
00042 record = 0;
00043 strecord = 0;
00044 emrecord = 0;
00045 mcrecord = 0;
00046 threcord = 0;
00047 Clear();
00048 isMC = true;
00049 isTH = true;
00050 isEM = true;
00051 Nentries = entries;
00052 whichSource = 1;
00053 jcPath = path;
00054 fJC = j;
00055 fChain = NULL;
00056
00057 }
00058
00059
00060 MadEdAnalysis::~MadEdAnalysis()
00061 {
00062 }
00063
00064
00065
00066
00067 Bool_t MadEdAnalysis::MyIsFidVtx(Int_t itrk){
00068 if(!LoadTrack(itrk)) return false;
00069
00070
00071
00072
00073
00074 Float_t trkendx=0;
00075 Float_t trkendy=0;
00076 Float_t trkendz=0;
00077
00078 trkendx=ntpTrack->end.x;
00079 trkendy=ntpTrack->end.y;
00080 trkendz=ntpTrack->end.z;
00081
00082
00083 if (!(trkendz<16 && sqrt(pow(trkendx,2)+pow(trkendy,2))>0.4 &&
00084 trkendx<2.7 && trkendx>-1.65 && trkendy<1.65 && trkendy>-1.65 &&
00085 trkendy>(-trkendx)-1.65 && trkendy<trkendx+1.65 &&
00086 trkendy<(-trkendx)+3.55 && trkendy>trkendx-3.55) ) return false;
00087
00088
00089 return true;
00090 }
00091
00092
00093
00094 Bool_t MadEdAnalysis::MyIsFidVtxrz(Int_t itrk){
00095 if(!LoadTrack(itrk)) return false;
00096
00097 if(ntpTrack->fidvtx.dr<0.5||ntpTrack->fidvtx.dz<0.5) return false;
00098 return true;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 Bool_t MadEdAnalysis::PassAnalysisCuts(Int_t event)
00110 {
00111 if(!LoadEvent(event)) return false;
00112
00113 Float_t pidcut=-0.4;
00114 if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
00115
00116 if(PID(event,0)>pidcut) return true;
00117 return false;
00118 }
00119
00120
00121 Bool_t MadEdAnalysis::PassBasicCuts(Int_t event)
00122 {
00123
00124 if(!LoadEvent(event)) return false;
00125 if(ntpEvent->ntrack==0) return false;
00126 Int_t track = -1;
00127 LoadLargestTrackFromEvent(event,track);
00128 if(!PassTrackCuts(track)) return false;
00129
00130 return true;
00131 }
00132
00133
00134
00135
00136
00137
00138 Float_t MadEdAnalysis::RecoMuDCosY(Int_t itrk){
00139 if(!LoadTrack(itrk)) return 0.;
00140 return ntpTrack->vtx.dcosy;
00141 }
00142
00143
00144
00145
00146 Float_t MadEdAnalysis::RecoMuZn(Int_t itrk){
00147 if(!LoadTrack(itrk)) return 0.;
00148 return ntpTrack->cr.zenith;
00149 }
00150
00151
00152
00153
00154
00155
00156 Float_t MadEdAnalysis::RecoMuAZM(Int_t itrk){
00157 if(!LoadTrack(itrk)) return 0.;
00158 return ntpTrack->cr.azimuth;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168 Float_t MadEdAnalysis::PID(Int_t event,Int_t method){
00169
00170
00171 if(method==0) {
00172
00173 if(!LoadEvent(event)) return -10;
00174 Int_t track = -1;
00175 if(!LoadLargestTrackFromEvent(event,track)) return -10;
00176 if(method!=0) return -10;
00177
00178 Float_t ProbNC = 1.;
00179 Float_t ProbCC = 1.;
00180 Float_t PidCC = 0.;
00181
00182 if(!MyIsFidVtx()) return false;
00183
00184 Float_t var1=eventSummary->plane.n;
00185 Float_t var2=0;
00186 Float_t var3=0;
00187
00188
00189
00190 if (eventSummary->ph.sigcor!=0) var2 = float(ntpTrack->ph.sigcor)
00191 /float(eventSummary->ph.sigcor);
00192
00193
00194 if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor)
00195 /float(ntpTrack->plane.n);
00196
00197
00198
00199
00200
00201 Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00202 Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00203 Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00204
00205 if (prob1==0) {prob1=0.0001;}
00206 if (prob2==0) {prob2=0.0001;}
00207 if (prob3==0) {prob3=0.0001;}
00208
00209 ProbCC=prob1*prob2*prob3;
00210
00211
00212 prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00213 prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00214 prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00215
00216 if (prob1==0) {prob1=0.0001;}
00217 if (prob2==0) {prob2=0.0001;}
00218 if (prob3==0) {prob3=0.0001;}
00219
00220 ProbNC=prob1*prob2*prob3;
00221
00222
00223 PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00224 return PidCC;
00225
00226
00227
00228 }
00229
00230
00231
00232 if(method==1){
00233 if(!LoadEvent(event)) return -10;
00234 return NeuNetEval(event);
00235 }
00236
00237 return 0;
00238
00239
00240 }
00241
00242
00243
00244
00245
00246 void MadEdAnalysis::MyCreatePAN(std::string tag){
00247
00248 std::string savename = "PAN_" + tag + ".root";
00249 TFile save(savename.c_str(),"RECREATE");
00250 save.cd();
00251
00252
00253
00254 ofstream ofs("text.txt");
00255
00256
00257 GnumiInterface gnumi;
00258 if(!gnumi.Status()) {
00259 cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00260 << " Will not be filling NuParent info"
00261 << endl;
00262 cout << "Set environmental variable $GNUMIAUX to point to the "
00263 << "directory containing the gnumi files to fix this!"
00264 << endl;
00265
00266 }
00267
00268
00269
00270 Float_t true_enu;
00271 Float_t true_pxnu;
00272 Float_t true_pynu;
00273 Float_t true_pznu;
00274 Float_t true_etgt;
00275 Float_t true_pxtgt;
00276 Float_t true_pytgt;
00277 Float_t true_pztgt;
00278 Float_t true_emu;
00279 Float_t true_eshw;
00280 Float_t true_x;
00281 Float_t true_y;
00282 Float_t true_q2;
00283 Float_t true_w2;
00284 Float_t true_dircosneu;
00285 Float_t true_dircosz;
00286 Float_t true_vtxx;
00287 Float_t true_vtxy;
00288 Float_t true_vtxz;
00289
00290
00291 Int_t flavour;
00292 Int_t process;
00293 Int_t nucleus;
00294 Int_t cc_nc;
00295 Int_t initial_state;
00296 Int_t had_fs;
00297
00298
00299 NuParent *nuparent = new NuParent();
00300
00301
00302 Int_t detector;
00303 Int_t run;
00304 Int_t snarl;
00305 Int_t nevent;
00306 Int_t event;
00307 Int_t mcevent;
00308 Int_t ntrack;
00309 Int_t nshower;
00310
00311 Int_t is_fid;
00312 Int_t is_fidrk;
00313
00314
00315 Int_t is_pcev;
00316 Int_t is_cev;
00317 Int_t is_mc;
00318 Int_t pass_basic;
00319 Float_t pid0;
00320 Float_t pid1;
00321 Float_t pid2;
00322 Float_t Pidvalsk;
00323
00324 Int_t pass;
00325
00326 Float_t reco_enu;
00327 Float_t reco_emu;
00328 Float_t reco_eshw;
00329 Float_t reco_eshw_sqrt;
00330 Float_t reco_qe_enu;
00331 Float_t reco_qe_q2;
00332 Float_t reco_y;
00333 Float_t reco_dircosneu;
00334 Float_t reco_dircosz;
00335 Float_t reco_dircosy;
00336 Float_t reco_azimuth;
00337 Float_t reco_vtxx;
00338 Float_t reco_vtxy;
00339 Float_t reco_vtxz;
00340
00341 Float_t evtlength;
00342 Float_t trklength;
00343 Float_t trkmomrange;
00344 Float_t trkqp;
00345 Float_t trkeqp;
00346 Float_t trkphfrac;
00347 Float_t trkphplane;
00348
00349 Float_t trkvtxx=0;
00350 Float_t trkvtxy=0;
00351 Float_t trkvtxz=0;
00352
00353 Float_t trkendx=0;
00354 Float_t trkendy=0;
00355 Float_t trkendz=0;
00356
00357
00358
00359
00360
00361 TTree *tree = new TTree("pan","pan");
00362
00363 tree->Branch("true_enu",&true_enu,"true_enu/F",32000);
00364 tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",32000);
00365 tree->Branch("true_pynu",&true_pynu,"true_pynu/F",32000);
00366 tree->Branch("true_pznu",&true_pznu,"true_pznu/F",32000);
00367 tree->Branch("true_etgt",&true_etgt,"true_etgt/F",32000);
00368 tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",32000);
00369 tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",32000);
00370 tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",32000);
00371 tree->Branch("true_emu",&true_emu,"true_emu/F",32000);
00372 tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000);
00373 tree->Branch("true_x",&true_x,"true_x/F",32000);
00374 tree->Branch("true_y",&true_y,"true_y/F",32000);
00375 tree->Branch("true_q2",&true_q2,"true_q2/F",32000);
00376 tree->Branch("true_w2",&true_w2,"true_w2/F",32000);
00377 tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000);
00378 tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000);
00379 tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000);
00380 tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000);
00381 tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000);
00382
00383 tree->Branch("flavour",&flavour,"flavour/F",32000);
00384 tree->Branch("process",&process,"process/I",32000);
00385 tree->Branch("nucleus",&nucleus,"nucleus/I",32000);
00386 tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000);
00387 tree->Branch("initial_state",&initial_state,"initial_state/I",32000);
00388 tree->Branch("had_fs",&had_fs,"had_fs/I",32000);
00389
00390 tree->Branch("nuparent","NuParent",&nuparent,8000,1);
00391
00392 tree->Branch("detector",&detector,"detector/I",32000);
00393 tree->Branch("run",&run,"run/I",32000);
00394 tree->Branch("snarl",&snarl,"snarl/I",32000);
00395 tree->Branch("event",&event,"event/I",32000);
00396 tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
00397 tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
00398 tree->Branch("nshower",&nshower,"nshower/I",32000);
00399
00400 tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
00401 tree->Branch("is_fidrk",&is_fidrk,"is_fidrk/I",32000);
00402 tree->Branch("is_pcev",&is_pcev,"is_pcev/I",32000);
00403
00404
00405 tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
00406 tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
00407 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
00408 tree->Branch("Pidvalsk",&Pidvalsk,"Pidvalsk/F",32000);
00409 tree->Branch("pid0",&pid0,"pid0/F",32000);
00410 tree->Branch("pass",&pass,"pass/I",32000);
00411
00412 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
00413 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
00414 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
00415 tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",32000);
00416 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
00417 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);
00418 tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
00419 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
00420 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);
00421 tree->Branch("reco_dircosy",&reco_dircosy,"reco_dircosy/F",32000);
00422 tree->Branch("reco_azimuth",&reco_azimuth,"reco_azimuth/F",32000);
00423
00424 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
00425 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
00426 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
00427
00428 tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
00429 tree->Branch("trklength",&trklength,"trklength/F",32000);
00430 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
00431 tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
00432 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
00433 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
00434 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
00435
00436 tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",32000);
00437 tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",32000);
00438 tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",32000);
00439
00440 tree->Branch("trkendx",&trkendx,"trkendx/F",32000);
00441 tree->Branch("trkendy",&trkendy,"trkendy/F",32000);
00442 tree->Branch("trkendz",&trkendz,"trkendz/F",32000);
00443
00444
00445
00446
00447
00448 for(int i=0;i<Nentries;i++){
00449 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
00450 << "% done" << std::endl;
00451
00452 if(!this->GetEntry(i)) continue;
00453
00454 nevent = eventSummary->nevent;
00455 if(nevent==0) continue;
00456
00457 run = ntpHeader->GetRun();
00458 snarl = ntpHeader->GetSnarl();
00459
00460 Int_t trkcount=0;
00461 Int_t shwcount=0;
00462 Float_t totph=0;
00463
00464
00465 for(int i=0;i<nevent;i++){
00466 if(!LoadEvent(i)) continue;
00467
00468 totph = 0;
00469
00470
00471 event = i;
00472 ntrack = ntpEvent->ntrack;
00473 nshower = ntpEvent->nshower;
00474
00475
00476 true_enu = 0; true_emu = 0; true_eshw = 0;
00477 true_pxnu = 0; true_pynu = 0; true_pznu = 0;
00478 true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
00479 true_dircosneu = 0; true_dircosz = 0;
00480 true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
00481 flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
00482 initial_state = 0; had_fs = 0;
00483 true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;
00484 nuparent->Zero();
00485 detector = -1; mcevent = -1;
00486 is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0;
00487 pid0 = 0; pid1 = 0; pid2 = 0; pass = 0;
00488 reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0;
00489 reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0;
00490 reco_dircosz = 0;
00491 reco_vtxx = 0; reco_vtxy = 0; reco_vtxz = 0;
00492 evtlength = 0; trklength = 0; trkphfrac = 0; trkphplane = 0;
00493 trkmomrange = 0; trkqp = 0; trkeqp = 0;
00494 reco_dircosy=0; reco_azimuth=0; is_pcev=0;
00495
00496
00497
00498 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00499 detector = 2;
00500 totph=eventSummary->ph.sigcor;
00501
00502 is_fid = 1;
00503 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||
00504 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||
00505 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
00506 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00507 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
00508 +(ntpEvent->vtx.y*ntpEvent->vtx.y))<0.4) is_fid = 0;
00509 }
00510 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00511 detector = 1;
00512
00513
00514 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00515 LoadTrack(ii);
00516 totph+=ntpTrack->ph.sigcor;
00517 }
00518 trkcount+=ntrack;
00519
00520 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00521 LoadShower(ii);
00522 totph+=ntpShower->ph.sigcor;
00523 }
00524 shwcount+=nshower;
00525
00526
00527 is_fid = 1;
00528 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>6.5 ||
00529 sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3)) +
00530 (ntpEvent->vtx.y*ntpEvent->vtx.y))>1) is_fid = 0;
00531 }
00532
00533
00534
00535
00536
00537 is_fidrk=0;
00538
00539 if(LoadTrack(event)) {
00540
00541
00542 trkvtxx=ntpTrack->vtx.x;
00543 trkvtxy=ntpTrack->vtx.y;
00544 trkvtxz=ntpTrack->vtx.z;
00545
00546
00547 trkendx=ntpTrack->end.x;
00548 trkendy=ntpTrack->end.y;
00549 trkendz=ntpTrack->end.z;
00550
00551
00552
00553
00554 if (trkendz<16 && sqrt(pow(trkendx,2)+pow(trkendy,2))>0.4 &&
00555 trkendx<2.7 && trkendx>-1.65 && trkendy<1.65 && trkendy>-1.65 &&
00556 trkendy>(-trkendx)-1.65 && trkendy<trkendx+1.65 &&
00557 trkendy<(-trkendx)+3.55 && trkendy>trkendx-3.55 ) is_fidrk=1;
00558
00559
00560
00561 }
00562
00563
00564 if(LoadTruthForRecoTH(i,mcevent)) {
00565 Float_t *vtx = TrueVtx(mcevent);
00566 Float_t *nu3mom = TrueNuMom(mcevent);
00567 Float_t *targ4vec = Target4Vector(mcevent);
00568
00569 is_mc = 1;
00570 nucleus = Nucleus(mcevent);
00571 flavour = Flavour(mcevent);
00572 initial_state = Initial_state(mcevent);
00573 had_fs = HadronicFinalState(mcevent);
00574 process = IResonance(mcevent);
00575 cc_nc = IAction(mcevent);
00576 true_enu = TrueNuEnergy(mcevent);
00577 true_pxnu = nu3mom[0];
00578 true_pynu = nu3mom[1];
00579 true_pznu = nu3mom[2];
00580 true_emu = TrueMuEnergy(mcevent);
00581 true_eshw = TrueShwEnergy(mcevent);
00582 true_x = X(mcevent);
00583 true_y = Y(mcevent);
00584 true_q2 = Q2(mcevent);
00585 true_w2 = W2(mcevent);
00586 true_dircosz = TrueMuDCosZ(mcevent);
00587 true_dircosneu = TrueMuDCosNeu(mcevent);
00588 true_vtxx = vtx[0];
00589 true_vtxy = vtx[1];
00590 true_vtxz = vtx[2];
00591 true_etgt = targ4vec[3];
00592 true_pxtgt = targ4vec[0];
00593 true_pytgt = targ4vec[1];
00594 true_pztgt = targ4vec[2];
00595
00596 if(gnumi.Status()) {
00597 if(isST) ;
00598 else gnumi.GetParent(mcrecord,mcevent,*nuparent);
00599 }
00600
00601 delete [] vtx;
00602 delete [] nu3mom;
00603 delete [] targ4vec;
00604 }
00605
00606
00607
00608
00609
00610
00611
00612 if(PassBasicCuts(event)) {
00613
00614
00615
00616
00617
00618 pass_basic = 1;
00619 reco_vtxx = ntpEvent->vtx.x;
00620 reco_vtxy = ntpEvent->vtx.y;
00621 reco_vtxz = ntpEvent->vtx.z;
00622 evtlength = ntpEvent->plane.n;
00623
00624 int track_index = -1;
00625 LoadLargestTrackFromEvent(i,track_index);
00626 int shower_index = -1;
00627 LoadLargestShowerFromEvent(i,shower_index);
00628
00629
00630 reco_emu = RecoMuEnergy(0,track_index);
00631 reco_eshw = RecoShwEnergy(shower_index);
00632 reco_eshw_sqrt = RecoShwEnergySqrt(shower_index);
00633 reco_enu = reco_emu + reco_eshw;
00634 reco_qe_enu = RecoQENuEnergy(track_index);
00635 if (reco_enu>0) {
00636 reco_y = reco_eshw/reco_enu;
00637 }
00638 reco_qe_q2 = RecoQEQ2(track_index);
00639 reco_dircosz = RecoMuDCosZ(track_index);
00640 reco_dircosneu = RecoMuDCosNeu(track_index);
00641 is_cev = IsFidAll(track_index);
00642
00643 is_pcev = MyIsFidVtxrz(track_index);
00644
00645 if(track_index!=-1){
00646 trklength = ntpTrack->plane.n;
00647 trkmomrange = ntpTrack->momentum.range;
00648 trkqp = ntpTrack->momentum.qp;
00649 trkeqp = ntpTrack->momentum.eqp;
00650 if (totph>0) {
00651 trkphfrac = ntpTrack->ph.sigcor/totph;
00652 }
00653 if(ntpTrack->plane.n>0) {
00654 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n;
00655 }
00656 }
00657 else {
00658 trklength = 0;
00659 trkmomrange = 0;
00660 trkqp = 0;
00661 trkeqp = 0;
00662 trkphfrac = 0;
00663 trkphplane = 0;
00664 }
00665 pid0 = PID(event,0);
00666
00667
00668 reco_dircosy= RecoMuDCosY(track_index);
00669
00670 reco_azimuth=RecoMuAZM(track_index);
00671
00672
00673
00674
00675
00676
00677
00678 if(PassAnalysisCuts(event)) pass = 1;
00679 }
00680
00681
00682 if(is_fid==1&&is_pcev==1&&pass_basic==1){
00683 ofs << ntpHeader->GetSnarl() << endl;
00684 }
00685
00686
00687
00688
00689 tree->Fill();
00690 }
00691 }
00692 delete nuparent;
00693
00694 save.cd();
00695 tree->Write();
00696 save.Write();
00697 save.Close();
00698 }
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708 void MadEdAnalysis::MyCreatePANData(std::string tag){
00709
00710 std::string savename = "PANData_" + tag + ".root";
00711 TFile save(savename.c_str(),"RECREATE");
00712 save.cd();
00713
00714 GnumiInterface gnumi;
00715 if(!gnumi.Status()) {
00716 cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00717 << " Will not be filling NuParent info"
00718 << endl;
00719 cout << "Set environmental variable $GNUMIAUX to point to the "
00720 << "directory containing the gnumi files to fix this!"
00721 << endl;
00722
00723 }
00724
00725
00726
00727
00728
00729 Int_t detector;
00730 Int_t run;
00731 Int_t snarl;
00732 Int_t nevent;
00733 Int_t event;
00734 Int_t mcevent;
00735 Int_t ntrack;
00736 Int_t nshower;
00737
00738 Int_t is_fid;
00739 Int_t is_fidrk;
00740
00741 Int_t is_cev;
00742 Int_t is_mc;
00743 Int_t pass_basic;
00744 Float_t pid0;
00745 Float_t pid1;
00746 Float_t pid2;
00747 Float_t Pidvalsk;
00748
00749 Int_t pass;
00750
00751 Float_t reco_enu;
00752 Float_t reco_emu;
00753 Float_t reco_eshw;
00754 Float_t reco_eshw_sqrt;
00755 Float_t reco_qe_enu;
00756 Float_t reco_qe_q2;
00757 Float_t reco_y;
00758 Float_t reco_dircosneu;
00759 Float_t reco_dircosz;
00760 Float_t reco_vtxx;
00761 Float_t reco_vtxy;
00762 Float_t reco_vtxz;
00763
00764 Float_t evtlength;
00765 Float_t trklength;
00766 Float_t trkmomrange;
00767 Float_t trkqp;
00768 Float_t trkeqp;
00769 Float_t trkphfrac;
00770 Float_t trkphplane;
00771
00772
00773 TTree *tree = new TTree("pan","pan");
00774
00775
00776 tree->Branch("detector",&detector,"detector/I",32000);
00777 tree->Branch("run",&run,"run/I",32000);
00778 tree->Branch("snarl",&snarl,"snarl/I",32000);
00779 tree->Branch("event",&event,"event/I",32000);
00780 tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
00781 tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
00782 tree->Branch("nshower",&nshower,"nshower/I",32000);
00783
00784 tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
00785 tree->Branch("is_fidrk",&is_fidrk,"is_fidrk/I",32000);
00786
00787 tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
00788 tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
00789 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
00790 tree->Branch("Pidvalsk",&Pidvalsk,"Pidvalsk/F",32000);
00791 tree->Branch("pid0",&pid0,"pid0/F",32000);
00792
00793
00794 tree->Branch("pass",&pass,"pass/I",32000);
00795
00796 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
00797 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
00798 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
00799 tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",32000);
00800 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
00801 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);
00802 tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
00803 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
00804 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);
00805 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
00806 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
00807 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
00808
00809 tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
00810 tree->Branch("trklength",&trklength,"trklength/F",32000);
00811 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
00812 tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
00813 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
00814 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
00815 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
00816
00817 for(int i=0;i<Nentries;i++){
00818 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
00819 << "% done" << std::endl;
00820
00821 if(!this->GetEntry(i)) continue;
00822
00823 nevent = eventSummary->nevent;
00824 if(nevent==0) continue;
00825
00826 run = ntpHeader->GetRun();
00827 snarl = ntpHeader->GetSnarl();
00828
00829 Int_t trkcount=0;
00830 Int_t shwcount=0;
00831 Float_t totph=0;
00832
00833
00834 for(int i=0;i<nevent;i++){
00835 if(!LoadEvent(i)) continue;
00836
00837 totph = 0;
00838
00839
00840 event = i;
00841 ntrack = ntpEvent->ntrack;
00842 nshower = ntpEvent->nshower;
00843
00844
00845
00846 detector = -1; mcevent = -1;
00847 is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0;
00848 pid0 = 0; pid1 = 0; pid2 = 0; pass = 0;
00849 reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0;
00850 reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0;
00851 reco_dircosz = 0;
00852 reco_vtxx = 0; reco_vtxy = 0; reco_vtxz = 0;
00853 evtlength = 0; trklength = 0; trkphfrac = 0; trkphplane = 0;
00854 trkmomrange = 0; trkqp = 0; trkeqp = 0;
00855
00856
00857 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00858 detector = 2;
00859 totph=eventSummary->ph.sigcor;
00860
00861 is_fid = 1;
00862 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||
00863 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||
00864 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
00865 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00866 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
00867 +(ntpEvent->vtx.y*ntpEvent->vtx.y))<0.4) is_fid = 0;
00868 }
00869 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00870 detector = 1;
00871
00872
00873 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00874 LoadTrack(ii);
00875 totph+=ntpTrack->ph.sigcor;
00876 }
00877 trkcount+=ntrack;
00878
00879 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00880 LoadShower(ii);
00881 totph+=ntpShower->ph.sigcor;
00882 }
00883 shwcount+=nshower;
00884
00885
00886 is_fid = 1;
00887 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>6.5 ||
00888 sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3)) +
00889 (ntpEvent->vtx.y*ntpEvent->vtx.y))>1) is_fid = 0;
00890 }
00891
00892
00893
00894 Float_t trkendx=0;
00895 Float_t trkendy=0;
00896 Float_t trkendz=0;
00897
00898
00899 is_fidrk=0;
00900
00901 if(LoadTrack(event)) {
00902 trkendx=ntpTrack->end.x;
00903 trkendy=ntpTrack->end.y;
00904 trkendz=ntpTrack->end.z;
00905
00906 if (trkendz<16 && sqrt(pow(trkendx,2)+pow(trkendy,2))>0.4 &&
00907 trkendx<2.7 && trkendx>-1.65 && trkendy<1.65 && trkendy>-1.65 &&
00908 trkendy>(-trkendx)-1.65 && trkendy<trkendx+1.65 &&
00909 trkendy<(-trkendx)+3.55 && trkendy>trkendx-3.55 ) is_fidrk=1;
00910
00911
00912
00913 }
00914
00915
00916
00917
00918
00919
00920 if(PassBasicCuts(event)) {
00921
00922
00923
00924
00925
00926 pass_basic = 1;
00927 reco_vtxx = ntpEvent->vtx.x;
00928 reco_vtxy = ntpEvent->vtx.y;
00929 reco_vtxz = ntpEvent->vtx.z;
00930 evtlength = ntpEvent->plane.n;
00931
00932 int track_index = -1;
00933 LoadLargestTrackFromEvent(i,track_index);
00934 int shower_index = -1;
00935 LoadLargestShowerFromEvent(i,shower_index);
00936
00937
00938 reco_emu = RecoMuEnergy(0,track_index);
00939 reco_eshw = RecoShwEnergy(shower_index);
00940
00941 reco_enu = reco_emu + reco_eshw;
00942 reco_qe_enu = RecoQENuEnergy(track_index);
00943 if (reco_enu>0) {
00944 reco_y = reco_eshw/reco_enu;
00945 }
00946 reco_qe_q2 = RecoQEQ2(track_index);
00947 reco_dircosz = RecoMuDCosZ(track_index);
00948 reco_dircosneu = RecoMuDCosNeu(track_index);
00949 is_cev = IsFidAll(track_index);
00950
00951 if(track_index!=-1){
00952 trklength = ntpTrack->plane.n;
00953 trkmomrange = ntpTrack->momentum.range;
00954 trkqp = ntpTrack->momentum.qp;
00955 trkeqp = ntpTrack->momentum.eqp;
00956 if (totph>0) {
00957 trkphfrac = ntpTrack->ph.sigcor/totph;
00958 }
00959 if(ntpTrack->plane.n>0) {
00960 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n;
00961 }
00962 }
00963 else {
00964 trklength = 0;
00965 trkmomrange = 0;
00966 trkqp = 0;
00967 trkeqp = 0;
00968 trkphfrac = 0;
00969 trkphplane = 0;
00970 }
00971 pid0 = PID(event,0);
00972
00973
00974
00975
00976 if(PassAnalysisCuts(event)) pass = 1;
00977 }
00978
00979 tree->Fill();
00980 }
00981 }
00982
00983 save.cd();
00984 tree->Write();
00985 save.Write();
00986 save.Close();
00987 }
00988
00989
00990
00991
00992 void MadEdAnalysis::MyMakeMyFile(std::string tag){
00993
00994 std::string savename = "DPHistos_" + tag + ".root";
00995 TFile save(savename.c_str(),"RECREATE");
00996 save.cd();
00997
00998
00999
01000
01001 TH1F *evlength_nc = new TH1F("Event length - nc",
01002 "Event length - nc",
01003 50,0,500);
01004 evlength_nc->SetXTitle("Event length (planes)");
01005 TH1F *evlength_cc = new TH1F("Event length - cc",
01006 "Event length - cc",
01007 50,0,500);
01008 evlength_cc->SetXTitle("Event length (planes)");
01009
01010
01011 TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
01012 "Track ph frac - nc",
01013 50,0,1);
01014 phfrac_nc->SetXTitle("pulse height fraction");
01015
01016
01017 TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
01018 "Track ph frac - cc",
01019 50,0,1);
01020 phfrac_cc->SetXTitle("pulse height fraction");
01021
01022
01023 TH1F *phplane_nc = new TH1F("ph per plane (nc)",
01024 "ph per plane (nc)",
01025 50,0,5000);
01026 phplane_nc->SetXTitle("pulse height per plane");
01027 TH1F *phplane_cc = new TH1F("ph per plane (cc)",
01028 "ph per plane (cc)",
01029 50,0,5000);
01030 phplane_cc->SetXTitle("pulse height per plane");
01031
01032
01033
01034
01035
01036
01037 for(int i=0;i<Nentries;i++){
01038
01039
01040 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
01041 << "% done" << std::endl;
01042 if(!this->GetEntry(i)) continue;
01043
01044 Int_t nevent = eventSummary->nevent;
01045
01046
01047
01048
01049 if(nevent==0) continue;
01050
01051 Int_t trkcount=0;
01052 Int_t shwcount=0;
01053
01054
01055
01056 for(int j=0;j<eventSummary->nevent;j++){
01057
01058 if(!LoadEvent(j)) continue;
01059
01060
01061
01062
01063 Int_t ntrack = 0;
01064 ntrack=ntpEvent->ntrack;
01065 Int_t nshower = 0;
01066 nshower=ntpEvent->nshower;
01067
01068 Float_t totph=0;
01069
01070
01071
01072 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
01073 LoadTrack(ii);
01074
01075 totph+=ntpTrack->ph.sigcor;
01076 LoadTHTrack(ii);
01077
01078 }
01079
01080 trkcount+=ntrack;
01081
01082
01083 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
01084 LoadShower(ii);
01085
01086 totph+=ntpShower->ph.sigcor;
01087 }
01088
01089 shwcount+=nshower;
01090
01091
01092
01093 Int_t cc_nc = 0;
01094 Int_t flavour = 0;
01095 Int_t detector = -1;
01096 Int_t is_fid = 0;
01097
01098
01099
01100 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01101 detector = 2;
01102
01103 is_fid = 1;
01104 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||
01105 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||
01106 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
01107 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0;
01108 }
01109 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01110 detector = 1;
01111
01112 is_fid = 1;
01113 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>6.5 ||
01114 sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3)) +
01115 (ntpEvent->vtx.y*ntpEvent->vtx.y))>1) is_fid = 0;
01116 }
01117
01118
01119
01120 Int_t mcevent=-1;
01121
01122 if(LoadTruthForRecoTH(j,mcevent)) {
01123
01124
01125 cc_nc = IAction(mcevent);
01126 flavour = Flavour(mcevent);
01127
01128 }
01129
01130 int track_index = -1;
01131 LoadLargestTrackFromEvent(j,track_index);
01132
01133 if (track_index!=-1) {
01134
01135
01136
01137 Int_t trkpass=ntpTrack->fit.pass;
01138
01139
01140 if (is_fid && trkpass) {
01141
01142
01143
01144
01145 Float_t evlength=0;
01146 Float_t phfrac=0;
01147 Float_t phplane=0;
01148
01149 evlength=ntpEvent->plane.n;
01150
01151 Float_t phtrack=ntpTrack->ph.sigcor;
01152 Float_t phevent=ntpEvent->ph.sigcor;
01153 if (totph>0) {phfrac=phtrack/phevent;}
01154
01155 Float_t trkplane=ntpTrack->plane.n;
01156 if (trkplane>0) {phplane=phtrack/trkplane;}
01157
01158
01159 if(flavour==2 && cc_nc==1) {
01160 evlength_cc->Fill(evlength);
01161 phfrac_cc->Fill(phfrac);
01162 phplane_cc->Fill(phplane);
01163 }
01164
01165 else if(cc_nc==0) {
01166
01167
01168
01169
01170 evlength_nc->Fill(evlength);
01171 phfrac_nc->Fill(phfrac);
01172 phplane_nc->Fill(phplane);
01173 }
01174
01175
01176 }
01177 }
01178 }
01179 }
01180
01181 save.Write();
01182 save.Close();
01183 }
01184
01185
01186
01187
01188 void MadEdAnalysis::MyReadPIDFile(std::string tag){
01189
01190 fLikeliFile = new TFile(tag.c_str(),"READ");
01191
01192 if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) {
01193
01194 fLikeliFile->cd();
01195
01196 fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
01197 fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
01198 fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
01199 fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
01200 fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
01201 fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
01202
01203 for(int i=0;i<6;i++){
01204 float integ = fLikeliHist[i]->Integral();
01205 fLikeliHist[i]->Scale(1./integ);
01206 }
01207 }
01208 else fLikeliFile = NULL;
01209 }
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219 Float_t MadEdAnalysis::EvtLength() {
01220 if(LoadEvent(0)) return ntpEvent->plane.n;
01221 return 0;
01222 }
01223
01224
01225 Float_t MadEdAnalysis::Trkphsig(Int_t itrk) {
01226 if(LoadTrack(itrk)) return ntpTrack->ph.sigcor;
01227 return 0.;
01228 }
01229
01230
01231 Float_t MadEdAnalysis::Trkplanes(Int_t itrk) {
01232 if(LoadTrack(itrk)) return ntpTrack->plane.n;
01233 return 0.;
01234 }
01235
01236
01237 Float_t MadEdAnalysis::Evtphsig() {
01238 if(LoadEvent(0)) return ntpEvent->ph.sigcor;
01239 return 0;
01240 }
01241
01242
01243 Float_t MadEdAnalysis::HitF(Int_t itr){
01244 if(!LoadTrack(itr)) return 0;
01245
01246 float ntpE = float(ntpEvent->nstrip);
01247 if(ntpE==0) ntpE=0.00001;
01248
01249 return float(ntpTrack->nstrip)/ntpE;
01250
01251 }
01252
01253
01254 Float_t MadEdAnalysis::ETrkF(Int_t itr){
01255 if(!LoadTrack(itr)) return 0.;
01256
01257 float ntpES = float(ntpEvent->ph.sigcor);
01258 if(ntpES==0) ntpES=0.00001;
01259
01260 return float(ntpTrack->ph.sigcor)/ntpES;
01261
01262 }
01263
01264
01265 Float_t MadEdAnalysis::DeDx(Int_t itr){
01266
01267 float ntpED = ntpTrack->ds;
01268 if(LoadTrack(itr)) return ntpTrack->ph.sigcor/(100.*ntpED);
01269 return 0.;
01270 }
01271
01272
01273 TMultiLayerPerceptron* MadEdAnalysis::NeuNetTrain(TFile* f) {
01274
01275 TTree *treeN = (TTree*) f->Get("NN");
01276
01277 TEventList *elist1= new TEventList("elist1","eventlist1",0,0);
01278 TEventList *elist2= new TEventList("elist2","eventlist3",0,0);
01279 TEventList *elist3= new TEventList("elist3","eventlist3",0,0);
01280 TEventList *elist4= new TEventList("elist4","eventlist4",0,0);
01281
01282
01283 treeN->Draw(">>elist1","flavor==1","",652,0);
01284 treeN->Draw(">>elist2","flavor==0","",2127,0);
01285
01286 elist1->Add(elist2);
01287
01288
01289 treeN->Draw(">>elist3","flavor==1","",659,652);
01290 treeN->Draw(">>elist4","flavor==0","",1656,2127);
01291 elist3->Add(elist4);
01292
01293
01294
01295 TMultiLayerPerceptron *neural;
01296 neural = new TMultiLayerPerceptron("evlength,phfrac,phplane:8:flavor",
01297 treeN,elist1,elist3);
01298 neural->Train(1000,"text, graph, update=100");
01299
01300 return neural;
01301
01302 }
01303
01304
01305 Float_t MadEdAnalysis::NeuNetEval(Int_t event) {
01306
01307 if(!fneural) return 0;
01308 if(!LoadEvent(event)) return 0;
01309 Int_t track = -1;
01310 if(!LoadLargestTrackFromEvent(event,track)) return 0;
01311
01312
01313 Float_t evtphsig=0;
01314 Float_t trkplanes=0;
01315
01316 evtphsig=this->Evtphsig();
01317 trkplanes=this->Trkplanes(track);
01318
01319 if(this->Evtphsig()==0) evtphsig=0.000001;
01320 if(this->Trkplanes(track)==0) trkplanes=0.00001;
01321
01322 Double_t params[3];
01323 params[0]=this->EvtLength();
01324 params[1]=this->Trkphsig(track)/evtphsig;
01325 params[2]=this->Trkphsig(track)/trkplanes;
01326
01327 return fneural->Evaluate(0,params);
01328 }
01329
01330
01331 void MadEdAnalysis::PIDHisto() {
01332
01333
01334
01335 TFile save("PidHisto","RECREATE");
01336 save.cd();
01337
01338
01339 TH1F *EvtlengthCC = new TH1F("EvtlengthCC","No of planes in Event - #nu_{#mu} CC Events",100,0,500);
01340
01341 TH1F *EvtlengthNC = new TH1F("EvtlengthNC","No of planes in Event - #nu_{#mu} NC Events",100,0,500);
01342
01343 TH1F *FracTrkphEvtphCC = new TH1F("FracTrkphEvtphCC","Pulse height of track/Pulse height of Event - #nu_{#mu} CC Events",100,0,1.);
01344
01345 TH1F *FracTrkphEvtphNC = new TH1F("FracTrkphEvtphNC","Pulse height of track/Pulse height of Event - #nu_{#mu} NC Events",100,0,1.);
01346
01347 TH1F *TrkphperplaneCC = new TH1F("TrkphperplaneCC","Pulse height of track/Number of planes in track - #nu_{#mu} CC Events",500,0.,2000.);
01348
01349 TH1F *TrkphperplaneNC = new TH1F("TrkphperplaneNC","Pulse height of track/Number of planes in track - #nu_{#mu} NC Events",500,0.,2000.);
01350
01351
01352 for(int i=0;i<Nentries;i++){
01353
01354 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
01355 << "% done" << std::endl;
01356
01357 this->GetEntry(i);
01358
01359
01360
01361 if(eventSummary->nevent!=1)continue;
01362
01363 for(int i=0;i<eventSummary->nevent;i++){
01364
01365
01366
01367 if(!LoadEvent(i)) continue;
01368 int mc_entry = 0;
01369 if(!LoadTruthForRecoTH(i,mc_entry)) continue;
01370
01371
01372
01373 int inu = INu(mc_entry);
01374 int iaction = IAction(mc_entry);
01375
01376
01377
01378
01379
01380
01381
01382 if(abs(inu)==14&&iaction==1&&ntpEvent->ntrack==1&&this->PassCuts(0)){
01383
01384
01385 EvtlengthCC->Fill(this->EvtLength());
01386 if(this->Evtphsig()==0.) {FracTrkphEvtphCC->Fill(0.);}
01387 else {FracTrkphEvtphCC->Fill(this->Trkphsig(0)/this->Evtphsig());}
01388 if(this->Trkplanes(0)==0.) { TrkphperplaneCC->Fill(0.);}
01389 else {TrkphperplaneCC->Fill(this->Trkphsig(0)/this->Trkplanes(0));}
01390 }
01391
01392 else if(iaction==0){
01393
01394 EvtlengthNC->Fill(this->EvtLength());
01395 if(this->Evtphsig()==0.) {FracTrkphEvtphNC->Fill(0.);}
01396 else {FracTrkphEvtphNC->Fill(this->Trkphsig(0)/this->Evtphsig());}
01397 if(this->Trkplanes(0)==0.) { TrkphperplaneNC->Fill(0.);}
01398 else {TrkphperplaneNC->Fill(this->Trkphsig(0)/this->Trkplanes(0));}
01399 }
01400
01401 }
01402
01403 }
01404
01405 save.Write();
01406 save.Close();
01407
01408 }
01409
01410
01411 Float_t* MadEdAnalysis::LikeliPID() {
01412
01413 TFile pfile("Pid_R1.9n.root","READ");
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428 Double_t P_fraccc, P_trkphcc;
01429
01430 Float_t *Pid;
01431 Float_t PD[2];
01432 PD[0]=0;
01433 PD[1]=0;
01434 Pid=PD;
01435
01436 TH1F *nEvtlengthCC = (TH1F*) pfile.Get("EvtlengthCC");
01437 TH1F *nEvtlengthNC = (TH1F*) pfile.Get("EvtlengthNC");
01438
01439 TH1F *nFracTrkphEvtphCC = (TH1F*) pfile.Get("FracTrkphEvtphCC");
01440 TH1F *nFracTrkphEvtphNC = (TH1F*) pfile.Get("FracTrkphEvtphNC");
01441
01442 TH1F *nTrkphperplaneCC = (TH1F*) pfile.Get("TrkphperplaneCC");
01443 TH1F *nTrkphperplaneNC = (TH1F*) pfile.Get("TrkphperplaneNC");
01444
01445
01446 Int_t bin_evtcc = nEvtlengthCC->FindBin(this->EvtLength());
01447 Double_t P_evtcc = nEvtlengthCC->GetBinContent(bin_evtcc);
01448
01449 Float_t Evtphsig=this->Evtphsig();
01450 if (Evtphsig==0) Evtphsig=1e-04;
01451
01452 Int_t bin_fraccc = nFracTrkphEvtphCC->FindBin(this->Trkphsig(0)/Evtphsig);
01453 P_fraccc = nFracTrkphEvtphCC->GetBinContent(bin_fraccc);
01454
01455
01456 Float_t Trkplanes =this->Trkplanes(0);
01457 if (Trkplanes==0) Trkplanes=1e-04;
01458
01459 Int_t bin_trkphcc = nTrkphperplaneCC->FindBin(this->Trkphsig(0)/Trkplanes);
01460 P_trkphcc = nTrkphperplaneCC->GetBinContent(bin_trkphcc);
01461
01462 Double_t P_cc=(P_evtcc/nEvtlengthCC->Integral())*(P_fraccc/nFracTrkphEvtphCC->Integral())*(P_trkphcc/nTrkphperplaneCC->Integral());
01463
01464 Double_t P_fracnc, P_trkphnc;
01465
01466 Int_t bin_evtnc = nEvtlengthNC->FindBin(this->EvtLength());
01467 Double_t P_evtnc = nEvtlengthNC->GetBinContent(bin_evtnc);
01468
01469
01470 Int_t bin_fracnc = nFracTrkphEvtphNC->FindBin(this->Trkphsig(0)/Evtphsig);
01471 P_fracnc = nFracTrkphEvtphNC->GetBinContent(bin_fracnc);
01472
01473
01474 Int_t bin_trkphnc = nTrkphperplaneNC->FindBin(this->Trkphsig(0)/Trkplanes);
01475 P_trkphnc = nTrkphperplaneNC->GetBinContent(bin_trkphnc);
01476
01477 Double_t P_nc=(P_evtnc/nEvtlengthNC->Integral())*(P_fracnc/nFracTrkphEvtphNC->Integral())*(P_trkphnc/nTrkphperplaneNC->Integral());
01478
01479
01480
01481 if (P_cc==0) P_cc=1e-04;
01482 if (P_nc==0) P_nc=1e-04;
01483
01484
01485 Double_t Pcc = sqrt(-log(P_cc));
01486 Double_t Pnc = sqrt(-(log(P_nc)));
01487
01488
01489 PD[0] = Pnc-Pcc;
01490 PD[1] = P_cc/(P_cc+P_nc);
01491
01492 return Pid;
01493
01494 }
01495
01496
01497 Float_t* MadEdAnalysis::MyLikeliQE(TH1F **hist){
01498
01499 Float_t *PQE;
01500 Float_t PidParQE[5];
01501 Float_t ProbNC = 1.;
01502 Float_t ProbCC = 1.;
01503 Float_t ProbQE = 1.;
01504 Float_t PidParQECC = 0;
01505 Float_t PidParQENC = 0;
01506
01507 PQE=PidParQE;
01508
01509 Float_t *CCNCVars = this->CCNCSepVars();
01510
01511 Float_t nstp=ntpEvent->nstrip;
01512 Float_t ntds=ntpTrack->ds;
01513 if(nstp==0) nstp=0.00001;
01514 if(ntds==0) ntds=0.00001;
01515
01516
01517 int bin1 = hist[0]->FindBin(float(ntpTrack->nstrip)
01518 /nstp);
01519 ProbNC *= hist[0]->GetBinContent(bin1);
01520
01521 bin1 = hist[1]->FindBin(ntpTrack->ph.sigcor/(100.*ntds));
01522 ProbNC *= hist[1]->GetBinContent(bin1);
01523
01524 bin1 = hist[2]->FindBin(CCNCVars[2]);
01525 ProbNC *= hist[2]->GetBinContent(bin1);
01526
01527
01528
01529 int bin2 = hist[3]->FindBin(float(ntpTrack->nstrip)
01530 /nstp);
01531 ProbCC *= hist[3]->GetBinContent(bin2);
01532
01533 bin2 = hist[4]->FindBin(ntpTrack->ph.sigcor/(100.*ntds));
01534 ProbCC *= hist[4]->GetBinContent(bin2);
01535
01536 bin2 = hist[5]->FindBin(CCNCVars[2]);
01537 ProbCC *= hist[5]->GetBinContent(bin2);
01538
01539
01540
01541 int bin3 = hist[6]->FindBin(float(ntpTrack->nstrip)
01542 /nstp);
01543 ProbQE *= hist[6]->GetBinContent(bin3);
01544
01545 bin3 = hist[7]->FindBin(ntpTrack->ph.sigcor/(100.*ntds));
01546 ProbQE *= hist[7]->GetBinContent(bin3);
01547
01548 bin3 = hist[8]->FindBin(CCNCVars[2]);
01549 ProbQE *= hist[8]->GetBinContent(bin3);
01550
01551
01552 delete CCNCVars;
01553
01554
01555 if(ProbQE!=0&&ProbCC!=0)
01556 PidParQECC = sqrt(-TMath::Log(ProbQE))
01557 -sqrt(-TMath::Log(ProbCC));
01558 else if(ProbQE==0&&ProbCC==0) PidParQECC=10;
01559 else if(ProbQE==0) PidParQECC = 10.-sqrt(-TMath::Log(ProbCC));
01560 else if(ProbCC==0) PidParQECC = sqrt(-TMath::Log(ProbQE))-10.;
01561
01562 if(ProbQE!=0&&ProbNC!=0)
01563 PidParQENC = sqrt(-TMath::Log(ProbQE))
01564 -sqrt(-TMath::Log(ProbNC));
01565 else if(ProbQE==0&&ProbNC==0) PidParQENC=10;
01566 else if(ProbQE==0) PidParQENC = 10.-sqrt(-TMath::Log(ProbNC));
01567 else if(ProbNC==0) PidParQENC = sqrt(-TMath::Log(ProbQE))-10.;
01568
01569
01570 PidParQE[0]=PidParQECC;
01571 PidParQE[1]=PidParQENC;
01572
01573 PidParQE[2]=ProbNC;
01574 PidParQE[3]=ProbCC;
01575 PidParQE[4]=ProbQE;
01576
01577
01578 return PQE;
01579 }
01580
01581
01582
01583
01584
01585
01586 void MadEdAnalysis::MyMakeQEFile(std::string tag){
01587
01588 std::string savename = "PidQE_" + tag + ".root";
01589 TFile save(savename.c_str(),"RECREATE");
01590 save.cd();
01591
01592
01593 TH1F *NEvent_NC = new TH1F("NEvent_NC",
01594 "Number of Reco'd Events - NC",
01595 20,0,20);
01596 NEvent_NC->SetXTitle("Number of Events");
01597 TH1F *NEvent_CC = new TH1F("NEvent_CC",
01598 "Number of Reco'd Events - CC",
01599 20,0,20);
01600 NEvent_CC->SetXTitle("Number of Events");
01601 TH1F *NEvent_QE = new TH1F("NEvent_QE",
01602 "Number of Reco'd Events - QE",
01603 20,0,20);
01604 NEvent_QE->SetXTitle("Number of Events");
01605
01606
01607 TH1F *NShower_NC = new TH1F("NShower_NC",
01608 "Number of Reco'd Showers - NC",
01609 20,0,20);
01610 NShower_NC->SetXTitle("Number of Showers");
01611 TH1F *NShower_CC = new TH1F("NShower_CC",
01612 "Number of Reco'd Showers - CC",
01613 20,0,20);
01614 NShower_CC->SetXTitle("Number of Showers");
01615 TH1F *NShower_QE = new TH1F("NShower_QE",
01616 "Number of Reco'd Showers - QE",
01617 20,0,20);
01618 NShower_QE->SetXTitle("Number of Showers");
01619
01620
01621 TH1F *NTrack_NC = new TH1F("NTrack_NC",
01622 "Number of Reco'd Tracks - NC",
01623 20,0,20);
01624 NTrack_NC->SetXTitle("Number of Tracks");
01625 TH1F *NTrack_CC = new TH1F("NTrack_CC",
01626 "Number of Reco'd Tracks - CC",
01627 20,0,20);
01628 NTrack_CC->SetXTitle("Number of Tracks");
01629 TH1F *NTrack_QE = new TH1F("NTrack_QE",
01630 "Number of Reco'd Tracks - QE",
01631 20,0,20);
01632 NTrack_QE->SetXTitle("Number of Tracks");
01633
01634
01635 TH2F *NShower_Vs_NTrack_NC =
01636 new TH2F("NShower_Vs_NTrack_NC",
01637 "Number of Reco'd Showers Vs Tracks- NC",
01638 20,0,5,20,0,5);
01639 NShower_Vs_NTrack_NC->SetXTitle("Number of Tracks");
01640 NShower_Vs_NTrack_NC->SetYTitle("Number of Showers");
01641 TH2F *NShower_Vs_NTrack_CC =
01642 new TH2F("NShower_Vs_NTrack_CC",
01643 "Number of Reco'd Showers Vs Tracks- CC",
01644 20,0,5,20,0,5);
01645 NShower_Vs_NTrack_CC->SetXTitle("Number of Tracks");
01646 NShower_Vs_NTrack_CC->SetYTitle("Number of Showers");
01647 TH2F *NShower_Vs_NTrack_QE =
01648 new TH2F("NShower_Vs_NTrack_QE",
01649 "Number of Reco'd Showers Vs Tracks- QE",
01650 20,0,5,20,0,5);
01651 NShower_Vs_NTrack_QE->SetXTitle("Number of Tracks");
01652 NShower_Vs_NTrack_QE->SetYTitle("Number of Showers");
01653
01654
01655 TH1F *NHitTrkFrac_NC=new TH1F("NHitTrkFrac_NC",
01656 "Fraction of Event Hits in Primary Track - NC",
01657 100,0,1);
01658 NHitTrkFrac_NC->SetXTitle("Hit Fraction");
01659 TH1F *NHitTrkFrac_CC=new TH1F("NHitTrkFrac_CC",
01660 "Fraction of Event Hits in Primary Track - CC",
01661 100,0,1);
01662 NHitTrkFrac_CC->SetXTitle("Hit Fraction");
01663 TH1F *NHitTrkFrac_QE=new TH1F("NHitTrkFrac_QE",
01664 "Fraction of Event Hits in Primary Track - QE",
01665 100,0,1);
01666 NHitTrkFrac_QE->SetXTitle("Hit Fraction");
01667
01668
01669
01670
01671
01672
01673 TH1F *ETrkFrac_NC=new TH1F("ETrkFrac_NC",
01674 "Fraction of Event SigCor in Primary Track - NC",
01675 100,0,1);
01676 ETrkFrac_NC->SetXTitle("SigCor Fraction");
01677 TH1F *ETrkFrac_CC=new TH1F("ETrkFrac_CC",
01678 "Fraction of Event SigCor in Primary Track - CC",
01679 100,0,1);
01680 ETrkFrac_CC->SetXTitle("SigCor Fraction");
01681 TH1F *ETrkFrac_QE=new TH1F("ETrkFrac_QE",
01682 "Fraction of Event SigCor in Primary Track - QE",
01683 100,0,1);
01684 ETrkFrac_QE->SetXTitle("SigCor Fraction");
01685
01686
01687 TH1F *TrkMomFrac_NC = new TH1F("TrkMomFrac_NC",
01688 "Event Momentum Fraction in Track - NC",
01689 100,0,1);
01690 TrkMomFrac_NC->SetXTitle("Momentum (GeV)");
01691 TH1F *TrkMomFrac_CC = new TH1F("TrkMomFrac_CC",
01692 "Event Momentum Fraction in Track - CC",
01693 100,0,1);
01694 TrkMomFrac_CC->SetXTitle("Momentum (GeV)");
01695 TH1F *TrkMomFrac_QE = new TH1F("TrkMomFrac_QE",
01696 "Event Momentum Fraction in Track - QE",
01697 100,0,1);
01698 TrkMomFrac_QE->SetXTitle("Momentum (GeV)");
01699
01700
01701 TH1F *TrkMomFrac0_NC = new TH1F("TrkMomFrac0_NC",
01702 "Event Momentum Fraction in Track - NC",
01703 100,0,1);
01704 TrkMomFrac0_NC->SetXTitle("Momentum (GeV)");
01705 TH1F *TrkMomFrac0_CC = new TH1F("TrkMomFrac0_CC",
01706 "Event Momentum Fraction in Track - CC",
01707 100,0,1);
01708 TrkMomFrac0_CC->SetXTitle("Momentum (GeV)");
01709 TH1F *TrkMomFrac0_QE = new TH1F("TrkMomFrac0_QE",
01710 "Event Momentum Fraction in Track - QE",
01711 100,0,1);
01712 TrkMomFrac0_QE->SetXTitle("Momentum (GeV)");
01713
01714
01715 TH1F *VtxShwEnergy_NC = new TH1F("VtxShwEnergy_NC",
01716 "Vertex Shower Energy - NC",
01717 1000,0,100);
01718 VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
01719 TH1F *VtxShwEnergy_CC = new TH1F("VtxShwEnergy_CC",
01720 "Vertex Shower Energy - CC",
01721 1000,0,100);
01722 VtxShwEnergy_CC->SetXTitle("Energy (GeV)");
01723 TH1F *VtxShwEnergy_QE = new TH1F("VtxShwEnergy_QE",
01724 "Vertex Shower Energy - QE",
01725 1000,0,100);
01726 VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
01727
01728
01729 TH1F *dsdz_NC = new TH1F("dsdz_NC","Primary Track dsdz - NC",100,-1,20);
01730 dsdz_NC->SetXTitle("dsdz");
01731 TH1F *dsdz_CC = new TH1F("dsdz_CC","Primary Track dsdz - CC",100,-1,20);
01732 dsdz_CC->SetXTitle("dsdz");
01733 TH1F *dsdz_QE = new TH1F("dsdz_QE","Primary Track dsdz - QE",100,-1,20);
01734 dsdz_QE->SetXTitle("dsdz");
01735
01736 TH1F *NHitPlanes_NC = new TH1F("NHitPlanes_NC","Number of Hit Planes - NC",
01737 500,-0.5,499.5);
01738 NHitPlanes_NC->SetXTitle("Number of Planes");
01739 TH1F *NHitPlanes_CC = new TH1F("NHitPlanes_CC","Number of Hit Planes - CC",
01740 500,-0.5,499.5);
01741 NHitPlanes_NC->SetXTitle("Number of Planes");
01742 TH1F *NHitPlanes_QE = new TH1F("NHitPlanes_QE","Number of Hit Planes - QE",
01743 500,-0.5,499.5);
01744 NHitPlanes_QE->SetXTitle("Number of Planes");
01745
01746 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_NC
01747 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_NC",
01748 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - NC",1000,0,50,100,0,1);
01749 NHitTrkFrac_Vs_VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
01750 NHitTrkFrac_Vs_VtxShwEnergy_NC->SetYTitle("Hit Fraction");
01751 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_CC
01752 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_CC",
01753 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - CC",1000,0,50,100,0,1);
01754 NHitTrkFrac_Vs_VtxShwEnergy_CC->SetXTitle("Energy (GeV)");
01755 NHitTrkFrac_Vs_VtxShwEnergy_CC->SetYTitle("Hit Fraction");
01756 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_QE
01757 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_QE",
01758 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - QE",1000,0,50,100,0,1);
01759 NHitTrkFrac_Vs_VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
01760 NHitTrkFrac_Vs_VtxShwEnergy_QE->SetYTitle("Hit Fraction");
01761
01762
01763 TH2F *NTrack_Vs_EvSC_NC
01764 = new TH2F("NTrack_Vs_EvSC_NC",
01765 "Number of Reco'd Tracks vs Event SigCor - NC",
01766 1000,0,100000,20,0,20);
01767 NTrack_Vs_EvSC_NC->SetXTitle("Event SigCor");
01768 NTrack_Vs_EvSC_NC->SetYTitle("Number of Tracks");
01769
01770 TH2F *NTrack_Vs_EvSC_CC
01771 = new TH2F("NTrack_Vs_EvSC_CC",
01772 "Number of Reco'd Tracks vs Event SigCor - CC",
01773 1000,0,100000,20,0,20);
01774 NTrack_Vs_EvSC_CC->SetXTitle("Event SigCor");
01775 NTrack_Vs_EvSC_CC->SetYTitle("Number of Tracks");
01776
01777 TH2F *NTrack_Vs_EvSC_QE
01778 = new TH2F("NTrack_Vs_EvSC_QE",
01779 "Number of Reco'd Tracks vs Event SigCor - QE",
01780 1000,0,100000,20,0,20);
01781 NTrack_Vs_EvSC_QE->SetXTitle("Event SigCor");
01782 NTrack_Vs_EvSC_QE->SetYTitle("Number of Tracks");
01783
01784
01785 TH2F *NHitTrkFrac_Vs_Y_CC
01786 = new TH2F("NHitTrkFrac_Vs_Y_CC",
01787 "Fraction of Event Hits in Primary Track Vs Y - CC",
01788 100,0,1,100,0,1);
01789 NHitTrkFrac_Vs_Y_CC->SetXTitle("Y");
01790 NHitTrkFrac_Vs_Y_CC->SetYTitle("Hit Fraction");
01791
01792
01793 TH1F *TrkLen_NC = new TH1F("TrkLen_NC","Track length - NC",50,0,50);
01794 TH1F *TrkLen_CC = new TH1F("TrkLen_CC","Track length - CC",50,0,50);
01795 TH1F *TrkLen_QE = new TH1F("TrkLen_QE","Track length - CC",50,0,50);
01796
01797
01798 TH1F *dedx_NC = new TH1F("dedx_NC","Average dEdx - NC",1000,0,2000);
01799 TH1F *dedx_CC = new TH1F("dedx_CC","Average dEdx - CC",1000,0,2000);
01800 TH1F *dedx_QE = new TH1F("dedx_QE","Average dEdx - QE",1000,0,2000);
01801
01802
01803 TH1F *HalfRatio_NC
01804 = new TH1F("HalfRatio_NC",
01805 "Charge Ratio: First Half/Second Half of Track - NC",
01806 150,-1,14);
01807 TH1F *HalfRatio_CC
01808 = new TH1F("HalfRatio_CC",
01809 "Charge Ratio: First Half/Second Half of Track - CC",
01810 150,-1,14);
01811 TH1F *HalfRatio_QE
01812 = new TH1F("HalfRatio_QE",
01813 "Charge Ratio: First Half/Second Half of Track - QE",
01814 150,-1,14);
01815
01816
01817 TH2F *RangeCurvDiff_Vs_TrkLen_NC
01818 = new TH2F("RangeCurvDiff_Vs_TrkLen_NC",
01819 "Diff in Momentum from Range and Curv vs Track Length - NC",
01820 100,0,30,200,-3,17);
01821 TH2F *RangeCurvDiff_Vs_TrkLen_CC
01822 = new TH2F("RangeCurvDiff_Vs_TrkLen_CC",
01823 "Diff in Momentum from Range and Curv vs Track Length - CC",
01824 100,0,30,200,-3,17);
01825 TH2F *RangeCurvDiff_Vs_TrkLen_QE
01826 = new TH2F("RangeCurvDiff_Vs_TrkLen_QE",
01827 "Diff in Momentum from Range and Curv vs Track Length - QE",
01828 100,0,30,200,-3,17);
01829
01830
01831
01832
01833
01834 Float_t true_enu;
01835 Float_t true_emu;
01836 Float_t true_eshw;
01837 Float_t true_x;
01838 Float_t true_y;
01839 Float_t true_q2;
01840 Float_t true_w2;
01841 Float_t true_dircosneu;
01842 Float_t true_dircosz;
01843 Float_t true_vtxx;
01844 Float_t true_vtxy;
01845 Float_t true_vtxz;
01846 Int_t process = 0;
01847
01848
01849 Int_t detector;
01850 Int_t run;
01851 Int_t snarl;
01852 Int_t nevent;
01853 Int_t event;
01854 Int_t mcevent;
01855 Int_t ntrack;
01856 Int_t nshower;
01857 Int_t inu = 0;
01858 Int_t iaction = 0;
01859 Int_t is_fid = 0;
01860
01861 Int_t is_mc;
01862 Int_t pass_basic;
01863
01864
01865
01866
01867
01868 Float_t reco_emu = 0.;
01869 Float_t reco_eshw = 0.;
01870
01871
01872
01873
01874 Float_t reco_dircosz = 0.;
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887 Int_t mcindex=0;
01888 Int_t recoindex=0;
01889 Int_t trkindex=0;
01890 Int_t shwindex=0;
01891
01892
01893
01894
01895
01896
01897 for(int i=0;i<Nentries;i++){
01898
01899 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
01900 << "% done" << std::endl;
01901 if(!this->GetEntry(i)) continue;
01902
01903
01904
01905 nevent = eventSummary->nevent;
01906 if(nevent==0) continue;
01907
01908 run = ntpHeader->GetRun();
01909 snarl = ntpHeader->GetSnarl();
01910
01911
01912 for(int i=0;i<eventSummary->nevent;i++){
01913 if(!LoadEvent(i)) continue;
01914
01915
01916 event = i;
01917 ntrack = ntpEvent->ntrack;
01918 nshower = ntpEvent->nshower;
01919
01920
01921
01922 detector = -1; mcevent = -1;
01923
01924
01925
01926 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01927 detector = 2;
01928
01929 is_fid = 1;
01930 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||
01931 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||
01932 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
01933 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0;
01934 }
01935 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01936 detector = 1;
01937
01938 is_fid = 1;
01939 if(ntpEvent->vtx.z<0.5 || ntpEvent->end.z>6.5 ||
01940 sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3)) +
01941 (ntpEvent->vtx.y*ntpEvent->vtx.y))>1) is_fid = 0;
01942 }
01943
01944
01945
01946 if(LoadTruthForRecoTH(i,mcevent)) {
01947
01948 mcindex = mcevent;
01949 Float_t *vtx = TrueVtx(mcevent);
01950
01951 is_mc = 1;
01952
01953 process = IResonance(mcevent);
01954 iaction = IAction(mcevent);
01955 inu = INu(mcevent);
01956 true_enu = TrueNuEnergy(mcevent);
01957 true_emu = TrueMuEnergy(mcevent);
01958 true_eshw = TrueShwEnergy(mcevent);
01959 true_x = X(mcevent);
01960 true_y = Y(mcevent);
01961 true_q2 = Q2(mcevent);
01962 true_w2 = W2(mcevent);
01963 true_dircosz = TrueMuDCosZ(mcevent);
01964 true_dircosneu = TrueMuDCosNeu(mcevent);
01965 true_vtxx = vtx[0];
01966 true_vtxy = vtx[1];
01967 true_vtxz = vtx[2];
01968 delete [] vtx;
01969 }
01970
01971
01972 if(PassBasicCuts(event)&&is_fid) {
01973
01974 recoindex = event;
01975 pass_basic = 1;
01976
01977
01978
01979
01980 int track_index = -1;
01981
01982 trkindex = track_index;
01983 LoadLargestTrackFromEvent(i,track_index);
01984 int shower_index = -1;
01985
01986 shwindex = shower_index;
01987
01988 LoadLargestShowerFromEvent(i,shower_index);
01989
01990
01991
01992
01993
01994
01995 if(abs(inu)==14&&iaction==1){
01996
01997
01998 if(process!=1001) {
01999 NEvent_CC->Fill(eventSummary->nevent);
02000
02001
02002
02003 NShower_CC->Fill(ntpEvent->nshower);
02004 NTrack_CC->Fill(ntpEvent->ntrack);
02005 NShower_Vs_NTrack_CC->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02006 NHitPlanes_CC->Fill(ntpEvent->plane.n);
02007 NTrack_Vs_EvSC_CC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
02008 if(PassTrackCuts(0)) {
02009 NHitTrkFrac_CC->Fill(float(ntpTrack->nstrip)/float(ntpEvent->nstrip));
02010
02011 if(ntpEvent->ph.sigcor>0) ETrkFrac_CC->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
02012
02013 if((reco_emu+reco_eshw)>0) TrkMomFrac_CC->Fill(reco_emu/(reco_emu+reco_eshw));
02014
02015
02016
02017 if(reco_dircosz>0.) dsdz_CC->Fill(1./reco_dircosz);
02018
02019
02020
02021
02022
02023
02024
02025 TrkLen_CC->Fill(ntpTrack->ds);
02026 dedx_CC->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
02027 Float_t *CCNCVars = this->CCNCSepVars();
02028 HalfRatio_CC->Fill(CCNCVars[2]);
02029 delete CCNCVars;
02030 if(this->IsFidAll(0)){
02031 RangeCurvDiff_Vs_TrkLen_CC->
02032 Fill(ntpTrack->ds,
02033 2.*(this->RecoMuEnergy(1,0)-this->RecoMuEnergy(2,0))
02034 /(this->RecoMuEnergy(1,0)+this->RecoMuEnergy(2,0)));
02035 }
02036 }
02037 VtxShwEnergy_CC->Fill(reco_eshw);
02038
02039 }
02040 else {
02041 NEvent_QE->Fill(eventSummary->nevent);
02042
02043 NShower_QE->Fill(ntpEvent->nshower);
02044 NTrack_QE->Fill(ntpEvent->ntrack);
02045 NShower_Vs_NTrack_QE->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02046 NHitPlanes_QE->Fill(ntpEvent->plane.n);
02047 NTrack_Vs_EvSC_QE->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
02048 if(PassTrackCuts(0)) {
02049 NHitTrkFrac_QE->Fill(float(ntpTrack->nstrip)
02050 /float(ntpEvent->nstrip));
02051 ETrkFrac_QE->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
02052 TrkMomFrac_QE->Fill(reco_emu/(reco_emu+reco_eshw));
02053 TrkMomFrac0_QE->Fill(reco_emu/(reco_emu+reco_eshw));
02054 dsdz_QE->Fill(1./reco_dircosz);
02055 NHitTrkFrac_Vs_VtxShwEnergy_QE->Fill(reco_eshw,
02056 float(ntpTrack->nstrip)
02057 /float(ntpEvent->nstrip));
02058
02059 TrkLen_QE->Fill(ntpTrack->ds);
02060 dedx_QE->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
02061 Float_t *CCNCVars = this->CCNCSepVars();
02062 HalfRatio_QE->Fill(CCNCVars[2]);
02063 delete CCNCVars;
02064 if(this->IsFidAll(0)){
02065 RangeCurvDiff_Vs_TrkLen_QE->
02066 Fill(ntpTrack->ds,
02067 2.*(this->RecoMuEnergy(1,0)-this->RecoMuEnergy(2,0))
02068 /(this->RecoMuEnergy(1,0)+this->RecoMuEnergy(2,0)));
02069 }
02070 }
02071 else {
02072 VtxShwEnergy_QE->Fill(reco_eshw);
02073 }
02074
02075 }
02076 }
02077
02078
02079 if(iaction==0){
02080 NEvent_NC->Fill(eventSummary->nevent);
02081
02082 NShower_NC->Fill(ntpEvent->nshower);
02083 NTrack_NC->Fill(ntpEvent->ntrack);
02084 NShower_Vs_NTrack_NC->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02085 NHitPlanes_NC->Fill(ntpEvent->plane.n);
02086 NTrack_Vs_EvSC_NC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
02087 if(PassTrackCuts(0)) {
02088 NHitTrkFrac_NC->Fill(float(ntpTrack->nstrip)
02089 /float(ntpEvent->nstrip));
02090 ETrkFrac_NC->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
02091
02092
02093
02094 TrkMomFrac_NC->Fill(reco_emu/(reco_emu+reco_eshw));
02095 TrkMomFrac0_NC->Fill(reco_emu/(reco_emu+reco_eshw));
02096 dsdz_NC->Fill(1./reco_dircosz);
02097 NHitTrkFrac_Vs_VtxShwEnergy_NC->Fill(reco_eshw,
02098 float(ntpTrack->nstrip)
02099 /float(ntpEvent->nstrip));
02100 TrkLen_NC->Fill(ntpTrack->ds);
02101 dedx_NC->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
02102 Float_t *CCNCVars = this->CCNCSepVars();
02103 HalfRatio_NC->Fill(CCNCVars[2]);
02104 delete CCNCVars;
02105 if(this->IsFidAll(0)){
02106 RangeCurvDiff_Vs_TrkLen_NC->
02107 Fill(ntpTrack->ds,
02108 2.*(this->RecoMuEnergy(1,0)-this->RecoMuEnergy(2,0))
02109 /(this->RecoMuEnergy(1,0)+this->RecoMuEnergy(2,0)));
02110 }
02111 }
02112 VtxShwEnergy_NC->Fill(reco_eshw);
02113
02114 }
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128 }
02129
02130 }
02131
02132 }
02133 save.Write();
02134 save.Close();
02135
02136 }
02137
02138
02139
02140
02141 void MadEdAnalysis::MakeEff(std::string tag){
02142
02143 std::string savename = "EffHist_" + tag + ".root";
02144 TFile save(savename.c_str(),"RECREATE");
02145 save.cd();
02146
02147
02148
02149 TProfile *EvtEff = new TProfile("EvtEff","Event is reconstructed",200,0.,20.,0.,1.);
02150
02151 TProfile *MuEff = new TProfile("MuEff","Track is found and reconstructed",200,0.,20.,0.,1.);
02152
02153 TProfile *MuEffS = new TProfile("MuEffS","Single Track is reconstructed",200,0.,20.,0.,1.);
02154
02155
02156 TProfile *MuEffT = new TProfile("MuEffT","Good Track, fit.pass=1",200,0.,20.,0.,1.);
02157
02158
02159
02160 TProfile *MuEffA = new TProfile("MuEffA","All cuts",200,0.,20.,0.,1.);
02161
02162
02163 TProfile *MuEffFd = new TProfile("MuEffFd","Fuducial cut on vtx",200,0.,20.,0.,1.);
02164
02165 TProfile *MuEffQE = new TProfile("MuEffQE","QE Event",200,0.,20.,0.,1.);
02166
02167
02168
02169
02170 Int_t flavor;
02171
02172
02173
02174
02175 Float_t* mypid;
02176 Float_t Pidvalsk;
02177 Float_t Pidval2;
02178
02179
02180 for(int i=0;i<Nentries;i++){
02181
02182 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
02183 << "% done" << std::endl;
02184 this->GetEntry(i);
02185
02186
02187 if(eventSummary->nevent!=1) continue;
02188
02189
02190
02191 for(int i=0;i<eventSummary->nevent;i++){
02192
02193
02194 if(!LoadEvent(i)) continue;
02195
02196
02197
02198 int mc_entry = 0;
02199 if(!LoadTruthForRecoTH(i,mc_entry)) continue;
02200
02201
02202
02203
02204 int inu = INu(mc_entry);
02205 int iaction = IAction(mc_entry);
02206 int IRes=this->IResonance(mc_entry);
02207
02208 Float_t NuEn=this->TrueNuEnergy(mc_entry);
02209
02210 if(abs(inu)==14&&iaction==1) flavor=1;
02211 if(iaction==0) flavor =0;
02212
02213
02214
02215
02216
02217
02218
02219
02220 if(abs(inu)==14&&iaction==1){
02221
02222 Int_t evt=0;
02223
02224 if( LoadEvent(i)) { evt=1; EvtEff->Fill(this->TrueNuEnergy(),evt);}
02225 else {EvtEff->Fill(NuEn,evt);}
02226
02227
02228 if(IRes==1001){
02229 Int_t truQE=0;
02230
02231 if( LoadEvent(i)) {
02232 truQE=1;
02233 MuEffQE->Fill(NuEn,truQE);
02234 }
02235 else { MuEffQE->Fill(NuEn,truQE);}
02236
02237 }
02238
02239
02240 }
02241
02242
02243
02244
02245
02246
02247
02248 mypid=this->LikeliPID();
02249
02250 Pidvalsk=mypid[0];
02251 Pidval2=mypid[1];
02252
02253
02254
02255
02256
02257 int *tracks = ntpEvent->trk;
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268 if(abs(inu)==14&&iaction==1) {
02269
02270
02271
02272
02273 Int_t trk=0;
02274 if(ntpEvent->ntrack>0.) {
02275 Int_t trk=1;
02276
02277 MuEff->Fill(NuEn,trk);
02278
02279
02280
02281 Int_t trks=0;
02282 if (ntpEvent->ntrack==1){
02283
02284 trks=1;
02285 for(int k=0;k<ntpEvent->ntrack;k++){
02286 int index = tracks[k];
02287 MuEffS->Fill(NuEn,trks);
02288
02289
02290
02291
02292
02293 Int_t muon_countrk=0;
02294 if(PassCuts(index)) {
02295
02296 muon_countrk=1;
02297
02298 MuEffT->Fill(NuEn,muon_countrk);
02299
02300 }
02301
02302 else {
02303
02304 MuEffT->Fill(this->TrueNuEnergy(),muon_countrk);
02305
02306 }
02307
02308
02309 Int_t fud=0;
02310 if(PassCuts(index)) {
02311 if(ntpEvent->vtx.z<0.5 || ntpEvent->end.z>6.5
02312 || sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3))+(ntpEvent->vtx.y*ntpEvent->vtx.y))>1) {
02313 fud=0;
02314 MuEffFd->Fill(NuEn,fud);}
02315 else { fud=1;
02316 MuEffFd->Fill(NuEn,fud);}
02317 }
02318
02319
02320
02321
02322
02323
02324
02325
02326 }
02327 }
02328 else { MuEffS->Fill(NuEn,trks);}
02329
02330
02331 }
02332 else {
02333 MuEff->Fill(NuEn,trk);}
02334
02335
02336
02337
02338
02339
02340 }
02341
02342
02343
02344
02345
02346
02347
02348 Int_t ndcut=0;
02349 Int_t fdcut=0;
02350
02351 if(ntpEvent->vtx.z<0.5||ntpEvent->vtx.z>29.4
02352 ||(ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5)
02353 ||fabs(ntpEvent->vtx.x)>3.5||fabs(ntpEvent->vtx.y)>3.5) fdcut=1;
02354
02355
02356
02357 if(ntpEvent->vtx.z>0.5&&ntpEvent->end.z<6.5
02358 &&sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3))+(ntpEvent->vtx.y*ntpEvent->vtx.y))<1) ndcut=1;
02359
02360
02361
02362
02363
02364
02365 if(ntpEvent->ntrack!=1) continue;
02366
02367
02368 for(int j=0;j<ntpEvent->ntrack;j++){
02369 int index = tracks[j];
02370
02371
02372
02373
02374
02375 if(abs(inu)==14&&iaction==1) {
02376
02377
02378
02379
02380
02381
02382
02383 Int_t muon_countall=0;
02384 if(ntpEvent->vtx.z>0.5&&ntpEvent->end.z<6.5
02385 &&sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3))+(ntpEvent->vtx.y*ntpEvent->vtx.y))<1&&PassCuts(index)&&ntpEvent->ntrack==1) {
02386
02387 muon_countall=1;
02388
02389 MuEffA->Fill(NuEn,muon_countall);
02390 }
02391
02392
02393 else { MuEffA->Fill(NuEn,muon_countall);}
02394
02395
02396
02397
02398
02399
02400 }
02401
02402
02403 }
02404
02405
02406
02407 }
02408
02409
02410
02411
02412 }
02413
02414 save.Write();
02415 save.Close();
02416 }
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427 void MadEdAnalysis::DataHist(std::string tag){
02428
02429
02430 Float_t *PdQE;
02431 Float_t PidQECC;
02432 Float_t PidQENC;
02433 Float_t ProbNC;
02434 Float_t ProbCC;
02435 Float_t ProbQE;
02436
02437 Int_t PassTCut;
02438 Int_t ntrack;
02439 Int_t EvtSum;
02440
02441
02442 Float_t evlength;
02443 Float_t phfrac;
02444 Float_t phplane;
02445 Float_t RecoQEQ2;
02446 Float_t RecoMuAngN;
02447 Float_t QENuEn;
02448 Float_t ETrkFracN;
02449 Float_t nnval;
02450 Float_t cosz;
02451 Float_t lcosz;
02452 Float_t RecoY;
02453 Float_t RecoNuEn;
02454 Float_t RecoMuEn;
02455 Float_t RecoShwEn;
02456 Float_t HitFrac;
02457 Float_t EvtLen;
02458
02459
02460 std::string savename = "DataHist_" + tag + ".root";
02461
02462 TFile *save = new TFile(savename.c_str(),"RECREATE");
02463
02464
02465
02466 TH1F *NEvent = new TH1F("NEvent",
02467 "Number of Reco'd Events",
02468 20,0,20);
02469 NEvent->SetXTitle("Number of Events");
02470
02471
02472
02473 TH1F *NShower = new TH1F("NShower",
02474 "Number of Reco'd Showers ",
02475 20,0,20);
02476 NShower->SetXTitle("Number of Showers");
02477
02478
02479 TH1F *NTrack = new TH1F("NTrack",
02480 "Number of Reco'd Tracks",
02481 20,0,20);
02482 NTrack->SetXTitle("Number of Tracks");
02483
02484
02485 TH2F *NShower_Vs_NTrack =
02486 new TH2F("NShower_Vs_NTrack",
02487 "Number of Reco'd Showers Vs Tracks",
02488 20,0,5,20,0,5);
02489
02490 NShower_Vs_NTrack->SetXTitle("Number of Tracks");
02491 NShower_Vs_NTrack->SetYTitle("Number of Showers");
02492
02493
02494 TH1F *NHitTrkFrac=new TH1F("NHitTrkFrac",
02495 "Fraction of Event Hits in Primary Track ",
02496 100,0,1);
02497 NHitTrkFrac->SetXTitle("Hit Fraction");
02498
02499 TH1F *ETrkFrac=new TH1F("ETrkFrac",
02500 "Fraction of Event SigCor in Primary Track",
02501 100,0,1);
02502 ETrkFrac->SetXTitle("SigCor Fraction");
02503
02504
02505 TH1F *TrkMomFrac = new TH1F("TrkMomFrac",
02506 "Event Momentum Fraction in Track",
02507 100,0,1);
02508 TrkMomFrac->SetXTitle("Momentum (GeV)");
02509
02510
02511
02512 TH1F *VtxShwEnergy = new TH1F("VtxShwEnergy",
02513 "Vertex Shower Energy",
02514 1000,0,100);
02515 VtxShwEnergy->SetXTitle("Energy (GeV)");
02516
02517 TH1F *dsdz = new TH1F("dsdz","Primary Track dsdz",100,-1,20);
02518 dsdz->SetXTitle("dsdz");
02519
02520 TH1F *NHitPlanes = new TH1F("NHitPlanes","Number of Hit Planes",
02521 500,-0.5,499.5);
02522 NHitPlanes->SetXTitle("Number of Planes");
02523
02524
02525 TH2F *NHitTrkFrac_Vs_VtxShwEnergy = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy","Fraction of Event Hits in Primary Track Vs Vertex Shower Energy",1000,0,50,100,0,1);
02526 NHitTrkFrac_Vs_VtxShwEnergy->SetXTitle("Energy (GeV)");
02527 NHitTrkFrac_Vs_VtxShwEnergy->SetYTitle("Hit Fraction");
02528
02529 TH2F *NTrack_Vs_EvSC = new TH2F("NTrack_Vs_EvSC","Number of Reco'd Tracks vs Event SigCor",1000,0,100000,20,0,20);
02530 NTrack_Vs_EvSC->SetXTitle("Event SigCor");
02531 NTrack_Vs_EvSC->SetYTitle("Number of Tracks");
02532
02533
02534
02535
02536 TH1F *TrkLen = new TH1F("TrkLen","Track length",50,0,50);
02537
02538
02539 TH1F *dedx = new TH1F("dedx","Average dEdx",1000,0,500);
02540
02541
02542 TH1F *HalfRatio = new TH1F("HalfRatio","Charge Ratio: First Half/Second Half of Track",100,0,2);
02543
02544
02545
02546 TH2F *RangeCurvDiff_Vs_TrkLen = new TH2F("RangeCurvDiff_Vs_TrkLen","Diff in Momentum from Range and Curv vs Track Length",100,0,30,200,-3,17);
02547
02548 TTree *tree = new TTree("mytree","mytree");
02549
02550
02551 tree->Branch("evlength",&evlength,"evlength/F",32000);
02552
02553 tree->Branch("phfrac",&phfrac,"phfrac/F",32000);
02554
02555 tree->Branch("phplane",&phplane,"phplane/F",32000);
02556
02557
02558 tree->Branch("PidQECC",&PidQECC,"PidQECC/F",32000);
02559
02560 tree->Branch("PidQENC",&PidQECC,"PidQENC/F",32000);
02561
02562 tree->Branch("ProbCC",&ProbCC,"ProbCC/F",32000);
02563
02564 tree->Branch("ProbNC",&ProbCC,"ProbNC/F",32000);
02565
02566 tree->Branch("ProbQE",&ProbQE,"ProbQE/F",32000);
02567
02568 tree->Branch("nnval",&nnval,"nnval/F",32000);
02569 tree->Branch("RecoQEQ2",&RecoQEQ2,"RecoQEQ2/F",32000);
02570 tree->Branch("RecoMuAngN",&RecoMuAngN,"RecoMuAngN/F",32000);
02571
02572 tree->Branch("QENuEn",&QENuEn,"QENuEn/F",32000);
02573
02574 tree->Branch("PassTCut",&PassTCut,"PassTCut/I",32000);
02575
02576 tree->Branch("RecoY",&RecoY,"RecoY/F",32000);
02577
02578 tree->Branch("RecoNuEn",&RecoNuEn,"RecoNuEn/F",32000);
02579 tree->Branch("RecoMuEn",&RecoMuEn,"RecoMuEn/F",32000);
02580 tree->Branch("RecoShwEn",&RecoShwEn,"RecoShwEn/F",32000);
02581
02582 tree->Branch("HitFrac",&HitFrac,"HitFrac/F",32000);
02583 tree->Branch("EvtLen",&EvtLen,"EvtLen/F",32000);
02584 tree->Branch("cosz",&cosz,"cosz/F",32000);
02585 tree->Branch("lcosz",&lcosz,"lcosz/F",32000);
02586 tree->Branch("EvtSum",&EvtSum,"EvtSum/I",32000);
02587 tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
02588
02589
02590
02591
02592 TFile *likeliFile = new TFile("PidQE_R1.9n.root","READ");
02593 TH1F *myhist[9];
02594
02595
02596 if(likeliFile) {
02597
02598 cout << "PidQE_R1.9n.root found" << endl;
02599 cout << "Analysing MC" << endl;
02600
02601 myhist[0] = (TH1F*) likeliFile->Get("NHitTrkFrac_NC");
02602 myhist[1] = (TH1F*) likeliFile->Get("dedx_NC");
02603 myhist[2] = (TH1F*) likeliFile->Get("HalfRatio_NC");
02604 myhist[3] = (TH1F*) likeliFile->Get("NHitTrkFrac_CC");
02605 myhist[4] = (TH1F*) likeliFile->Get("dedx_CC");
02606 myhist[5] = (TH1F*) likeliFile->Get("HalfRatio_CC");
02607 myhist[6] = (TH1F*) likeliFile->Get("NHitTrkFrac_QE");
02608 myhist[7] = (TH1F*) likeliFile->Get("dedx_QE");
02609 myhist[8] = (TH1F*) likeliFile->Get("HalfRatio_QE");
02610 for(int i=0;i<9;i++){
02611 float integ = myhist[i]->Integral(); myhist[i]->Scale(1./integ);
02612 }
02613
02614
02615
02616 }
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626 for(int i=0;i<Nentries;i++){
02627
02628
02629
02630 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
02631 << "% done" << std::endl;
02632 if(!this->GetEntry(i)) continue;
02633
02634
02635
02636
02637
02638 if(eventSummary->nevent==0) continue;
02639
02640
02641 for(int i=0;i<eventSummary->nevent;i++){
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657 if(ntpEvent->vtx.z<0.5||ntpEvent->end.z>6.5
02658 ||sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3))+(ntpEvent->vtx.y*ntpEvent->vtx.y))>1) continue;
02659
02660
02661 if(!LoadEvent(i)) continue;
02662
02663
02664
02665
02666
02667 int *tracks = ntpEvent->trk;
02668
02669 EvtSum=eventSummary->nevent;
02670
02671 for(int j=0;j<ntpEvent->ntrack;j++){
02672
02673 int index = tracks[j];
02674
02675
02676
02677 NEvent->Fill(eventSummary->nevent);
02678
02679 NShower->Fill(ntpEvent->nshower);
02680 NTrack->Fill(ntpEvent->ntrack);
02681 NShower_Vs_NTrack->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02682 NHitPlanes->Fill(ntpEvent->plane.n);
02683 NTrack_Vs_EvSC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
02684
02685 if(PassTrackCuts(index)) {
02686 NHitTrkFrac->Fill(HitF(index));
02687 ETrkFrac->Fill(ETrkF(index));
02688 float trkenergy = this->RecoMuEnergy(0,index);
02689 float shwenergy = this->RecoShwEnergy(index);
02690
02691 TrkMomFrac->Fill(trkenergy/(trkenergy+shwenergy));
02692
02693 dsdz->Fill(1./ntpTrack->vtx.dcosz);
02694 NHitTrkFrac_Vs_VtxShwEnergy->Fill(this->RecoShwEnergy(index),HitF(index));
02695 TrkLen->Fill(ntpTrack->ds);
02696 dedx->Fill(ETrkF(index));
02697 Float_t *CCNCVars = this->CCNCSepVars();
02698 HalfRatio->Fill(CCNCVars[2]);
02699 delete CCNCVars;
02700 if(this->IsFidAll(0)){
02701 RangeCurvDiff_Vs_TrkLen->
02702 Fill(ntpTrack->ds,
02703 2.*(this->RecoMuEnergy(1,index)-this->RecoMuEnergy(2,index))
02704 /(this->RecoMuEnergy(1,index)+this->RecoMuEnergy(2,index)));
02705 }
02706 }
02707 VtxShwEnergy->Fill(this->RecoShwEnergy(index));
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723 if(PassTrackCuts(index)) PassTCut=1;
02724 else PassTCut=0;
02725
02726
02727
02728
02729 if(PassTrackCuts(index)) {
02730
02731
02732 ntrack = ntpEvent->ntrack;
02733
02734
02735
02736
02737 nnval=0;
02738 nnval=NeuNetEval(0);
02739
02740
02741
02742
02743
02744
02745 PdQE=this->MyLikeliQE(myhist);
02746
02747
02748 PidQECC=PdQE[0];
02749 PidQENC=PdQE[1];
02750 ProbNC=PdQE[2];
02751 ProbCC=PdQE[3];
02752 ProbQE=PdQE[4];
02753
02754
02755 Float_t evtphsig=0;
02756 Float_t trkplanes=0;
02757
02758 evtphsig=this->Evtphsig();
02759 trkplanes=this->Trkplanes(index);
02760
02761 if(this->Evtphsig()==0) evtphsig=0.000001;
02762 if(this->Trkplanes(index)==0) trkplanes=0.00001;
02763
02764 evlength=this->EvtLength();
02765 phfrac=this->Trkphsig(index)/evtphsig;
02766 phplane=this->Trkphsig(index)/trkplanes;
02767
02768
02769
02770 HitFrac=HitF(index);
02771
02772 EvtLen=ntpEvent->plane.n;
02773 ETrkFracN=ETrkF(index);
02774
02775 cosz = ntpTrack->vtx.dcosz;
02776 lcosz= ntpTrack->lin.dcosz;
02777
02778
02779
02780
02781 RecoNuEn = (RecoShwEnergy(index)/1.23) + RecoMuEnergy(0,index);
02782
02783 RecoMuEn = RecoMuEnergy(0,index);
02784 RecoShwEn = (RecoShwEnergy(index)/1.23);
02785
02786 QENuEn = RecoQENuEnergy(index);
02787
02788 RecoNuEn = (RecoShwEnergy(index)/1.23) + RecoMuEnergy(0,index);
02789
02790 RecoY = (RecoShwEnergy(index)/1.23)/RecoNuEn;
02791
02792
02793 RecoQEQ2=this->RecoQEQ2(index);
02794 RecoMuAngN=this->RecoMuDCosNeu(index);
02795
02796
02797
02798 tree->Fill();
02799
02800 }
02801
02802
02803
02804
02805
02806
02807
02808 }
02809
02810
02811 }
02812
02813
02814 }
02815
02816 tree->Write();
02817 save->Write();
02818 save->Close();
02819
02820
02821 }
02822
02823
02824 void MadEdAnalysis::MCHist(std::string tag){
02825
02826
02827
02828
02829 Int_t IAction;
02830 Int_t INu;
02831 Int_t flavor;
02832
02833 Float_t Y;
02834 Float_t W2;
02835 Float_t Q2;
02836 Int_t IRes;
02837 Float_t NuEn;
02838 Float_t TrueMuAngN;
02839 Float_t MuEn;
02840 Float_t ShwEn;
02841
02842
02843
02844
02845 Int_t PassTCut;
02846 Int_t ntrack;
02847 Int_t EvtSum;
02848
02849
02850 Float_t *PdQE;
02851 Float_t PidQECC;
02852 Float_t PidQENC;
02853 Float_t ProbNC;
02854 Float_t ProbCC;
02855 Float_t ProbQE;
02856 Float_t evlength;
02857 Float_t phfrac;
02858 Float_t phplane;
02859 Float_t RecoQEQ2;
02860 Float_t QEAngDiff;
02861 Float_t cosz;
02862 Float_t lcosz;
02863 Float_t RecoY;
02864 Float_t RecoNuEn;
02865 Float_t RecoMuEn;
02866 Float_t RecoShwEn;
02867 Float_t RecoMuAngN;
02868
02869 Float_t QENuEn;
02870 Float_t RecoEnDiff;
02871 Float_t QEEnDiff;
02872 Float_t QEQ2Diff;
02873 Float_t YDiff;
02874 Float_t HitFrac;
02875 Float_t EvtLen;
02876 Float_t ETrkFracCC;
02877 Float_t nnval;
02878
02879
02880 std::string savename = "MCHist_" + tag + ".root";
02881
02882 TFile *save = new TFile(savename.c_str(),"RECREATE");
02883
02884
02885
02886
02887
02888
02889
02890 TH1F *NEvent_NC = new TH1F("NEvent_NC",
02891 "Number of Reco'd Events - NC",
02892 100,0,20);
02893 NEvent_NC->SetXTitle("Number of Events");
02894 TH1F *NEvent_CC = new TH1F("NEvent_CC",
02895 "Number of Reco'd Events - CC",
02896 100,0,20);
02897 NEvent_CC->SetXTitle("Number of Events");
02898 TH1F *NEvent_QE = new TH1F("NEvent_QE",
02899 "Number of Reco'd Events - QE",
02900 100,0,20);
02901 NEvent_QE->SetXTitle("Number of Events");
02902
02903
02904 TH1F *NShower_NC = new TH1F("NShower_NC",
02905 "Number of Reco'd Showers - NC",
02906 100,0,20);
02907 NShower_NC->SetXTitle("Number of Showers");
02908 TH1F *NShower_CC = new TH1F("NShower_CC",
02909 "Number of Reco'd Showers - CC",
02910 100,0,20);
02911 NShower_CC->SetXTitle("Number of Showers");
02912 TH1F *NShower_QE = new TH1F("NShower_QE",
02913 "Number of Reco'd Showers - QE",
02914 100,0,20);
02915 NShower_QE->SetXTitle("Number of Showers");
02916
02917
02918 TH1F *NTrack_NC = new TH1F("NTrack_NC",
02919 "Number of Reco'd Tracks - NC",
02920 100,0,20);
02921 NTrack_NC->SetXTitle("Number of Tracks");
02922 TH1F *NTrack_CC = new TH1F("NTrack_CC",
02923 "Number of Reco'd Tracks - CC",
02924 100,0,20);
02925 NTrack_CC->SetXTitle("Number of Tracks");
02926 TH1F *NTrack_QE = new TH1F("NTrack_QE",
02927 "Number of Reco'd Tracks - QE",
02928 100,0,20);
02929 NTrack_QE->SetXTitle("Number of Tracks");
02930
02931
02932
02933
02934 TH2F *NShower_Vs_NTrack =
02935 new TH2F("NShower_Vs_NTrack",
02936 "Number of Reco'd Showers Vs Tracks- PassTcut",
02937 20,0,5,20,0,5);
02938 NShower_Vs_NTrack->SetXTitle("Number of Tracks");
02939 NShower_Vs_NTrack->SetYTitle("Number of Showers");
02940
02941
02942 TH2F *NHitTrkFrac_Vs_VtxShwEnergy
02943 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy",
02944 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - CC",1000,0,50,100,0,1);
02945 NHitTrkFrac_Vs_VtxShwEnergy->SetXTitle("Energy (GeV)");
02946 NHitTrkFrac_Vs_VtxShwEnergy->SetYTitle("Hit Fraction");
02947
02948
02949 TH1F *ETrkFrac=new TH1F("ETrkFrac",
02950 "Fraction of Event SigCor in Primary Track - PassTCut",
02951 100,0,1);
02952 ETrkFrac->SetXTitle("SigCor Fraction");
02953
02954
02955
02956
02957 TH2F *NShower_Vs_NTrack_NC =
02958 new TH2F("NShower_Vs_NTrack_NC",
02959 "Number of Reco'd Showers Vs Tracks- NC",
02960 100,0,5,100,0,5);
02961 NShower_Vs_NTrack_NC->SetXTitle("Number of Tracks");
02962 NShower_Vs_NTrack_NC->SetYTitle("Number of Showers");
02963 TH2F *NShower_Vs_NTrack_CC =
02964 new TH2F("NShower_Vs_NTrack_CC",
02965 "Number of Reco'd Showers Vs Tracks- CC",
02966 100,0,5,100,0,5);
02967 NShower_Vs_NTrack_CC->SetXTitle("Number of Tracks");
02968 NShower_Vs_NTrack_CC->SetYTitle("Number of Showers");
02969 TH2F *NShower_Vs_NTrack_QE =
02970 new TH2F("NShower_Vs_NTrack_QE",
02971 "Number of Reco'd Showers Vs Tracks- QE",
02972 100,0,5,100,0,5);
02973 NShower_Vs_NTrack_QE->SetXTitle("Number of Tracks");
02974 NShower_Vs_NTrack_QE->SetYTitle("Number of Showers");
02975
02976
02977 TH1F *NHitTrkFrac_NC=new TH1F("NHitTrkFrac_NC",
02978 "Fraction of Event Hits in Primary Track - NC",
02979 100,0,1);
02980 NHitTrkFrac_NC->SetXTitle("Hit Fraction");
02981 TH1F *NHitTrkFrac_CC=new TH1F("NHitTrkFrac_CC",
02982 "Fraction of Event Hits in Primary Track - CC",
02983 100,0,1);
02984 NHitTrkFrac_CC->SetXTitle("Hit Fraction");
02985 TH1F *NHitTrkFrac_QE=new TH1F("NHitTrkFrac_QE",
02986 "Fraction of Event Hits in Primary Track - QE",
02987 100,0,1);
02988 NHitTrkFrac_QE->SetXTitle("Hit Fraction");
02989
02990
02991
02992
02993
02994
02995
02996 TH1F *ETrkFrac_NC=new TH1F("ETrkFrac_NC",
02997 "Fraction of Event SigCor in Primary Track - NC",
02998 100,0,1);
02999 ETrkFrac_NC->SetXTitle("SigCor Fraction");
03000 TH1F *ETrkFrac_CC=new TH1F("ETrkFrac_CC",
03001 "Fraction of Event SigCor in Primary Track - CC",
03002 100,0,1);
03003 ETrkFrac_CC->SetXTitle("SigCor Fraction");
03004 TH1F *ETrkFrac_QE=new TH1F("ETrkFrac_QE",
03005 "Fraction of Event SigCor in Primary Track - QE",
03006 100,0,1);
03007 ETrkFrac_QE->SetXTitle("SigCor Fraction");
03008
03009
03010 TH1F *TrkMomFrac_NC = new TH1F("TrkMomFrac_NC",
03011 "Event Momentum Fraction in Track - NC",
03012 100,0,1);
03013 TrkMomFrac_NC->SetXTitle("Momentum (GeV)");
03014 TH1F *TrkMomFrac_CC = new TH1F("TrkMomFrac_CC",
03015 "Event Momentum Fraction in Track - CC",
03016 100,0,1);
03017 TrkMomFrac_CC->SetXTitle("Momentum (GeV)");
03018 TH1F *TrkMomFrac_QE = new TH1F("TrkMomFrac_QE",
03019 "Event Momentum Fraction in Track - QE",
03020 100,0,1);
03021 TrkMomFrac_QE->SetXTitle("Momentum (GeV)");
03022
03023
03024
03025 TH1F *VtxShwEnergy_NC = new TH1F("VtxShwEnergy_NC",
03026 "Vertex Shower Energy - NC",
03027 1000,0,100);
03028 VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
03029 TH1F *VtxShwEnergy_CC = new TH1F("VtxShwEnergy_CC",
03030 "Vertex Shower Energy - CC",
03031 1000,0,100);
03032 VtxShwEnergy_CC->SetXTitle("Energy (GeV)");
03033 TH1F *VtxShwEnergy_QE = new TH1F("VtxShwEnergy_QE",
03034 "Vertex Shower Energy - QE",
03035 1000,0,100);
03036 VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
03037
03038
03039 TH1F *dsdz_NC = new TH1F("dsdz_NC","Primary Track dsdz - NC",100,-1,20);
03040 dsdz_NC->SetXTitle("dsdz");
03041 TH1F *dsdz_CC = new TH1F("dsdz_CC","Primary Track dsdz - CC",100,-1,20);
03042 dsdz_CC->SetXTitle("dsdz");
03043 TH1F *dsdz_QE = new TH1F("dsdz_QE","Primary Track dsdz - QE",100,-1,20);
03044 dsdz_QE->SetXTitle("dsdz");
03045
03046 TH1F *NHitPlanes_NC = new TH1F("NHitPlanes_NC","Number of Hit Planes - NC",
03047 500,-0.5,499.5);
03048 NHitPlanes_NC->SetXTitle("Number of Planes");
03049 TH1F *NHitPlanes_CC = new TH1F("NHitPlanes_CC","Number of Hit Planes - CC",
03050 500,-0.5,499.5);
03051 NHitPlanes_NC->SetXTitle("Number of Planes");
03052 TH1F *NHitPlanes_QE = new TH1F("NHitPlanes_QE","Number of Hit Planes - QE",
03053 500,-0.5,499.5);
03054 NHitPlanes_QE->SetXTitle("Number of Planes");
03055
03056 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_NC
03057 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_NC",
03058 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - NC",1000,0,50,100,0,1);
03059 NHitTrkFrac_Vs_VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
03060 NHitTrkFrac_Vs_VtxShwEnergy_NC->SetYTitle("Hit Fraction");
03061 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_CC
03062 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_CC",
03063 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - CC",1000,0,50,100,0,1);
03064 NHitTrkFrac_Vs_VtxShwEnergy_CC->SetXTitle("Energy (GeV)");
03065 NHitTrkFrac_Vs_VtxShwEnergy_CC->SetYTitle("Hit Fraction");
03066 TH2F *NHitTrkFrac_Vs_VtxShwEnergy_QE
03067 = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_QE",
03068 "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - QE",1000,0,50,100,0,1);
03069 NHitTrkFrac_Vs_VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
03070 NHitTrkFrac_Vs_VtxShwEnergy_QE->SetYTitle("Hit Fraction");
03071
03072
03073 TH2F *NTrack_Vs_EvSC_NC
03074 = new TH2F("NTrack_Vs_EvSC_NC",
03075 "Number of Reco'd Tracks vs Event SigCor - NC",
03076 1000,0,100000,20,0,20);
03077 NTrack_Vs_EvSC_NC->SetXTitle("Event SigCor");
03078 NTrack_Vs_EvSC_NC->SetYTitle("Number of Tracks");
03079
03080 TH2F *NTrack_Vs_EvSC_CC
03081 = new TH2F("NTrack_Vs_EvSC_CC",
03082 "Number of Reco'd Tracks vs Event SigCor - CC",
03083 1000,0,100000,20,0,20);
03084 NTrack_Vs_EvSC_CC->SetXTitle("Event SigCor");
03085 NTrack_Vs_EvSC_CC->SetYTitle("Number of Tracks");
03086
03087 TH2F *NTrack_Vs_EvSC_QE
03088 = new TH2F("NTrack_Vs_EvSC_QE",
03089 "Number of Reco'd Tracks vs Event SigCor - QE",
03090 1000,0,100000,20,0,20);
03091 NTrack_Vs_EvSC_QE->SetXTitle("Event SigCor");
03092 NTrack_Vs_EvSC_QE->SetYTitle("Number of Tracks");
03093
03094
03095 TH2F *NHitTrkFrac_Vs_Y_CC
03096 = new TH2F("NHitTrkFrac_Vs_Y_CC",
03097 "Fraction of Event Hits in Primary Track Vs Y - CC",
03098 100,0,1,100,0,1);
03099 NHitTrkFrac_Vs_Y_CC->SetXTitle("Y");
03100 NHitTrkFrac_Vs_Y_CC->SetYTitle("Hit Fraction");
03101
03102
03103 TH1F *TrkLen_NC = new TH1F("TrkLen_NC","Track length - NC",100,0,50);
03104 TH1F *TrkLen_CC = new TH1F("TrkLen_CC","Track length - CC",100,0,50);
03105 TH1F *TrkLen_QE = new TH1F("TrkLen_QE","Track length - CC",100,0,50);
03106
03107
03108 TH1F *dedx_NC = new TH1F("dedx_NC","Average dEdx - NC",100,0,500);
03109 TH1F *dedx_CC = new TH1F("dedx_CC","Average dEdx - CC",100,0,500);
03110 TH1F *dedx_QE = new TH1F("dedx_QE","Average dEdx - QE",100,0,500);
03111
03112
03113 TH1F *HalfRatio_NC
03114 = new TH1F("HalfRatio_NC",
03115 "Charge Ratio: First Half/Second Half of Track - NC",
03116 150,-1,14);
03117 TH1F *HalfRatio_CC
03118 = new TH1F("HalfRatio_CC",
03119 "Charge Ratio: First Half/Second Half of Track - CC",
03120 150,-1,14);
03121 TH1F *HalfRatio_QE
03122 = new TH1F("HalfRatio_QE",
03123 "Charge Ratio: First Half/Second Half of Track - QE",
03124 150,-1,14);
03125
03126
03127 TH2F *RangeCurvDiff_Vs_TrkLen_NC
03128 = new TH2F("RangeCurvDiff_Vs_TrkLen_NC",
03129 "Diff in Momentum from Range and Curv vs Track Length - NC",
03130 100,0,30,200,-3,17);
03131 TH2F *RangeCurvDiff_Vs_TrkLen_CC
03132 = new TH2F("RangeCurvDiff_Vs_TrkLen_CC",
03133 "Diff in Momentum from Range and Curv vs Track Length - CC",
03134 100,0,30,200,-3,17);
03135 TH2F *RangeCurvDiff_Vs_TrkLen_QE
03136 = new TH2F("RangeCurvDiff_Vs_TrkLen_QE",
03137 "Diff in Momentum from Range and Curv vs Track Length - QE",
03138 100,0,30,200,-3,17);
03139
03140
03141
03142
03143
03144 TTree *tree = new TTree("mytree","mytree");
03145
03146 tree->Branch("evlength",&evlength,"evlength/F",32000);
03147
03148 tree->Branch("phfrac",&phfrac,"phfrac/F",32000);
03149
03150 tree->Branch("phplane",&phplane,"phplane/F",32000);
03151 tree->Branch("flavor",&flavor,"flavor/I",32000);
03152
03153 tree->Branch("PidQECC",&PidQECC,"PidQECC/F",32000);
03154
03155 tree->Branch("PidQENC",&PidQECC,"PidQENC/F",32000);
03156
03157 tree->Branch("ProbCC",&ProbCC,"ProbCC/F",32000);
03158
03159 tree->Branch("ProbNC",&ProbCC,"ProbNC/F",32000);
03160
03161 tree->Branch("ProbQE",&ProbQE,"ProbQE/F",32000);
03162
03163 tree->Branch("nnval",&nnval,"nnval/F",32000);
03164 tree->Branch("RecoQEQ2",&RecoQEQ2,"RecoQEQ2/F",32000);
03165 tree->Branch("RecoMuAngN",&RecoMuAngN,"RecoMuAngN/F",32000);
03166
03167 tree->Branch("RecoNuEn",&RecoNuEn,"RecoNuEn/F",32000);
03168 tree->Branch("RecoMuEn",&RecoMuEn,"RecoMuEn/F",32000);
03169 tree->Branch("RecoShwEn",&RecoShwEn,"RecoShwEn/F",32000);
03170
03171 tree->Branch("RecoY",&RecoY,"RecoY/F",32000);
03172 tree->Branch("QENuEn",&QENuEn,"QENuEn/F",32000);
03173
03174 tree->Branch("QEAngDiff",&QEAngDiff,"QEAngDiff/F",32000);
03175 tree->Branch("QEQ2Diff",&QEQ2Diff,"QEQ2Diff/F",32000);
03176 tree->Branch("QEEnDiff",&QEEnDiff,"QEEnDiff/F",32000);
03177
03178 tree->Branch("HitFrac",&HitFrac,"HitFrac/F",32000);
03179 tree->Branch("EvtLen",&EvtLen,"EvtLen/F",32000);
03180
03181 tree->Branch("ETrkFracCC",&ETrkFracCC,"ETrkFracCC/F",32000);
03182 tree->Branch("PassTCut",&PassTCut,"PassTCut/I",32000);
03183 tree->Branch("cosz",&cosz,"cosz/F",32000);
03184 tree->Branch("lcosz",&lcosz,"lcosz/F",32000);
03185 tree->Branch("EvtSum",&EvtSum,"EvtSum/I",32000);
03186
03187 tree->Branch("YDiff",&YDiff,"YDiff/F",32000);
03188
03189 tree->Branch("RecoEnDiff",&RecoEnDiff,"RecoEnDiff/F",32000);
03190
03191
03192
03193
03194
03195 tree->Branch("IRes",&IRes,"IRes/I",32000);
03196
03197 tree->Branch("Y",&Y,"Y/F",32000);
03198
03199 tree->Branch("NuEn",&NuEn,"NuEn/F",32000);
03200
03201 tree->Branch("TrueMuAngN",&TrueMuAngN,"TrueMuAngN/F",32000);
03202
03203 tree->Branch("Q2",&Q2,"Q2/F",32000);
03204
03205 tree->Branch("W2",&W2,"W2/F",32000);
03206
03207 tree->Branch("MuEn",&MuEn,"MuEn/F",32000);
03208 tree->Branch("ShwEn",&ShwEn,"ShwEn/F",32000);
03209
03210
03211
03212 ETrkFrac_CC->SetXTitle("SigCor Fraction");
03213
03214
03215
03216 TFile *likeliFile = new TFile("PidQE_R1.9n.root","READ");
03217 TH1F *myhist[9];
03218
03219
03220 if(likeliFile) {
03221
03222 cout << "PidQE_R1.9n.root found" << endl;
03223 cout << "Analysing MC" << endl;
03224
03225 myhist[0] = (TH1F*) likeliFile->Get("NHitTrkFrac_NC");
03226 myhist[1] = (TH1F*) likeliFile->Get("dedx_NC");
03227 myhist[2] = (TH1F*) likeliFile->Get("HalfRatio_NC");
03228 myhist[3] = (TH1F*) likeliFile->Get("NHitTrkFrac_CC");
03229 myhist[4] = (TH1F*) likeliFile->Get("dedx_CC");
03230 myhist[5] = (TH1F*) likeliFile->Get("HalfRatio_CC");
03231 myhist[6] = (TH1F*) likeliFile->Get("NHitTrkFrac_QE");
03232 myhist[7] = (TH1F*) likeliFile->Get("dedx_QE");
03233 myhist[8] = (TH1F*) likeliFile->Get("HalfRatio_QE");
03234 for(int i=0;i<9;i++){
03235 float integ = myhist[i]->Integral(); myhist[i]->Scale(1./integ);
03236 }
03237
03238
03239
03240 }
03241
03242
03243
03244
03245
03246
03247 for(int i=0;i<Nentries;i++){
03248
03249
03250
03251 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
03252 << "% done" << std::endl;
03253 if(!this->GetEntry(i)) continue;
03254
03255
03256
03257
03258
03259
03260 if(eventSummary->nevent==0) continue;
03261
03262
03263 for(int i=0;i<eventSummary->nevent;i++){
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279 if(ntpEvent->vtx.z<0.5||ntpEvent->end.z>6.5
03280 ||sqrt(((ntpEvent->vtx.x-1.3)*(ntpEvent->vtx.x-1.3))+(ntpEvent->vtx.y*ntpEvent->vtx.y))>1) continue;
03281
03282
03283
03284
03285
03286
03287
03288
03289 if(!LoadEvent(i)) continue;
03290 int mc_entry = 0;
03291
03292
03293 if(!LoadTruthForRecoTH(i,mc_entry)) continue;
03294
03295
03296 int inu = this->INu(mc_entry);
03297 int iaction = this->IAction(mc_entry);
03298
03299
03300 Y = this->Y(mc_entry);
03301 W2 = this->W2(mc_entry);
03302 IRes = this->IResonance(mc_entry);
03303 Q2 = this->Q2(mc_entry);
03304 NuEn =this->TrueNuEnergy(mc_entry);
03305 TrueMuAngN=this->TrueMuDCosNeu(mc_entry);
03306 MuEn = this->TrueMuEnergy(mc_entry);
03307 ShwEn = this->TrueShwEnergy(mc_entry);
03308
03309
03310
03311 if(abs(inu)==14&&iaction==1) flavor=1;
03312 if(iaction==0) flavor =0;
03313
03314 if(ntpEvent->ntrack!=1) continue;
03315
03316 int *tracks = ntpEvent->trk;
03317
03318 EvtSum=eventSummary->nevent;
03319
03320 for(int j=0;j<ntpEvent->ntrack;j++){
03321 int index = tracks[j];
03322
03323
03324
03325
03326
03327
03328
03329 if(abs(inu)==14&&iaction==1){
03330
03331 if(IRes!=1001) {
03332 NEvent_CC->Fill(eventSummary->nevent);
03333
03334
03335
03336 NShower_CC->Fill(ntpEvent->nshower);
03337 NTrack_CC->Fill(ntpEvent->ntrack);
03338 NShower_Vs_NTrack_CC->Fill(ntpEvent->ntrack,ntpEvent->nshower);
03339 NHitPlanes_CC->Fill(ntpEvent->plane.n);
03340 NTrack_Vs_EvSC_CC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
03341
03342
03343 if(PassTrackCuts(index)) {
03344 NHitTrkFrac_CC->Fill(HitF(index));
03345 ETrkFrac_CC->Fill(ETrkF(index));
03346
03347
03348 float trkenergy = this->RecoMuEnergy(0,index);
03349 float shwenergy = this->RecoShwEnergy(index);
03350
03351 TrkMomFrac_CC->Fill(trkenergy/(trkenergy+shwenergy));
03352
03353 dsdz_CC->Fill(1./ntpTrack->vtx.dcosz);
03354 NHitTrkFrac_Vs_VtxShwEnergy_CC->Fill(this->RecoShwEnergy(index),HitF(index));
03355 NHitTrkFrac_Vs_Y_CC->Fill(ntpTruth->y,HitF(index));
03356
03357 TrkLen_CC->Fill(ntpTrack->ds);
03358 dedx_CC->Fill(DeDx(index));
03359 Float_t *CCNCVars = this->CCNCSepVars();
03360 HalfRatio_CC->Fill(CCNCVars[2]);
03361 delete CCNCVars;
03362 if(this->IsFidAll(index)){
03363 RangeCurvDiff_Vs_TrkLen_CC->Fill(ntpTrack->ds,2.*(this->RecoMuEnergy(1,index)-this->RecoMuEnergy(2,index))/(this->RecoMuEnergy(1,index)+this->RecoMuEnergy(2,index)));
03364 }
03365
03366
03367 }
03368 VtxShwEnergy_CC->Fill(this->RecoShwEnergy(index));
03369
03370
03371 }
03372
03373
03374
03375 else if (IRes==1001) {
03376
03377
03378
03379
03380 NEvent_QE->Fill(eventSummary->nevent);
03381
03382 NShower_QE->Fill(ntpEvent->nshower);
03383 NTrack_QE->Fill(ntpEvent->ntrack);
03384 NShower_Vs_NTrack_QE->Fill(ntpEvent->ntrack,ntpEvent->nshower);
03385 NHitPlanes_QE->Fill(ntpEvent->plane.n);
03386 NTrack_Vs_EvSC_QE->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
03387
03388 if(PassTrackCuts(index)) {
03389
03390 NHitTrkFrac_QE->Fill(HitF(index));
03391
03392
03393
03394 ETrkFrac_QE->Fill(ETrkF(index));
03395 float trkenergy = this->RecoMuEnergy(0,index);
03396
03397 float shwenergy = this->RecoShwEnergy(index);
03398 TrkMomFrac_QE->Fill(trkenergy/(trkenergy+shwenergy));
03399
03400 dsdz_QE->Fill(1./ntpTrack->vtx.dcosz);
03401 NHitTrkFrac_Vs_VtxShwEnergy_QE->Fill(this->RecoShwEnergy(index),HitF(index));
03402
03403 TrkLen_QE->Fill(ntpTrack->ds);
03404 dedx_QE->Fill(DeDx(index));
03405 Float_t *CCNCVars = this->CCNCSepVars();
03406 HalfRatio_QE->Fill(CCNCVars[2]);
03407 delete CCNCVars;
03408 if(this->IsFidAll(index)){
03409 RangeCurvDiff_Vs_TrkLen_QE->
03410 Fill(ntpTrack->ds,
03411 2.*(this->RecoMuEnergy(1,index)-this->RecoMuEnergy(2,index))
03412 /(this->RecoMuEnergy(1,index)+this->RecoMuEnergy(2,index)));
03413 }
03414 }
03415 else {
03416 VtxShwEnergy_QE->Fill(this->RecoShwEnergy(index));
03417 }
03418
03419
03420 }
03421
03422 }
03423
03424
03425 if(iaction==0){
03426 NEvent_NC->Fill(eventSummary->nevent);
03427
03428 NShower_NC->Fill(ntpEvent->nshower);
03429 NTrack_NC->Fill(ntpEvent->ntrack);
03430 NShower_Vs_NTrack_NC->Fill(ntpEvent->ntrack,ntpEvent->nshower);
03431 NHitPlanes_NC->Fill(ntpEvent->plane.n);
03432 NTrack_Vs_EvSC_NC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
03433
03434 if(PassTrackCuts(index)) {
03435 NHitTrkFrac_NC->Fill(HitF(index));
03436 ETrkFrac_NC->Fill(ETrkF(index));
03437 float trkenergy = this->RecoMuEnergy(0,index);
03438 float shwenergy = this->RecoShwEnergy(index);
03439
03440 TrkMomFrac_NC->Fill(trkenergy/(trkenergy+shwenergy));
03441
03442 dsdz_NC->Fill(1./ntpTrack->vtx.dcosz);
03443 NHitTrkFrac_Vs_VtxShwEnergy_NC->Fill(this->RecoShwEnergy(index),HitF(index));
03444 TrkLen_NC->Fill(ntpTrack->ds);
03445 dedx_NC->Fill(DeDx(index));
03446 Float_t *CCNCVars = this->CCNCSepVars();
03447 HalfRatio_NC->Fill(CCNCVars[2]);
03448 delete CCNCVars;
03449 if(this->IsFidAll(0)){
03450 RangeCurvDiff_Vs_TrkLen_NC->
03451 Fill(ntpTrack->ds,
03452 2.*(this->RecoMuEnergy(1,index)-this->RecoMuEnergy(2,index))
03453 /(this->RecoMuEnergy(1,0)+this->RecoMuEnergy(2,0)));
03454 }
03455 }
03456 VtxShwEnergy_NC->Fill(this->RecoShwEnergy(index));
03457
03458 }
03459
03460
03461
03462
03463
03464
03465
03466
03467
03468
03469
03470
03471
03472 if(PassTrackCuts(index)) PassTCut=1;
03473 else PassTCut=0;
03474
03475
03476
03477
03478
03479 if(PassTrackCuts(index)) {
03480
03481
03482
03483 INu=inu;
03484 IAction=iaction;
03485
03486
03487
03488
03489
03490
03491
03492
03493 nnval=0;
03494
03495
03496
03497 ntrack = ntpEvent->ntrack;
03498
03499
03500
03501 PdQE=this->MyLikeliQE(myhist);
03502
03503
03504 PidQECC=PdQE[0];
03505 PidQENC=PdQE[1];
03506 ProbNC=PdQE[2];
03507 ProbCC=PdQE[3];
03508 ProbQE=PdQE[4];
03509
03510
03511
03512
03513 Float_t evtphsig=0;
03514 Float_t trkplanes=0;
03515
03516 evtphsig=this->Evtphsig();
03517 trkplanes=this->Trkplanes(index);
03518
03519 if(this->Evtphsig()==0) evtphsig=0.000001;
03520 if(this->Trkplanes(index)==0) trkplanes=0.00001;
03521
03522 evlength=this->EvtLength();
03523 phfrac=this->Trkphsig(index)/evtphsig;
03524 phplane=this->Trkphsig(index)/trkplanes;
03525
03526
03527
03528 cosz = ntpTrack->vtx.dcosz;
03529 lcosz= ntpTrack->lin.dcosz;
03530
03531 QENuEn = RecoQENuEnergy(index);
03532
03533 RecoNuEn = (RecoShwEnergy(index)/1.23) + RecoMuEnergy(0,index);
03534
03535 RecoMuEn= RecoMuEnergy(0,index);
03536 RecoShwEn = (RecoShwEnergy(index)/1.23);
03537
03538 RecoY = (RecoShwEnergy(index)/1.23)/RecoNuEn;
03539 RecoEnDiff = (RecoNuEn - NuEn)/NuEn;
03540
03541
03542
03543
03544 HitFrac=HitF(index);
03545
03546 EvtLen=ntpEvent->plane.n;
03547 ETrkFracCC=ETrkF(index);
03548
03549 ETrkFrac->Fill(ETrkF(index));
03550 NShower_Vs_NTrack->Fill(ntpEvent->ntrack,ntpEvent->nshower);
03551
03552 NHitTrkFrac_Vs_VtxShwEnergy->Fill(this->RecoShwEnergy(index),HitF(index));
03553
03554
03555
03556
03557
03558
03559
03560 if(Y==0) Y=0.000001;
03561 YDiff = (RecoY - Y)/Y;
03562
03563
03564
03565
03566 RecoQEQ2=this->RecoQEQ2(index);
03567 RecoMuAngN=this->RecoMuDCosNeu(index);
03568
03569 QEEnDiff = (QENuEn - NuEn)/NuEn;
03570
03571 if(TrueMuAngN==0) TrueMuAngN=0.00001;
03572
03573 QEAngDiff = (RecoMuAngN-TrueMuAngN)/TrueMuAngN;
03574
03575 if(Q2==0) Q2=0.000001;
03576 QEQ2Diff = (RecoQEQ2-Q2)/Q2;
03577
03578
03579
03580 tree->Fill();
03581
03582
03583
03584
03585
03586
03587
03588
03589 }
03590
03591
03592 }
03593
03594
03595 }
03596
03597
03598 }
03599
03600 tree->Write();
03601 save->Write();
03602 save->Close();
03603 }
03604
03605
03606
03607
03608 #endif // #ifdef madedanalysis_cxx
