loon: MadCBSQEAnalysis.cxx Source File

00001 #ifndef madcbsqeanalysis_cxx
00002 #define madcbsqeanalysis_cxx
00003 #include <iostream>
00004 #include "Mad/MadCBSQEAnalysis.h"
00005 
00006 using namespace std;
00007 
00008 //********************************************************
00009 MadCBSQEAnalysis::MadCBSQEAnalysis(TChain *chainSR,TChain *chainMC,
00010                                    TChain *chainTH,TChain *chainEM)
00011 {
00012 
00013   if(!chainSR) {
00014     record = 0;
00015     strecord = 0;
00016     emrecord = 0;
00017     mcrecord = 0;
00018     threcord = 0;
00019     Clear();
00020     whichSource = -1;
00021     isMC = true;
00022     isTH = true;
00023     isEM = true;
00024     Nentries = -1;
00025     return;
00026   }
00027   
00028   InitChain(chainSR,chainMC,chainTH,chainEM);
00029   whichSource = 0;  
00030   fLikeliFile = NULL;
00031   for(int i=0;i<12;i++) fLikeliHist[i] = NULL;
00032   cbsqe_TrkFrac = 0;
00033 }
00034 
00035 //********************************************************
00036 MadCBSQEAnalysis::MadCBSQEAnalysis(JobC *j,string path,int entries)
00037 {
00038   record = 0;
00039   strecord = 0;
00040   emrecord = 0;
00041   mcrecord = 0;
00042   threcord = 0;
00043   Clear();
00044   isMC = true;
00045   isTH = true;
00046   isEM = true;
00047   Nentries = entries;
00048   whichSource = 1;
00049   jcPath = path;
00050   fJC = j;
00051   fChain = NULL;
00052   fLikeliFile = NULL;
00053   for(int i=0;i<12;i++) fLikeliHist[i] = NULL;
00054   cbsqe_TrkFrac = 0;
00055 }
00056 
00057 //********************************************************
00058 MadCBSQEAnalysis::~MadCBSQEAnalysis()
00059 {
00060   if(fLikeliFile) fLikeliFile->Close();
00061 }
00062 
00063 //*********************************************************
00064 Bool_t MadCBSQEAnalysis::PassTruthSignal(Int_t mcevent)
00065 {
00066   if(!LoadTruth(mcevent)) return false;
00067   if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
00068   return false;
00069 }
00070 
00071 //*********************************************************
00072 Bool_t MadCBSQEAnalysis::PassAnalysisCuts(Int_t event)
00073 {
00074   if(!LoadEvent(event)) return false;
00075   if(PassBasicCuts(event)){
00076     int track = -1;
00077     LoadLargestTrackFromEvent(event,track);
00078     if(!IsFidVtx(track)) return false;
00079     if(PID(event,0)<-0.1) return true;
00080   }
00081   return false;
00082 }
00083 
00084 //********************************************************
00085 Bool_t MadCBSQEAnalysis::PassBasicCuts(Int_t event){
00086   if(!LoadEvent(event)) return false;  
00087   if(ntpEvent->ntrack!=1) return false;
00088   Int_t track = -1;
00089   LoadLargestTrackFromEvent(event,track);  
00090   if(!PassTrackCuts(track)) return false;
00091   return true;
00092 }
00093 
00094 //********************************************************
00095 Float_t MadCBSQEAnalysis::PID(Int_t event,Int_t method){
00096   
00097   if(!LoadEvent(event)) return 0;
00098   if(method==0){
00099     return LikeliQE(event);
00100   }
00101   return 0;
00102 
00103 }
00104 
00105 //********************************************************
00106 Float_t MadCBSQEAnalysis::RecoAnalysisEnergy(Int_t event){
00107 
00108   if(!LoadEvent(event)) return 0;
00109   int largestTrack = -1;
00110   LoadLargestTrackFromEvent(event,largestTrack);
00111   float emu = RecoMuEnergy(0,largestTrack);
00112   float ehad = RecoShwEnergy(event);
00113   float ereco = emu + ehad;
00114   float qe_ereco = RecoQENuEnergy(largestTrack);
00115   ereco = qe_ereco;
00116   return ereco;
00117 
00118 }
00119 
00120 //********************************************************
00121 Float_t MadCBSQEAnalysis::GetWeight(Int_t event)
00122 {  
00123   if(event>=0) return 1.0;
00124   return 0;
00125 }
00126 
00127 //********************************************************
00128 Bool_t MadCBSQEAnalysis::AddUserBranches(TTree *tree)
00129 {
00130   if(!tree) return false;
00131   tree->Branch("cbsqe_TrkFrac",&cbsqe_TrkFrac,"cbsqe_TrkFrac/D",32000);
00132   return true;
00133 }
00134 
00135 //********************************************************
00136 void MadCBSQEAnalysis::CallUserFunctions(Int_t event)
00137 {
00138   cbsqe_TrkFrac = 0;  
00139 
00140   if(!LoadEvent(event)) return;
00141   Int_t track = -1;
00142   if(!LoadLargestTrackFromEvent(event,track)) return;
00143   cbsqe_TrkFrac = double(ntpTrack->nstrip)/double(ntpEvent->nstrip);
00144 
00145 }
00146 
00147 //********************************************************
00148 Float_t MadCBSQEAnalysis::LikeliQE(Int_t event){
00149   
00150   if(!LoadEvent(event)) return 0;
00151   Int_t track = -1;
00152   if(!LoadLargestTrackFromEvent(event,track)) return 0;
00153 
00154   Float_t ProbNC = 1.;
00155   Float_t ProbDIS = 1.;
00156   Float_t ProbRES = 1.;
00157   Float_t ProbQE = 1.;
00158   Float_t PidParQEDIS = 0;
00159   Float_t PidParQERES = 0;
00160   Float_t PidParQENC = 0;  
00161   
00162   Float_t *CCNCVars = CCNCSepVars(track);
00163 
00164   //NC:
00165   int bin1 = fLikeliHist[0]->FindBin(float(ntpTrack->nstrip)
00166                                      /float(ntpEvent->nstrip));
00167   ProbNC *= fLikeliHist[0]->GetBinContent(bin1);
00168   
00169   bin1 = fLikeliHist[1]->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00170   ProbNC *= fLikeliHist[1]->GetBinContent(bin1);
00171   
00172   bin1 = fLikeliHist[2]->FindBin(CCNCVars[2]);
00173   ProbNC *= fLikeliHist[2]->GetBinContent(bin1);
00174   
00175   //DIS:
00176   int bin2 = fLikeliHist[3]->FindBin(float(ntpTrack->nstrip)
00177                                      /float(ntpEvent->nstrip));
00178   ProbDIS *= fLikeliHist[3]->GetBinContent(bin2);
00179   
00180   bin2 = fLikeliHist[4]->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00181   ProbDIS *= fLikeliHist[4]->GetBinContent(bin2);
00182   
00183   bin2 = fLikeliHist[5]->FindBin(CCNCVars[2]);
00184   ProbDIS *= fLikeliHist[5]->GetBinContent(bin2);
00185 
00186   //RES:
00187   int bin3 = fLikeliHist[6]->FindBin(float(ntpTrack->nstrip)
00188                                      /float(ntpEvent->nstrip));
00189   ProbRES *= fLikeliHist[6]->GetBinContent(bin3);
00190   
00191   bin3 = fLikeliHist[7]->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00192   ProbRES *= fLikeliHist[7]->GetBinContent(bin3);
00193   
00194   bin3 = fLikeliHist[8]->FindBin(CCNCVars[2]);
00195   ProbRES *= fLikeliHist[8]->GetBinContent(bin3);
00196   
00197   //QE:
00198   int bin4 = fLikeliHist[9]->FindBin(float(ntpTrack->nstrip)
00199                                      /float(ntpEvent->nstrip));
00200   ProbQE *= fLikeliHist[9]->GetBinContent(bin4);
00201   
00202   bin4 = fLikeliHist[10]->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00203   ProbQE *= fLikeliHist[10]->GetBinContent(bin4);
00204   
00205   bin4 = fLikeliHist[11]->FindBin(CCNCVars[2]);
00206   ProbQE *= fLikeliHist[11]->GetBinContent(bin4);
00207   
00208   delete CCNCVars;
00209 
00210   //pid pars
00211   if(ProbQE!=0&&ProbDIS!=0) 
00212     PidParQEDIS = sqrt(-TMath::Log(ProbQE))
00213       -sqrt(-TMath::Log(ProbDIS));
00214   else if(ProbQE==0&&ProbDIS==0) PidParQEDIS=10;
00215   else if(ProbQE==0) PidParQEDIS = 10.-sqrt(-TMath::Log(ProbDIS));
00216   else if(ProbDIS==0) PidParQEDIS = sqrt(-TMath::Log(ProbQE))-10.;
00217 
00218   if(ProbQE!=0&&ProbRES!=0) 
00219     PidParQERES = sqrt(-TMath::Log(ProbQE))
00220       -sqrt(-TMath::Log(ProbRES));
00221   else if(ProbQE==0&&ProbRES==0) PidParQERES=10;
00222   else if(ProbQE==0) PidParQERES = 10.-sqrt(-TMath::Log(ProbRES));
00223   else if(ProbRES==0) PidParQERES = sqrt(-TMath::Log(ProbQE))-10.;
00224   
00225   if(ProbQE!=0&&ProbNC!=0) 
00226     PidParQENC = sqrt(-TMath::Log(ProbQE))
00227       -sqrt(-TMath::Log(ProbNC));
00228   else if(ProbQE==0&&ProbNC==0) PidParQENC=10;
00229   else if(ProbQE==0) PidParQENC = 10.-sqrt(-TMath::Log(ProbNC));
00230   else if(ProbNC==0) PidParQENC = sqrt(-TMath::Log(ProbQE))-10.;
00231   
00232   return PidParQEDIS;
00233 }
00234 
00235 //******************************************
00236 void MadCBSQEAnalysis::MakeQEFile(std::string tag){
00237 
00238   std::string savename = "QEHistos_" + tag + ".root";
00239   TFile save(savename.c_str(),"RECREATE"); 
00240   save.cd();
00241   
00242   //#declare lots of histos:
00243   TH1F *NEvent_NC = new TH1F("NEvent_NC",
00244                               "Number of Reco'd Events - NC",
00245                               20,0,20);
00246   NEvent_NC->SetXTitle("Number of Events");
00247   TH1F *NEvent_DIS = new TH1F("NEvent_DIS",
00248                               "Number of Reco'd Events - DIS",
00249                               20,0,20);
00250   NEvent_DIS->SetXTitle("Number of Events");
00251   TH1F *NEvent_RES = new TH1F("NEvent_RES",
00252                               "Number of Reco'd Events - RES",
00253                               20,0,20);
00254   NEvent_RES->SetXTitle("Number of Events");
00255   TH1F *NEvent_QE = new TH1F("NEvent_QE",
00256                               "Number of Reco'd Events - QE",
00257                               20,0,20);
00258   NEvent_QE->SetXTitle("Number of Events");
00259 
00260 
00261   TH1F *NShower_NC = new TH1F("NShower_NC",
00262                               "Number of Reco'd Showers - NC",
00263                               20,0,20);
00264   NShower_NC->SetXTitle("Number of Showers");
00265   TH1F *NShower_DIS = new TH1F("NShower_DIS",
00266                               "Number of Reco'd Showers - DIS",
00267                               20,0,20);
00268   NShower_DIS->SetXTitle("Number of Showers");
00269   TH1F *NShower_RES = new TH1F("NShower_RES",
00270                               "Number of Reco'd Showers - RES",
00271                               20,0,20);
00272   NShower_RES->SetXTitle("Number of Showers");
00273   TH1F *NShower_QE = new TH1F("NShower_QE",
00274                               "Number of Reco'd Showers - QE",
00275                               20,0,20);
00276   NShower_QE->SetXTitle("Number of Showers");
00277 
00278 
00279   TH1F *NTrack_NC = new TH1F("NTrack_NC",
00280                              "Number of Reco'd Tracks - NC",
00281                              20,0,20);
00282   NTrack_NC->SetXTitle("Number of Tracks");
00283   TH1F *NTrack_DIS = new TH1F("NTrack_DIS",
00284                              "Number of Reco'd Tracks - DIS",
00285                              20,0,20);
00286   NTrack_DIS->SetXTitle("Number of Tracks");
00287   TH1F *NTrack_RES = new TH1F("NTrack_RES",
00288                              "Number of Reco'd Tracks - RES",
00289                              20,0,20);
00290   NTrack_RES->SetXTitle("Number of Tracks");
00291   TH1F *NTrack_QE = new TH1F("NTrack_QE",
00292                              "Number of Reco'd Tracks - QE",
00293                              20,0,20);
00294   NTrack_QE->SetXTitle("Number of Tracks");
00295 
00296   
00297   TH1F *FracSlcStpInEvt_NC = new TH1F("FracSlcStpInEvt_NC",
00298                                       "Fraction of Slice Strips in Event - NC",
00299                                       50,0,1);
00300   FracSlcStpInEvt_NC->SetXTitle("Fraction of Strips");
00301   TH1F *FracSlcStpInEvt_DIS = new TH1F("FracSlcStpInEvt_DIS",
00302                                       "Fraction of Slice Strips in Event - DIS",
00303                                       50,0,1);
00304   FracSlcStpInEvt_DIS->SetXTitle("Fraction of Strips");
00305   TH1F *FracSlcStpInEvt_RES = new TH1F("FracSlcStpInEvt_RES",
00306                                       "Fraction of Slice Strips in Event - RES",
00307                                       50,0,1);
00308   FracSlcStpInEvt_RES->SetXTitle("Fraction of Strips");
00309   TH1F *FracSlcStpInEvt_QE = new TH1F("FracSlcStpInEvt_QE",
00310                                       "Fraction of Slice Strips in Event - QE",
00311                                       50,0,1);
00312   FracSlcStpInEvt_QE->SetXTitle("Fraction of Strips");
00313   
00314 
00315   TH2F *NShower_Vs_NTrack_NC = 
00316     new TH2F("NShower_Vs_NTrack_NC",
00317              "Number of Reco'd Showers Vs Tracks- NC",
00318              20,0,20,20,0,20);
00319   NShower_Vs_NTrack_NC->SetXTitle("Number of Tracks");
00320   NShower_Vs_NTrack_NC->SetYTitle("Number of Showers");
00321   TH2F *NShower_Vs_NTrack_DIS = 
00322     new TH2F("NShower_Vs_NTrack_DIS",
00323              "Number of Reco'd Showers Vs Tracks- DIS",
00324              20,0,20,20,0,20);
00325   NShower_Vs_NTrack_DIS->SetXTitle("Number of Tracks");
00326   NShower_Vs_NTrack_DIS->SetYTitle("Number of Showers");
00327   TH2F *NShower_Vs_NTrack_RES = 
00328     new TH2F("NShower_Vs_NTrack_RES",
00329              "Number of Reco'd Showers Vs Tracks- RES",
00330              20,0,20,20,0,20);
00331   NShower_Vs_NTrack_RES->SetXTitle("Number of Tracks");
00332   NShower_Vs_NTrack_RES->SetYTitle("Number of Showers");
00333   TH2F *NShower_Vs_NTrack_QE = 
00334     new TH2F("NShower_Vs_NTrack_QE",
00335              "Number of Reco'd Showers Vs Tracks- QE",
00336              20,0,20,20,0,20);
00337   NShower_Vs_NTrack_QE->SetXTitle("Number of Tracks");
00338   NShower_Vs_NTrack_QE->SetYTitle("Number of Showers");
00339 
00340 
00341   TH1F *NHitTrkFrac_NC=new TH1F("NHitTrkFrac_NC",
00342                                 "Fraction of Event Hits in Primary Track - NC",
00343                                 100,0,1);
00344   NHitTrkFrac_NC->SetXTitle("Hit Fraction");
00345   TH1F *NHitTrkFrac_DIS=new TH1F("NHitTrkFrac_DIS",
00346                                 "Fraction of Event Hits in Primary Track - DIS",
00347                                 100,0,1);
00348   NHitTrkFrac_DIS->SetXTitle("Hit Fraction");
00349   TH1F *NHitTrkFrac_RES=new TH1F("NHitTrkFrac_RES",
00350                                 "Fraction of Event Hits in Primary Track - RES",
00351                                 100,0,1);
00352   NHitTrkFrac_RES->SetXTitle("Hit Fraction");
00353   TH1F *NHitTrkFrac_QE=new TH1F("NHitTrkFrac_QE",
00354                                 "Fraction of Event Hits in Primary Track - QE",
00355                                 100,0,1);
00356   NHitTrkFrac_QE->SetXTitle("Hit Fraction");
00357 
00358 
00359   TH1F *ETrkFrac_NC=new TH1F("ETrkFrac_NC",
00360                              "Fraction of Event SigCor in Primary Track - NC",
00361                              100,0,1);
00362   ETrkFrac_NC->SetXTitle("SigCor Fraction");
00363   TH1F *ETrkFrac_DIS=new TH1F("ETrkFrac_DIS",
00364                              "Fraction of Event SigCor in Primary Track - DIS",
00365                              100,0,1);
00366   ETrkFrac_DIS->SetXTitle("SigCor Fraction");
00367   TH1F *ETrkFrac_RES=new TH1F("ETrkFrac_RES",
00368                              "Fraction of Event SigCor in Primary Track - RES",
00369                              100,0,1);
00370   ETrkFrac_RES->SetXTitle("SigCor Fraction");
00371   TH1F *ETrkFrac_QE=new TH1F("ETrkFrac_QE",
00372                              "Fraction of Event SigCor in Primary Track - QE",
00373                              100,0,1);
00374   ETrkFrac_QE->SetXTitle("SigCor Fraction");
00375 
00376 
00377   TH1F *TrkMomFrac_NC = new TH1F("TrkMomFrac_NC",
00378                                  "Event Momentum Fraction in Track - NC",
00379                                  100,0,1);
00380   TrkMomFrac_NC->SetXTitle("Momentum (GeV)");
00381   TH1F *TrkMomFrac_DIS = new TH1F("TrkMomFrac_DIS",
00382                                  "Event Momentum Fraction in Track - DIS",
00383                                  100,0,1);
00384   TrkMomFrac_DIS->SetXTitle("Momentum (GeV)");
00385   TH1F *TrkMomFrac_RES = new TH1F("TrkMomFrac_RES",
00386                                  "Event Momentum Fraction in Track - RES",
00387                                  100,0,1);
00388   TrkMomFrac_RES->SetXTitle("Momentum (GeV)");
00389   TH1F *TrkMomFrac_QE = new TH1F("TrkMomFrac_QE",
00390                                  "Event Momentum Fraction in Track - QE",
00391                                  100,0,1);
00392   TrkMomFrac_QE->SetXTitle("Momentum (GeV)");
00393 
00394 
00395   TH1F *TrkMomFrac0_NC = new TH1F("TrkMomFrac0_NC",
00396                                  "Event Momentum Fraction in Track - NC",
00397                                  100,0,1);
00398   TrkMomFrac0_NC->SetXTitle("Momentum (GeV)");
00399   TH1F *TrkMomFrac0_DIS = new TH1F("TrkMomFrac0_DIS",
00400                                  "Event Momentum Fraction in Track - DIS",
00401                                  100,0,1);
00402   TrkMomFrac0_DIS->SetXTitle("Momentum (GeV)");
00403   TH1F *TrkMomFrac0_RES = new TH1F("TrkMomFrac0_RES",
00404                                  "Event Momentum Fraction in Track - RES",
00405                                  100,0,1);
00406   TrkMomFrac0_RES->SetXTitle("Momentum (GeV)");
00407   TH1F *TrkMomFrac0_QE = new TH1F("TrkMomFrac0_QE",
00408                                  "Event Momentum Fraction in Track - QE",
00409                                  100,0,1);
00410   TrkMomFrac0_QE->SetXTitle("Momentum (GeV)");
00411 
00412 
00413   TH1F *VtxShwEnergy_NC = new TH1F("VtxShwEnergy_NC",
00414                                    "Vertex Shower Energy - NC",
00415                                    1000,0,100);
00416   VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
00417   TH1F *VtxShwEnergy_DIS = new TH1F("VtxShwEnergy_DIS",
00418                                    "Vertex Shower Energy - DIS",
00419                                    1000,0,100);
00420   VtxShwEnergy_DIS->SetXTitle("Energy (GeV)");
00421   TH1F *VtxShwEnergy_RES = new TH1F("VtxShwEnergy_RES",
00422                                    "Vertex Shower Energy - RES",
00423                                    1000,0,100);
00424   VtxShwEnergy_RES->SetXTitle("Energy (GeV)");
00425   TH1F *VtxShwEnergy_QE = new TH1F("VtxShwEnergy_QE",
00426                                    "Vertex Shower Energy - QE",
00427                                    1000,0,100);
00428   VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
00429   
00430 
00431   TH1F *dsdz_NC = new TH1F("dsdz_NC","Primary Track dsdz - NC",100,-1,20);
00432   dsdz_NC->SetXTitle("dsdz");
00433   TH1F *dsdz_DIS = new TH1F("dsdz_DIS","Primary Track dsdz - DIS",100,-1,20);
00434   dsdz_DIS->SetXTitle("dsdz");
00435   TH1F *dsdz_RES = new TH1F("dsdz_RES","Primary Track dsdz - RES",100,-1,20);
00436   dsdz_RES->SetXTitle("dsdz");
00437   TH1F *dsdz_QE = new TH1F("dsdz_QE","Primary Track dsdz - QE",100,-1,20);
00438   dsdz_QE->SetXTitle("dsdz");
00439   
00440   TH1F *NHitPlanes_NC = new TH1F("NHitPlanes_NC","Number of Hit Planes - NC",
00441                                  500,-0.5,499.5);
00442   NHitPlanes_NC->SetXTitle("Number of Planes");
00443   TH1F *NHitPlanes_DIS = new TH1F("NHitPlanes_DIS","Number of Hit Planes - DIS",
00444                                  500,-0.5,499.5);
00445   NHitPlanes_DIS->SetXTitle("Number of Planes");
00446   TH1F *NHitPlanes_RES = new TH1F("NHitPlanes_RES","Number of Hit Planes - RES",
00447                                  500,-0.5,499.5);
00448   NHitPlanes_RES->SetXTitle("Number of Planes");
00449   TH1F *NHitPlanes_QE = new TH1F("NHitPlanes_QE","Number of Hit Planes - QE",
00450                                  500,-0.5,499.5);
00451   NHitPlanes_QE->SetXTitle("Number of Planes");
00452 
00453   TH2F *NHitTrkFrac_Vs_VtxShwEnergy_NC 
00454     = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_NC",
00455                "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - NC",1000,0,100,100,0,1);
00456   NHitTrkFrac_Vs_VtxShwEnergy_NC->SetXTitle("Energy (GeV)");
00457   NHitTrkFrac_Vs_VtxShwEnergy_NC->SetYTitle("Hit Fraction");
00458   TH2F *NHitTrkFrac_Vs_VtxShwEnergy_DIS 
00459     = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_DIS",
00460                "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - DIS",1000,0,100,100,0,1);
00461   NHitTrkFrac_Vs_VtxShwEnergy_DIS->SetXTitle("Energy (GeV)");
00462   NHitTrkFrac_Vs_VtxShwEnergy_DIS->SetYTitle("Hit Fraction");
00463   TH2F *NHitTrkFrac_Vs_VtxShwEnergy_RES 
00464     = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_RES",
00465                "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - RES",1000,0,100,100,0,1);
00466   NHitTrkFrac_Vs_VtxShwEnergy_RES->SetXTitle("Energy (GeV)");
00467   NHitTrkFrac_Vs_VtxShwEnergy_RES->SetYTitle("Hit Fraction");
00468   TH2F *NHitTrkFrac_Vs_VtxShwEnergy_QE 
00469     = new TH2F("NHitTrkFrac_Vs_VtxShwEnergy_QE",
00470                "Fraction of Event Hits in Primary Track Vs Vertex Shower Energy - QE",1000,0,100,100,0,1);
00471   NHitTrkFrac_Vs_VtxShwEnergy_QE->SetXTitle("Energy (GeV)");
00472   NHitTrkFrac_Vs_VtxShwEnergy_QE->SetYTitle("Hit Fraction");
00473 
00474 
00475   TH2F *NTrack_Vs_EvSC_NC 
00476     = new TH2F("NTrack_Vs_EvSC_NC",
00477                "Number of Reco'd Tracks vs Event SigCor - NC",
00478                1000,0,100000,20,0,20);
00479   NTrack_Vs_EvSC_NC->SetXTitle("Event SigCor");
00480   NTrack_Vs_EvSC_NC->SetYTitle("Number of Tracks");
00481 
00482   TH2F *NTrack_Vs_EvSC_DIS 
00483     = new TH2F("NTrack_Vs_EvSC_DIS",
00484                "Number of Reco'd Tracks vs Event SigCor - DIS",
00485                1000,0,100000,20,0,20);
00486   NTrack_Vs_EvSC_DIS->SetXTitle("Event SigCor");
00487   NTrack_Vs_EvSC_DIS->SetYTitle("Number of Tracks");  
00488 
00489   TH2F *NTrack_Vs_EvSC_RES 
00490     = new TH2F("NTrack_Vs_EvSC_RES",
00491                "Number of Reco'd Tracks vs Event SigCor - RES",
00492                1000,0,100000,20,0,20);
00493   NTrack_Vs_EvSC_RES->SetXTitle("Event SigCor");
00494   NTrack_Vs_EvSC_RES->SetYTitle("Number of Tracks");  
00495 
00496   TH2F *NTrack_Vs_EvSC_QE 
00497     = new TH2F("NTrack_Vs_EvSC_QE",
00498                "Number of Reco'd Tracks vs Event SigCor - QE",
00499                1000,0,100000,20,0,20);
00500   NTrack_Vs_EvSC_QE->SetXTitle("Event SigCor");
00501   NTrack_Vs_EvSC_QE->SetYTitle("Number of Tracks");  
00502 
00503   TH2F *NHitTrkFrac_Vs_Y_RES 
00504     = new TH2F("NHitTrkFrac_Vs_Y_RES",
00505                "Fraction of Event Hits in Primary Track Vs Y - RES",
00506                100,0,1,100,0,1);
00507   NHitTrkFrac_Vs_Y_RES->SetXTitle("Y");
00508   NHitTrkFrac_Vs_Y_RES->SetYTitle("Hit Fraction");
00509 
00510   TH2F *NHitTrkFrac_Vs_Y_DIS 
00511     = new TH2F("NHitTrkFrac_Vs_Y_DIS",
00512                "Fraction of Event Hits in Primary Track Vs Y - DIS",
00513                100,0,1,100,0,1);
00514   NHitTrkFrac_Vs_Y_DIS->SetXTitle("Y");
00515   NHitTrkFrac_Vs_Y_DIS->SetYTitle("Hit Fraction");
00516 
00517   //old discrim histos:
00518   TH1F *TrkLen_NC = new TH1F("TrkLen_NC","Track length - NC",50,0,50);
00519   TH1F *TrkLen_DIS = new TH1F("TrkLen_DIS","Track length - DIS",50,0,50);
00520   TH1F *TrkLen_RES = new TH1F("TrkLen_RES","Track length - RES",50,0,50);
00521   TH1F *TrkLen_QE = new TH1F("TrkLen_QE","Track length - QE",50,0,50);
00522 
00523   TH1F *dedx_NC = new TH1F("dedx_NC","Average dEdx - NC",1000,0,2000);  
00524   TH1F *dedx_DIS = new TH1F("dedx_DIS","Average dEdx - DIS",1000,0,2000);
00525   TH1F *dedx_RES = new TH1F("dedx_RES","Average dEdx - RES",1000,0,2000);
00526   TH1F *dedx_QE = new TH1F("dedx_QE","Average dEdx - QE",1000,0,2000);
00527 
00528   
00529   TH1F *HalfRatio_NC 
00530     = new TH1F("HalfRatio_NC",
00531                "Charge Ratio: First Half/Second Half of Track - NC",
00532                150,-1,14);
00533   TH1F *HalfRatio_DIS 
00534     = new TH1F("HalfRatio_DIS",
00535                "Charge Ratio: First Half/Second Half of Track - DIS",
00536                150,-1,14);
00537   TH1F *HalfRatio_RES 
00538     = new TH1F("HalfRatio_RES",
00539                "Charge Ratio: First Half/Second Half of Track - RES",
00540                150,-1,14);
00541   TH1F *HalfRatio_QE 
00542     = new TH1F("HalfRatio_QE",
00543                "Charge Ratio: First Half/Second Half of Track - QE",
00544                150,-1,14);
00545   
00546 
00547   TH2F *RangeCurvDiff_Vs_TrkLen_NC 
00548     = new TH2F("RangeCurvDiff_Vs_TrkLen_NC",
00549                "Diff in Momentum from Range and Curv vs Track Length - NC",
00550                100,0,30,200,-3,17);
00551   TH2F *RangeCurvDiff_Vs_TrkLen_DIS 
00552     = new TH2F("RangeCurvDiff_Vs_TrkLen_DIS",
00553                "Diff in Momentum from Range and Curv vs Track Length - DIS",
00554                100,0,30,200,-3,17);
00555   TH2F *RangeCurvDiff_Vs_TrkLen_RES 
00556     = new TH2F("RangeCurvDiff_Vs_TrkLen_RES",
00557                "Diff in Momentum from Range and Curv vs Track Length - RES",
00558                100,0,30,200,-3,17);
00559   TH2F *RangeCurvDiff_Vs_TrkLen_QE 
00560     = new TH2F("RangeCurvDiff_Vs_TrkLen_QE",
00561                "Diff in Momentum from Range and Curv vs Track Length - QE",
00562                100,0,30,200,-3,17);
00563 
00564   for(int i=0;i<Nentries;i++){
00565     
00566     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00567                             << "% done" << std::endl;
00568     if(!this->GetEntry(i)) continue;
00569     
00570     //get which detector this snarl is from
00571     int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
00572     
00573     //Calculate number of events per slice for ND:
00574     Int_t *nEventsPerSlice = new Int_t[eventSummary->nslice];
00575     for(int i=0;i<eventSummary->nslice;i++) nEventsPerSlice[i]= 0;
00576     if(theDetector==0) {
00577       for(int i=0;i<eventSummary->nslice;i++){
00578         nEventsPerSlice[i]= 0;
00579         if(!LoadEvent(i)) continue;
00580         nEventsPerSlice[ntpEvent->slc] += 1;
00581       }
00582     }
00583 
00584     for(int j=0;j<eventSummary->nevent;j++){
00585       
00586       if(!LoadEvent(j)) continue;  //does event exist
00587 
00588       //check for vertex in fiducial region
00589       if(!IsFidVtxEvt(j)) continue;
00590      
00591       //fraction of hits in slice used in event:
00592       float nSliceHitInEventFrac = 0;
00593       if(LoadSlice(ntpEvent->slc)) {
00594         nSliceHitInEventFrac = float(ntpEvent->nstrip)/float(ntpSlice->nstrip);
00595       }
00596       else cout << "Couldn't load slice!" << endl;
00597       
00598       //if no good truth candidate, we loose event
00599       int mcevent = -1;
00600       if(LoadTruthForRecoTH(j,mcevent)){ 
00601         
00602         //define shower energy per event here
00603         //sum energy from all showers in event
00604         Int_t *showers = ntpEvent->shw;
00605         float shwenergy = 0;
00606         for(int k=0;k<ntpEvent->nshower;k++){
00607           shwenergy += this->RecoShwEnergy(showers[k]);
00608         }
00609         
00610         if(PassTruthSignal(mcevent)){   //numu CC events
00611           
00612           if(this->IResonance(mcevent)==1003) {  //DIS
00613             
00614             if(theDetector==0) NEvent_DIS->Fill(nEventsPerSlice[ntpEvent->slc]);
00615             else if(theDetector==1) NEvent_DIS->Fill(eventSummary->nevent);
00616             
00617             if((theDetector==0 && nEventsPerSlice[ntpEvent->slc]==1) || 
00618                (theDetector==1 && eventSummary->nevent==1)){
00619               
00620               NShower_DIS->Fill(ntpEvent->nshower);
00621               NTrack_DIS->Fill(ntpEvent->ntrack);
00622               NShower_Vs_NTrack_DIS->Fill(ntpEvent->ntrack,
00623                                          ntpEvent->nshower);
00624               NHitPlanes_DIS->Fill(ntpEvent->plane.n);
00625               NTrack_Vs_EvSC_DIS->Fill(ntpEvent->ph.sigcor,
00626                                       ntpEvent->ntrack);
00627               
00628               if(PassBasicCuts(j)) { //pass basic cuts loads largest track
00629                 
00630                 Float_t extraTrkStp = 0;
00631 
00632                 FracSlcStpInEvt_DIS->Fill(nSliceHitInEventFrac);
00633 
00634                 NHitTrkFrac_DIS->Fill(float(ntpTrack->nstrip+extraTrkStp)
00635                                      /float(ntpEvent->nstrip+extraTrkStp));
00636                 ETrkFrac_DIS->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
00637                 float trkenergy = this->RecoMuEnergy(0,ntpTrack->index);
00638                 TrkMomFrac_DIS->Fill(trkenergy/(trkenergy+shwenergy));
00639                 dsdz_DIS->Fill(1./ntpTrack->vtx.dcosz);
00640                 NHitTrkFrac_Vs_VtxShwEnergy_DIS
00641                   ->Fill(shwenergy,(float(ntpTrack->nstrip)+extraTrkStp)
00642                          /(float(ntpEvent->nstrip)+extraTrkStp));
00643                 
00644                 NHitTrkFrac_Vs_Y_DIS
00645                   ->Fill(ntpTruth->y,(float(ntpTrack->nstrip)+extraTrkStp)
00646                          /(float(ntpEvent->nstrip)+extraTrkStp));
00647                 
00648                 TrkLen_DIS->Fill(ntpTrack->ds);
00649                 dedx_DIS->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00650                 Float_t *CCNCVars = this->CCNCSepVars(ntpTrack->index);
00651                 HalfRatio_DIS->Fill(CCNCVars[2]);
00652                 delete CCNCVars;
00653                 if(this->IsFidAll(ntpTrack->index)){
00654                   RangeCurvDiff_Vs_TrkLen_DIS
00655                     ->Fill(ntpTrack->ds,
00656                            2.*(this->RecoMuEnergy(1,ntpTrack->index) - 
00657                                this->RecoMuEnergy(2,ntpTrack->index))
00658                            /(this->RecoMuEnergy(1,ntpTrack->index) + 
00659                              this->RecoMuEnergy(2,ntpTrack->index)));
00660                 }
00661               }
00662               VtxShwEnergy_DIS->Fill(shwenergy);
00663             }
00664           }
00665 
00666           else if(this->IResonance(mcevent)==1002) {  // RES
00667             
00668             if(theDetector==0) NEvent_RES->Fill(nEventsPerSlice[ntpEvent->slc]);
00669             else if(theDetector==1) NEvent_RES->Fill(eventSummary->nevent);
00670             
00671             if((theDetector==0 && nEventsPerSlice[ntpEvent->slc]==1) || 
00672                (theDetector==1 && eventSummary->nevent==1)){
00673               
00674               NShower_RES->Fill(ntpEvent->nshower);
00675               NTrack_RES->Fill(ntpEvent->ntrack);
00676               NShower_Vs_NTrack_RES->Fill(ntpEvent->ntrack,
00677                                          ntpEvent->nshower);
00678               NHitPlanes_RES->Fill(ntpEvent->plane.n);
00679               NTrack_Vs_EvSC_RES->Fill(ntpEvent->ph.sigcor,
00680                                       ntpEvent->ntrack);
00681               
00682               if(PassBasicCuts(j)) { //pass basic cuts loads largest track
00683                 
00684                 Float_t extraTrkStp = 0;
00685 
00686                 FracSlcStpInEvt_RES->Fill(nSliceHitInEventFrac);
00687 
00688                 NHitTrkFrac_RES->Fill(float(ntpTrack->nstrip+extraTrkStp)
00689                                      /float(ntpEvent->nstrip+extraTrkStp));
00690                 ETrkFrac_RES->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
00691                 float trkenergy = this->RecoMuEnergy(0,ntpTrack->index);
00692                 TrkMomFrac_RES->Fill(trkenergy/(trkenergy+shwenergy));
00693                 dsdz_RES->Fill(1./ntpTrack->vtx.dcosz);
00694                 NHitTrkFrac_Vs_VtxShwEnergy_RES
00695                   ->Fill(shwenergy,(float(ntpTrack->nstrip)+extraTrkStp)
00696                          /(float(ntpEvent->nstrip)+extraTrkStp));
00697                 
00698                 NHitTrkFrac_Vs_Y_RES
00699                   ->Fill(ntpTruth->y,(float(ntpTrack->nstrip)+extraTrkStp)
00700                          /(float(ntpEvent->nstrip)+extraTrkStp));
00701                 
00702                 TrkLen_RES->Fill(ntpTrack->ds);
00703                 dedx_RES->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00704                 Float_t *CCNCVars = this->CCNCSepVars(ntpTrack->index);
00705                 HalfRatio_RES->Fill(CCNCVars[2]);
00706                 delete CCNCVars;
00707                 if(this->IsFidAll(ntpTrack->index)){
00708                   RangeCurvDiff_Vs_TrkLen_RES
00709                     ->Fill(ntpTrack->ds,
00710                            2.*(this->RecoMuEnergy(1,ntpTrack->index) - 
00711                                this->RecoMuEnergy(2,ntpTrack->index))
00712                            /(this->RecoMuEnergy(1,ntpTrack->index) + 
00713                              this->RecoMuEnergy(2,ntpTrack->index)));
00714                 }
00715               }
00716               VtxShwEnergy_RES->Fill(shwenergy);
00717             }
00718           }
00719 
00720           else if(this->IResonance(mcevent)==1001){  //QE events
00721             
00722             if(theDetector==0) NEvent_QE->Fill(nEventsPerSlice[ntpEvent->slc]);
00723             else if(theDetector==1) NEvent_QE->Fill(eventSummary->nevent);
00724             
00725             if((theDetector==0 && nEventsPerSlice[ntpEvent->slc]==1) ||
00726                (theDetector==1 && eventSummary->nevent==1)){
00727               
00728               NShower_QE->Fill(ntpEvent->nshower);
00729               NTrack_QE->Fill(ntpEvent->ntrack);
00730               NShower_Vs_NTrack_QE->Fill(ntpEvent->ntrack,ntpEvent->nshower);
00731               NHitPlanes_QE->Fill(ntpEvent->plane.n);
00732               NTrack_Vs_EvSC_QE->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
00733               
00734               if(PassBasicCuts(j)) { //pass basic cuts loads largest track
00735 
00736                 FracSlcStpInEvt_QE->Fill(nSliceHitInEventFrac);
00737 
00738                 NHitTrkFrac_QE
00739                   ->Fill((float(ntpTrack->nstrip))
00740                          /(float(ntpEvent->nstrip)));
00741                 ETrkFrac_QE->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
00742                 float trkenergy = this->RecoMuEnergy(0,ntpTrack->index);
00743                 TrkMomFrac_QE->Fill(trkenergy/(trkenergy+shwenergy));
00744                 dsdz_QE->Fill(1./ntpTrack->vtx.dcosz);
00745                 NHitTrkFrac_Vs_VtxShwEnergy_QE
00746                   ->Fill(shwenergy,
00747                          (float(ntpTrack->nstrip))
00748                          /(float(ntpEvent->nstrip)));
00749                 
00750                 TrkLen_QE->Fill(ntpTrack->ds);
00751                 dedx_QE->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00752                 
00753                 Float_t *CCNCVars = this->CCNCSepVars(ntpTrack->index);
00754                 HalfRatio_QE->Fill(CCNCVars[2]);
00755                 delete CCNCVars;
00756                 if(this->IsFidAll(ntpTrack->index)){
00757                   RangeCurvDiff_Vs_TrkLen_QE->Fill(ntpTrack->ds,
00758                          2.*(this->RecoMuEnergy(1,ntpTrack->index) - 
00759                              this->RecoMuEnergy(2,ntpTrack->index))
00760                          /(this->RecoMuEnergy(1,ntpTrack->index) + 
00761                            this->RecoMuEnergy(2,ntpTrack->index)));
00762                 }
00763               }
00764               else {  //no track
00765                 VtxShwEnergy_QE->Fill(shwenergy);
00766               }
00767             }
00768           }
00769         }
00770         else if(this->IAction(mcevent)==0){
00771           
00772           if(theDetector==0) NEvent_NC->Fill(nEventsPerSlice[ntpEvent->slc]);
00773           else if(theDetector==1) NEvent_NC->Fill(eventSummary->nevent);
00774           
00775           if((theDetector==0 && nEventsPerSlice[ntpEvent->slc]==1) ||
00776              (theDetector==1 && eventSummary->nevent==1)){
00777             
00778             NShower_NC->Fill(ntpEvent->nshower);
00779             NTrack_NC->Fill(ntpEvent->ntrack);
00780             NShower_Vs_NTrack_NC->Fill(ntpEvent->ntrack,ntpEvent->nshower);
00781             NHitPlanes_NC->Fill(ntpEvent->plane.n);
00782             NTrack_Vs_EvSC_NC->Fill(ntpEvent->ph.sigcor,ntpEvent->ntrack);
00783             
00784             if(PassBasicCuts(j)) {  //pass basic cuts loads largest track
00785 
00786               FracSlcStpInEvt_NC->Fill(nSliceHitInEventFrac);
00787 
00788               NHitTrkFrac_NC
00789                 ->Fill((float(ntpTrack->nstrip))
00790                        /(float(ntpEvent->nstrip)));
00791               ETrkFrac_NC->Fill(ntpTrack->ph.sigcor/ntpEvent->ph.sigcor);
00792               float trkenergy = this->RecoMuEnergy(0,ntpTrack->index);
00793               TrkMomFrac_NC->Fill(trkenergy/(trkenergy+shwenergy));
00794               dsdz_NC->Fill(1./ntpTrack->vtx.dcosz);
00795               NHitTrkFrac_Vs_VtxShwEnergy_NC
00796                 ->Fill(shwenergy,
00797                        (float(ntpTrack->nstrip))
00798                        /(float(ntpEvent->nstrip)));
00799               TrkLen_NC->Fill(ntpTrack->ds);
00800               dedx_NC->Fill(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
00801               Float_t *CCNCVars = this->CCNCSepVars(ntpTrack->index);
00802               HalfRatio_NC->Fill(CCNCVars[2]);
00803               delete CCNCVars;
00804               if(this->IsFidAll(ntpTrack->index)){
00805                 RangeCurvDiff_Vs_TrkLen_NC->Fill(ntpTrack->ds,
00806                        2.*(this->RecoMuEnergy(1,ntpTrack->index) - 
00807                            this->RecoMuEnergy(2,ntpTrack->index))
00808                        /(this->RecoMuEnergy(1,ntpTrack->index) + 
00809                          this->RecoMuEnergy(2,ntpTrack->index)));
00810               }
00811             }
00812             VtxShwEnergy_NC->Fill(shwenergy);
00813           }
00814         }
00815       }
00816     }
00817     delete [] nEventsPerSlice;
00818   }
00819 
00820   //#close file  
00821   save.Write();
00822   save.Close();
00823   
00824 }
00825 
00826 
00827 //*******************************************************
00828 void MadCBSQEAnalysis::TestQEDiscrim(std::string tag){
00829 //awoifids
00830 
00831   float integ = 0.;
00832   //read in input histos and normalise them
00833   char nam[256];
00834   sprintf(nam,"QEHistos_%s.root",tag.c_str());
00835   TFile *QEFile = new TFile(nam,"READ");
00836   TH2F *NShower_Vs_NTrack_NC = (TH2F*) QEFile->Get("NShower_Vs_NTrack_NC");
00837   TH2F *NShower_Vs_NTrack_DIS = (TH2F*) QEFile->Get("NShower_Vs_NTrack_DIS");
00838   TH2F *NShower_Vs_NTrack_RES = (TH2F*) QEFile->Get("NShower_Vs_NTrack_RES");
00839   TH2F *NShower_Vs_NTrack_QE = (TH2F*) QEFile->Get("NShower_Vs_NTrack_QE");
00840   integ = NShower_Vs_NTrack_NC->Integral(); 
00841   NShower_Vs_NTrack_NC->Scale(1./integ);
00842   integ = NShower_Vs_NTrack_DIS->Integral(); 
00843   NShower_Vs_NTrack_DIS->Scale(1./integ);
00844   integ = NShower_Vs_NTrack_RES->Integral(); 
00845   NShower_Vs_NTrack_RES->Scale(1./integ);
00846   integ = NShower_Vs_NTrack_QE->Integral(); 
00847   NShower_Vs_NTrack_QE->Scale(1./integ);
00848   
00849   TH1F *NHitPlanes_NC = (TH1F*) QEFile->Get("NHitPlanes_NC");
00850   TH1F *NHitPlanes_DIS = (TH1F*) QEFile->Get("NHitPlanes_DIS");
00851   TH1F *NHitPlanes_RES = (TH1F*) QEFile->Get("NHitPlanes_RES");
00852   TH1F *NHitPlanes_QE = (TH1F*) QEFile->Get("NHitPlanes_QE");
00853   integ = NHitPlanes_NC->Integral(); 
00854   NHitPlanes_NC->Scale(1./integ);
00855   integ = NHitPlanes_DIS->Integral(); 
00856   NHitPlanes_DIS->Scale(1./integ);
00857   integ = NHitPlanes_RES->Integral(); 
00858   NHitPlanes_RES->Scale(1./integ);
00859   integ = NHitPlanes_QE->Integral(); 
00860   NHitPlanes_QE->Scale(1./integ);
00861 
00862   TH1F *NHitTrkFrac_NC = (TH1F*) QEFile->Get("NHitTrkFrac_NC");
00863   TH1F *NHitTrkFrac_DIS = (TH1F*) QEFile->Get("NHitTrkFrac_DIS");
00864   TH1F *NHitTrkFrac_RES = (TH1F*) QEFile->Get("NHitTrkFrac_RES");
00865   TH1F *NHitTrkFrac_QE = (TH1F*) QEFile->Get("NHitTrkFrac_QE");
00866   integ = NHitTrkFrac_NC->Integral(); 
00867   NHitTrkFrac_NC->Scale(1./integ);
00868   integ = NHitTrkFrac_DIS->Integral(); 
00869   NHitTrkFrac_DIS->Scale(1./integ);
00870   integ = NHitTrkFrac_RES->Integral(); 
00871   NHitTrkFrac_RES->Scale(1./integ);
00872   integ = NHitTrkFrac_QE->Integral(); 
00873   NHitTrkFrac_QE->Scale(1./integ);
00874 
00875   TH1F *dsdz_NC = (TH1F*) QEFile->Get("dsdz_NC");
00876   TH1F *dsdz_DIS = (TH1F*) QEFile->Get("dsdz_DIS");
00877   TH1F *dsdz_RES = (TH1F*) QEFile->Get("dsdz_RES");
00878   TH1F *dsdz_QE = (TH1F*) QEFile->Get("dsdz_QE");
00879   integ = dsdz_NC->Integral(); 
00880   dsdz_NC->Scale(1./integ);
00881   integ = dsdz_DIS->Integral(); 
00882   dsdz_DIS->Scale(1./integ);
00883   integ = dsdz_RES->Integral(); 
00884   dsdz_RES->Scale(1./integ);
00885   integ = dsdz_QE->Integral(); 
00886   dsdz_QE->Scale(1./integ);
00887 
00888   TH1F *dedx_NC = (TH1F*) QEFile->Get("dedx_NC");
00889   TH1F *dedx_DIS = (TH1F*) QEFile->Get("dedx_DIS");
00890   TH1F *dedx_RES = (TH1F*) QEFile->Get("dedx_RES");
00891   TH1F *dedx_QE = (TH1F*) QEFile->Get("dedx_QE");
00892   integ = dedx_NC->Integral(); 
00893   dedx_NC->Scale(1./integ);
00894   integ = dedx_DIS->Integral(); 
00895   dedx_DIS->Scale(1./integ);
00896   integ = dedx_RES->Integral(); 
00897   dedx_RES->Scale(1./integ);
00898   integ = dedx_QE->Integral(); 
00899   dedx_QE->Scale(1./integ);
00900 
00901   TH1F *HalfRatio_NC = (TH1F*) QEFile->Get("HalfRatio_NC");
00902   TH1F *HalfRatio_DIS = (TH1F*) QEFile->Get("HalfRatio_DIS");
00903   TH1F *HalfRatio_RES = (TH1F*) QEFile->Get("HalfRatio_RES");
00904   TH1F *HalfRatio_QE = (TH1F*) QEFile->Get("HalfRatio_QE");
00905   integ = HalfRatio_NC->Integral(); 
00906   HalfRatio_NC->Scale(1./integ);
00907   integ = HalfRatio_DIS->Integral(); 
00908   HalfRatio_DIS->Scale(1./integ);
00909   integ = HalfRatio_RES->Integral(); 
00910   HalfRatio_RES->Scale(1./integ);
00911   integ = HalfRatio_QE->Integral(); 
00912   HalfRatio_QE->Scale(1./integ);
00913 
00914   TH2F *RangeCurvDiff_Vs_TrkLen_NC = (TH2F*) 
00915     QEFile->Get("RangeCurvDiff_Vs_TrkLen_NC");
00916   TH1D *RangeCurvDiff_NC = RangeCurvDiff_Vs_TrkLen_NC->ProjectionY();
00917   TH2F *RangeCurvDiff_Vs_TrkLen_DIS = (TH2F*) 
00918     QEFile->Get("RangeCurvDiff_Vs_TrkLen_DIS");
00919   TH1D *RangeCurvDiff_DIS = RangeCurvDiff_Vs_TrkLen_DIS->ProjectionY();
00920   TH2F *RangeCurvDiff_Vs_TrkLen_RES = (TH2F*) 
00921     QEFile->Get("RangeCurvDiff_Vs_TrkLen_RES");
00922   TH1D *RangeCurvDiff_RES = RangeCurvDiff_Vs_TrkLen_RES->ProjectionY();
00923   TH2F *RangeCurvDiff_Vs_TrkLen_QE = (TH2F*) 
00924     QEFile->Get("RangeCurvDiff_Vs_TrkLen_QE");
00925   TH1D *RangeCurvDiff_QE = RangeCurvDiff_Vs_TrkLen_QE->ProjectionY();
00926   integ = RangeCurvDiff_NC->Integral(); 
00927   RangeCurvDiff_NC->Scale(1./integ);
00928   integ = RangeCurvDiff_DIS->Integral(); 
00929   RangeCurvDiff_DIS->Scale(1./integ);
00930   integ = RangeCurvDiff_RES->Integral(); 
00931   RangeCurvDiff_RES->Scale(1./integ);
00932   integ = RangeCurvDiff_QE->Integral(); 
00933   RangeCurvDiff_QE->Scale(1./integ);
00934 
00935   //output tree
00936   Float_t ProbNC[8];
00937   Float_t ProbDIS[8];
00938   Float_t ProbRES[8];
00939   Float_t ProbQE[8];
00940   Float_t PidParQEDIS[8];
00941   Float_t PidParQERES[8];
00942   Float_t PidParQENC[8];
00943   Float_t PidParDISNC[8];
00944   Int_t IAction;
00945   Int_t INu;
00946   Int_t Process;
00947   Float_t Y;
00948   Float_t W2;
00949   Float_t NuEn;
00950   Float_t QENuEn;
00951 
00952   std::string savename = "TestQEDiscrim_" + tag + ".root";
00953   TFile *save = new TFile(savename.c_str(),"RECREATE"); 
00954   TTree *tree = new TTree("QETree","QETree");
00955   tree->Branch("ProbNC",&ProbNC,"ProbNC[8]/F",32000);
00956   tree->Branch("ProbDIS",&ProbDIS,"ProbDIS[8]/F",32000);
00957   tree->Branch("ProbRES",&ProbRES,"ProbRES[8]/F",32000);
00958   tree->Branch("ProbQE",&ProbQE,"ProbQE[8]/F",32000);
00959   tree->Branch("PidParQEDIS",&PidParQEDIS,"PidParQEDIS[8]/F",32000);
00960   tree->Branch("PidParQERES",&PidParQERES,"PidParQERES[8]/F",32000);
00961   tree->Branch("PidParQENC",&PidParQENC,"PidParQENC[8]/F",32000);
00962   tree->Branch("PidParDISNC",&PidParDISNC,"PidParDISNC[8]/F",32000);
00963   tree->Branch("IAction",&IAction,"IAction/I",32000);
00964   tree->Branch("INu",&INu,"INu/I",32000);
00965   tree->Branch("Process",&Process,"Process/I",32000);
00966   tree->Branch("Y",&Y,"Y/F",32000);
00967   tree->Branch("W2",&W2,"W2/F",32000);
00968   tree->Branch("NuEn",&NuEn,"NuEn/F",32000);
00969   tree->Branch("QENuEn",&QENuEn,"QENuEn/F",32000);
00970   
00971   for(int i=0;i<Nentries;i++){ //Nentries
00972     
00973     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00974                             << "% done" << std::endl;
00975     if(!this->GetEntry(i)) continue;
00976 
00977     //get which detector this snarl is from
00978     int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
00979 
00980     //Calculate number of events per slice for ND:
00981     Int_t *nEventsPerSlice = new Int_t[eventSummary->nslice];
00982     for(int i=0;i<eventSummary->nslice;i++) nEventsPerSlice[i]= 0;
00983     if(theDetector==0) {
00984       for(int i=0;i<eventSummary->nslice;i++){
00985         nEventsPerSlice[i]= 0;
00986         if(!LoadEvent(i)) continue;
00987         nEventsPerSlice[ntpEvent->slc] += 1;
00988       }
00989     }
00990   
00991     for(int j=0;j<eventSummary->nevent;j++){ //events
00992       
00993       if(!LoadEvent(j)) continue;
00994       if(!IsFidVtxEvt(j)) continue;
00995 
00996       //if no good truth candidate, we loose event
00997       int mcevent = -1;
00998       if(LoadTruthForRecoTH(j,mcevent)){  //load truth
00999         
01000         if((theDetector==0 && nEventsPerSlice[ntpEvent->slc]==1) || 
01001            (theDetector==1 && eventSummary->nevent==1)){  //1 event only
01002         
01003           //reset probs
01004           for(int k=0;k<8;k++){
01005             ProbNC[k] = 1;
01006             ProbDIS[k] = 1;
01007             ProbRES[k] = 1;
01008             ProbQE[k] = 1;
01009           }
01010         
01011           //NC:
01012           Int_t bin1 = NShower_Vs_NTrack_NC->FindBin(ntpEvent->ntrack,
01013                                                      ntpEvent->nshower);
01014           // ProbNC[0] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01015           // ProbNC[2] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01016           // ProbNC[3] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01017           // ProbNC[4] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01018           // ProbNC[5] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01019           // ProbNC[6] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01020           // ProbNC[7] *= NShower_Vs_NTrack_NC->GetBinContent(bin1);
01021           
01022           // bin1 = NHitPlanes_NC->FindBin(ntpEvent->plane.n);
01023           // ProbNC[0] *= NHitPlanes_NC->GetBinContent(bin1);
01024           // ProbNC[1] *= NHitPlanes_NC->GetBinContent(bin1);
01025           // ProbNC[3] *= NHitPlanes_NC->GetBinContent(bin1);
01026           // ProbNC[4] *= NHitPlanes_NC->GetBinContent(bin1);
01027           // ProbNC[5] *= NHitPlanes_NC->GetBinContent(bin1);
01028           // ProbNC[6] *= NHitPlanes_NC->GetBinContent(bin1);
01029           // ProbNC[7] *= NHitPlanes_NC->GetBinContent(bin1);
01030           
01031           //DIS:    
01032           Int_t bin2 = NShower_Vs_NTrack_DIS->FindBin(ntpEvent->ntrack,
01033                                                       ntpEvent->nshower);
01034           // ProbDIS[0] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01035           // ProbDIS[2] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01036           // ProbDIS[3] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01037           // ProbDIS[4] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01038           // ProbDIS[5] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01039           // ProbDIS[6] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01040           // ProbDIS[7] *= NShower_Vs_NTrack_DIS->GetBinContent(bin2);
01041           
01042           // bin2 = NHitPlanes_DIS->FindBin(ntpEvent->plane.n);
01043           // ProbDIS[0] *= NHitPlanes_DIS->GetBinContent(bin2);
01044           // ProbDIS[1] *= NHitPlanes_DIS->GetBinContent(bin2);
01045           // ProbDIS[3] *= NHitPlanes_DIS->GetBinContent(bin2);
01046           // ProbDIS[4] *= NHitPlanes_DIS->GetBinContent(bin2);
01047           // ProbDIS[5] *= NHitPlanes_DIS->GetBinContent(bin2);
01048           // ProbDIS[6] *= NHitPlanes_DIS->GetBinContent(bin2);
01049           // ProbDIS[7] *= NHitPlanes_DIS->GetBinContent(bin2);
01050 
01051           //RES:
01052           Int_t bin3 = NShower_Vs_NTrack_RES->FindBin(ntpEvent->ntrack,
01053                                                 ntpEvent->nshower);
01054           // ProbRES[0] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01055           // ProbRES[2] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01056           // ProbRES[3] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01057           // ProbRES[4] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01058           // ProbRES[5] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01059           // ProbRES[6] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01060           // ProbRES[7] *= NShower_Vs_NTrack_RES->GetBinContent(bin3);
01061           
01062           // bin3 = NHitPlanes_RES->FindBin(ntpEvent->plane.n);
01063           // ProbRES[0] *= NHitPlanes_RES->GetBinContent(bin3);
01064           // ProbRES[1] *= NHitPlanes_RES->GetBinContent(bin3);
01065           // ProbRES[3] *= NHitPlanes_RES->GetBinContent(bin3);
01066           // ProbRES[4] *= NHitPlanes_RES->GetBinContent(bin3);
01067           // ProbRES[5] *= NHitPlanes_RES->GetBinContent(bin3);
01068           // ProbRES[6] *= NHitPlanes_RES->GetBinContent(bin3);
01069           // ProbRES[7] *= NHitPlanes_RES->GetBinContent(bin3);
01070           
01071           //QE:
01072           Int_t bin4 = NShower_Vs_NTrack_QE->FindBin(ntpEvent->ntrack,
01073                                                      ntpEvent->nshower);
01074           // ProbQE[0] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01075           // ProbQE[2] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01076           // ProbQE[3] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01077           // ProbQE[4] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01078           // ProbQE[5] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01079           // ProbQE[6] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01080           // ProbQE[7] *= NShower_Vs_NTrack_QE->GetBinContent(bin4);
01081                 
01082           // bin4 = NHitPlanes_QE->FindBin(ntpEvent->plane.n);
01083           // ProbQE[0] *= NHitPlanes_QE->GetBinContent(bin4);
01084           // ProbQE[1] *= NHitPlanes_QE->GetBinContent(bin4);
01085           // ProbQE[3] *= NHitPlanes_QE->GetBinContent(bin4);
01086           // ProbQE[4] *= NHitPlanes_QE->GetBinContent(bin4);
01087           // ProbQE[5] *= NHitPlanes_QE->GetBinContent(bin4);
01088           // ProbQE[6] *= NHitPlanes_QE->GetBinContent(bin4);
01089           // ProbQE[7] *= NHitPlanes_QE->GetBinContent(bin4);
01090         
01091         
01092           if(PassBasicCuts(j)) {  //loads largest track
01093             Float_t *CCNCVars = this->CCNCSepVars(ntpTrack->index);      
01094           
01095             //NC:
01096             bin1 = NHitTrkFrac_NC->FindBin(float(ntpTrack->nstrip)
01097                                            /float(ntpEvent->nstrip));
01098             ProbNC[0] *= NHitTrkFrac_NC->GetBinContent(bin1);
01099             ProbNC[1] *= NHitTrkFrac_NC->GetBinContent(bin1);
01100             ProbNC[2] *= NHitTrkFrac_NC->GetBinContent(bin1);
01101             ProbNC[4] *= NHitTrkFrac_NC->GetBinContent(bin1);
01102             ProbNC[5] *= NHitTrkFrac_NC->GetBinContent(bin1);
01103             ProbNC[6] *= NHitTrkFrac_NC->GetBinContent(bin1);
01104             ProbNC[7] *= NHitTrkFrac_NC->GetBinContent(bin1);
01105           
01106             bin1 = dedx_NC->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
01107             ProbNC[0] *= dedx_NC->GetBinContent(bin1);
01108             ProbNC[1] *= dedx_NC->GetBinContent(bin1);
01109             ProbNC[2] *= dedx_NC->GetBinContent(bin1);
01110             ProbNC[3] *= dedx_NC->GetBinContent(bin1);
01111             ProbNC[5] *= dedx_NC->GetBinContent(bin1);
01112             ProbNC[6] *= dedx_NC->GetBinContent(bin1);
01113             ProbNC[7] *= dedx_NC->GetBinContent(bin1);
01114             
01115             /*
01116               bin1 = dsdz_NC->FindBin(1./ntpTrack->vtx.dcosz);
01117               ProbNC[0] *= dsdz_NC->GetBinContent(bin1);
01118               ProbNC[1] *= dsdz_NC->GetBinContent(bin1);
01119               ProbNC[2] *= dsdz_NC->GetBinContent(bin1);
01120               ProbNC[3] *= dsdz_NC->GetBinContent(bin1);
01121               ProbNC[4] *= dsdz_NC->GetBinContent(bin1);
01122               ProbNC[6] *= dsdz_NC->GetBinContent(bin1);
01123               ProbNC[7] *= dsdz_NC->GetBinContent(bin1);
01124             */
01125           
01126             bin1 = HalfRatio_NC->FindBin(CCNCVars[2]);
01127             ProbNC[0] *= HalfRatio_NC->GetBinContent(bin1);
01128             ProbNC[1] *= HalfRatio_NC->GetBinContent(bin1);
01129             ProbNC[2] *= HalfRatio_NC->GetBinContent(bin1);
01130             ProbNC[3] *= HalfRatio_NC->GetBinContent(bin1);
01131             ProbNC[4] *= HalfRatio_NC->GetBinContent(bin1);
01132             ProbNC[5] *= HalfRatio_NC->GetBinContent(bin1);
01133             ProbNC[7] *= HalfRatio_NC->GetBinContent(bin1);
01134             
01135             //DIS:
01136             bin2 = NHitTrkFrac_DIS->FindBin(float(ntpTrack->nstrip)
01137                                             /float(ntpEvent->nstrip));
01138             ProbDIS[0] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01139             ProbDIS[1] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01140             ProbDIS[2] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01141             ProbDIS[4] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01142             ProbDIS[5] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01143             ProbDIS[6] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01144             ProbDIS[7] *= NHitTrkFrac_DIS->GetBinContent(bin2);
01145             
01146             bin2 = dedx_DIS->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
01147             ProbDIS[0] *= dedx_DIS->GetBinContent(bin2);
01148             ProbDIS[1] *= dedx_DIS->GetBinContent(bin2);
01149             ProbDIS[2] *= dedx_DIS->GetBinContent(bin2);
01150             ProbDIS[3] *= dedx_DIS->GetBinContent(bin2);
01151             ProbDIS[5] *= dedx_DIS->GetBinContent(bin2);
01152             ProbDIS[6] *= dedx_DIS->GetBinContent(bin2);
01153             ProbDIS[7] *= dedx_DIS->GetBinContent(bin2);
01154             
01155             // bin2 = dsdz_DIS->FindBin(1./ntpTrack->vtx.dcosz);
01156             // ProbDIS[0] *= dsdz_DIS->GetBinContent(bin2);
01157             // ProbDIS[1] *= dsdz_DIS->GetBinContent(bin2);
01158             // ProbDIS[2] *= dsdz_DIS->GetBinContent(bin2);
01159             // ProbDIS[3] *= dsdz_DIS->GetBinContent(bin2);
01160             // ProbDIS[4] *= dsdz_DIS->GetBinContent(bin2);
01161             // ProbDIS[6] *= dsdz_DIS->GetBinContent(bin2);
01162             // ProbDIS[7] *= dsdz_DIS->GetBinContent(bin2);
01163             
01164             bin2 = HalfRatio_DIS->FindBin(CCNCVars[2]);
01165             ProbDIS[0] *= HalfRatio_DIS->GetBinContent(bin2);
01166             ProbDIS[1] *= HalfRatio_DIS->GetBinContent(bin2);
01167             ProbDIS[2] *= HalfRatio_DIS->GetBinContent(bin2);
01168             ProbDIS[3] *= HalfRatio_DIS->GetBinContent(bin2);
01169             ProbDIS[4] *= HalfRatio_DIS->GetBinContent(bin2);
01170             ProbDIS[5] *= HalfRatio_DIS->GetBinContent(bin2);
01171             ProbDIS[7] *= HalfRatio_DIS->GetBinContent(bin2);
01172             
01173             //RES:
01174             bin3 = NHitTrkFrac_RES->FindBin(float(ntpTrack->nstrip)
01175                                             /float(ntpEvent->nstrip));
01176             ProbRES[0] *= NHitTrkFrac_RES->GetBinContent(bin3);
01177             ProbRES[1] *= NHitTrkFrac_RES->GetBinContent(bin3);
01178             ProbRES[2] *= NHitTrkFrac_RES->GetBinContent(bin3);
01179             ProbRES[4] *= NHitTrkFrac_RES->GetBinContent(bin3);
01180             ProbRES[5] *= NHitTrkFrac_RES->GetBinContent(bin3);
01181             ProbRES[6] *= NHitTrkFrac_RES->GetBinContent(bin3);
01182             ProbRES[7] *= NHitTrkFrac_RES->GetBinContent(bin3);
01183             
01184             bin3 = dedx_RES->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
01185             ProbRES[0] *= dedx_RES->GetBinContent(bin3);
01186             ProbRES[1] *= dedx_RES->GetBinContent(bin3);
01187             ProbRES[2] *= dedx_RES->GetBinContent(bin3);
01188             ProbRES[3] *= dedx_RES->GetBinContent(bin3);
01189             ProbRES[5] *= dedx_RES->GetBinContent(bin3);
01190             ProbRES[6] *= dedx_RES->GetBinContent(bin3);
01191             ProbRES[7] *= dedx_RES->GetBinContent(bin3);
01192             
01193             // bin3 = dsdz_RES->FindBin(1./ntpTrack->vtx.dcosz);
01194             // ProbRES[0] *= dsdz_RES->GetBinContent(bin3);
01195             // ProbRES[1] *= dsdz_RES->GetBinContent(bin3);
01196             // ProbRES[2] *= dsdz_RES->GetBinContent(bin3);
01197             // ProbRES[3] *= dsdz_RES->GetBinContent(bin3);
01198             // ProbRES[4] *= dsdz_RES->GetBinContent(bin3);
01199             // ProbRES[6] *= dsdz_RES->GetBinContent(bin3);
01200             // ProbRES[7] *= dsdz_RES->GetBinContent(bin3);
01201             
01202             bin3 = HalfRatio_RES->FindBin(CCNCVars[2]);
01203             ProbRES[0] *= HalfRatio_RES->GetBinContent(bin3);
01204             ProbRES[1] *= HalfRatio_RES->GetBinContent(bin3);
01205             ProbRES[2] *= HalfRatio_RES->GetBinContent(bin3);
01206             ProbRES[3] *= HalfRatio_RES->GetBinContent(bin3);
01207             ProbRES[4] *= HalfRatio_RES->GetBinContent(bin3);
01208             ProbRES[5] *= HalfRatio_RES->GetBinContent(bin3);
01209             ProbRES[7] *= HalfRatio_RES->GetBinContent(bin3);
01210             
01211             //QE:
01212             bin4 = NHitTrkFrac_QE->FindBin(float(ntpTrack->nstrip)
01213                                            /float(ntpEvent->nstrip));
01214             ProbQE[0] *= NHitTrkFrac_QE->GetBinContent(bin4);
01215             ProbQE[1] *= NHitTrkFrac_QE->GetBinContent(bin4);
01216             ProbQE[2] *= NHitTrkFrac_QE->GetBinContent(bin4);
01217             ProbQE[4] *= NHitTrkFrac_QE->GetBinContent(bin4);
01218             ProbQE[5] *= NHitTrkFrac_QE->GetBinContent(bin4);
01219             ProbQE[6] *= NHitTrkFrac_QE->GetBinContent(bin4);
01220             ProbQE[7] *= NHitTrkFrac_QE->GetBinContent(bin4);
01221             
01222             bin4 = dedx_QE->FindBin(ntpTrack->ph.sigcor/(100.*ntpTrack->ds));
01223             ProbQE[0] *= dedx_QE->GetBinContent(bin4);
01224             ProbQE[1] *= dedx_QE->GetBinContent(bin4);
01225             ProbQE[2] *= dedx_QE->GetBinContent(bin4);
01226             ProbQE[3] *= dedx_QE->GetBinContent(bin4);
01227             ProbQE[5] *= dedx_QE->GetBinContent(bin4);
01228             ProbQE[6] *= dedx_QE->GetBinContent(bin4);
01229             ProbQE[7] *= dedx_QE->GetBinContent(bin4);
01230             
01231             // bin4 = dsdz_QE->FindBin(1./ntpTrack->vtx.dcosz);
01232             // ProbQE[0] *= dsdz_QE->GetBinContent(bin4);
01233             // ProbQE[1] *= dsdz_QE->GetBinContent(bin4);
01234             // ProbQE[2] *= dsdz_QE->GetBinContent(bin4);
01235             // ProbQE[3] *= dsdz_QE->GetBinContent(bin4);
01236             // ProbQE[4] *= dsdz_QE->GetBinContent(bin4);
01237             // ProbQE[6] *= dsdz_QE->GetBinContent(bin4);
01238             // ProbQE[7] *= dsdz_QE->GetBinContent(bin4);
01239             
01240             bin4 = HalfRatio_QE->FindBin(CCNCVars[2]);
01241             ProbQE[0] *= HalfRatio_QE->GetBinContent(bin4);
01242             ProbQE[1] *= HalfRatio_QE->GetBinContent(bin4);
01243             ProbQE[2] *= HalfRatio_QE->GetBinContent(bin4);
01244             ProbQE[3] *= HalfRatio_QE->GetBinContent(bin4);
01245             ProbQE[4] *= HalfRatio_QE->GetBinContent(bin4);
01246             ProbQE[5] *= HalfRatio_QE->GetBinContent(bin4);
01247             ProbQE[7] *= HalfRatio_QE->GetBinContent(bin4);
01248           
01249             if(false){
01250               if(this->IsFidAll(ntpTrack->index)){
01251                 Float_t fracDiff = 2.*(RecoMuEnergy(1,ntpTrack->index) - 
01252                                        RecoMuEnergy(2,ntpTrack->index))
01253                   /(RecoMuEnergy(1,ntpTrack->index) + 
01254                     RecoMuEnergy(2,ntpTrack->index));
01255                 
01256                 //NC:
01257                 bin1 = RangeCurvDiff_NC->FindBin(fracDiff);
01258                 ProbNC[0] *= RangeCurvDiff_NC->GetBinContent(bin1);
01259                 ProbNC[1] *= RangeCurvDiff_NC->GetBinContent(bin1);
01260                 ProbNC[2] *= RangeCurvDiff_NC->GetBinContent(bin1);
01261                 ProbNC[3] *= RangeCurvDiff_NC->GetBinContent(bin1);
01262                 ProbNC[4] *= RangeCurvDiff_NC->GetBinContent(bin1);
01263                 ProbNC[5] *= RangeCurvDiff_NC->GetBinContent(bin1);
01264                 ProbNC[6] *= RangeCurvDiff_NC->GetBinContent(bin1);
01265               
01266                 //DIS:
01267                 bin2 = RangeCurvDiff_DIS->FindBin(fracDiff);
01268                 ProbDIS[0] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01269                 ProbDIS[1] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01270                 ProbDIS[2] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01271                 ProbDIS[3] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01272                 ProbDIS[4] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01273                 ProbDIS[5] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01274                 ProbDIS[6] *= RangeCurvDiff_DIS->GetBinContent(bin2);
01275                 
01276                 //RES:
01277                 bin3 = RangeCurvDiff_RES->FindBin(fracDiff);
01278                 ProbRES[0] *= RangeCurvDiff_RES->GetBinContent(bin3);
01279                 ProbRES[1] *= RangeCurvDiff_RES->GetBinContent(bin3);
01280                 ProbRES[2] *= RangeCurvDiff_RES->GetBinContent(bin3);
01281                 ProbRES[3] *= RangeCurvDiff_RES->GetBinContent(bin3);
01282                 ProbRES[4] *= RangeCurvDiff_RES->GetBinContent(bin3);
01283                 ProbRES[5] *= RangeCurvDiff_RES->GetBinContent(bin3);
01284                 ProbRES[6] *= RangeCurvDiff_RES->GetBinContent(bin3);
01285                 
01286                 //QE:
01287                 bin4 = RangeCurvDiff_QE->FindBin(fracDiff);
01288                 ProbQE[0] *= RangeCurvDiff_QE->GetBinContent(bin4);
01289                 ProbQE[1] *= RangeCurvDiff_QE->GetBinContent(bin4);
01290                 ProbQE[2] *= RangeCurvDiff_QE->GetBinContent(bin4);
01291                 ProbQE[3] *= RangeCurvDiff_QE->GetBinContent(bin4);
01292                 ProbQE[4] *= RangeCurvDiff_QE->GetBinContent(bin4);
01293                 ProbQE[5] *= RangeCurvDiff_QE->GetBinContent(bin4);
01294                 ProbQE[6] *= RangeCurvDiff_QE->GetBinContent(bin4);
01295               }
01296             }
01297             delete CCNCVars;
01298           }
01299         }
01300         
01301         //pid pars
01302         for(int j=0;j<8;j++){
01303           if(ProbQE[j]!=0&&ProbDIS[j]!=0) 
01304             PidParQEDIS[j] = sqrt(-TMath::Log(ProbQE[j]))
01305               -sqrt(-TMath::Log(ProbDIS[j]));
01306           else if(ProbQE[j]==0&&ProbDIS[j]==0) PidParQEDIS[j]=10;
01307           else if(ProbQE[j]==0) PidParQEDIS[j] = 10.-sqrt(-TMath::Log(ProbDIS[j]));
01308           else if(ProbDIS[j]==0) PidParQEDIS[j] = sqrt(-TMath::Log(ProbQE[j]))-10.;
01309           
01310           if(ProbQE[j]!=0&&ProbRES[j]!=0) 
01311             PidParQERES[j] = sqrt(-TMath::Log(ProbQE[j]))
01312               -sqrt(-TMath::Log(ProbRES[j]));
01313           else if(ProbQE[j]==0&&ProbRES[j]==0) PidParQERES[j]=10;
01314           else if(ProbQE[j]==0) PidParQERES[j] = 10.-sqrt(-TMath::Log(ProbRES[j]));
01315           else if(ProbRES[j]==0) PidParQERES[j] = sqrt(-TMath::Log(ProbQE[j]))-10.;
01316           
01317           if(ProbQE[j]!=0&&ProbNC[j]!=0) 
01318             PidParQENC[j] = sqrt(-TMath::Log(ProbQE[j]))
01319               -sqrt(-TMath::Log(ProbNC[j]));
01320           else if(ProbQE[j]==0&&ProbNC[j]==0) PidParQENC[j]=10;
01321           else if(ProbQE[j]==0) PidParQENC[j] = 10.-sqrt(-TMath::Log(ProbNC[j]));
01322           else if(ProbNC[j]==0) PidParQENC[j] = sqrt(-TMath::Log(ProbQE[j]))-10.;
01323           
01324           if(ProbDIS[j]!=0&&ProbNC[j]!=0) 
01325             PidParDISNC[j] = sqrt(-TMath::Log(ProbDIS[j]))
01326               -sqrt(-TMath::Log(ProbNC[j]));
01327           else if(ProbDIS[j]==0&&ProbNC[j]==0) PidParDISNC[j]=10;
01328           else if(ProbDIS[j]==0) PidParDISNC[j] = 10.-sqrt(-TMath::Log(ProbNC[j]));
01329           else if(ProbNC[j]==0) PidParDISNC[j] = sqrt(-TMath::Log(ProbDIS[j]))-10.;
01330         }
01331         
01332         IAction = ntpTruth->iaction;
01333         INu = ntpTruth->inu;
01334         Y = ntpTruth->y;
01335         W2 = ntpTruth->w2;
01336         Process = ntpTruth->iresonance;
01337         NuEn = ntpTruth->p4neu[3];
01338         if(PassTrackCuts(ntpTrack->index)) QENuEn = RecoQENuEnergy(ntpTrack->index);
01339         else QENuEn = -1;
01340         tree->Fill();
01341       }
01342     }
01343     delete [] nEventsPerSlice;
01344   }
01345   save->cd();
01346   tree->Write();
01347   save->Close();
01348 }
01349 
01350 //********************************************************
01351 void MadCBSQEAnalysis::ReadPIDFile(std::string tag){
01352   fLikeliFile = new TFile(tag.c_str(),"READ");
01353   if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) {
01354     fLikeliFile->cd();
01355     fLikeliHist[0] = (TH1F*) fLikeliFile->Get("NHitTrkFrac_NC");
01356     fLikeliHist[1] = (TH1F*) fLikeliFile->Get("dedx_NC");
01357     fLikeliHist[2] = (TH1F*) fLikeliFile->Get("HalfRatio_NC");
01358     fLikeliHist[3] = (TH1F*) fLikeliFile->Get("NHitTrkFrac_DIS");
01359     fLikeliHist[4] = (TH1F*) fLikeliFile->Get("dedx_DIS");
01360     fLikeliHist[5] = (TH1F*) fLikeliFile->Get("HalfRatio_DIS");
01361     fLikeliHist[6] = (TH1F*) fLikeliFile->Get("NHitTrkFrac_RES");
01362     fLikeliHist[7] = (TH1F*) fLikeliFile->Get("dedx_RES");
01363     fLikeliHist[8] = (TH1F*) fLikeliFile->Get("HalfRatio_RES");
01364     fLikeliHist[9] = (TH1F*) fLikeliFile->Get("NHitTrkFrac_QE");
01365     fLikeliHist[10] = (TH1F*) fLikeliFile->Get("dedx_QE");
01366     fLikeliHist[11] = (TH1F*) fLikeliFile->Get("HalfRatio_QE");
01367     for(int i=0;i<12;i++){
01368       float integ = fLikeliHist[i]->Integral(); 
01369       fLikeliHist[i]->Scale(1./integ);
01370     }
01371   }
01372   else fLikeliFile = NULL;
01373 }
01374 
01375 #endif // #ifdef madcbsqeanalysis_cxx