MadEdAnalysis Class Reference

#include <MadEdAnalysis.h>

Inheritance diagram for MadEdAnalysis:
MadAnalysis MadQuantities MadBase

List of all members.

Public Member Functions

 MadEdAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadEdAnalysis (JobC *, string, int)
 ~MadEdAnalysis ()
Bool_t MyIsFidVtx (Int_t itrk=0)
Bool_t MyIsFidVtxrz (Int_t itrk=0)
Float_t PID (Int_t event=0, Int_t method=0)
Float_t NeuNetEval (Int_t event=0)
Float_t * LikeliPID ()
Float_t * MyLikeliQE (TH1F **)
TMultiLayerPerceptron * NeuNetTrain (TFile *)
void SetMLP (TMultiLayerPerceptron *neural)
void PIDHisto ()
void MyMakeQEFile (std::string)
void DataHist (std::string)
void MCHist (std::string)
void MakeEff (std::string)
void MyMakeMyFile (std::string)
void MyReadPIDFile (std::string)
void MyCreatePAN (std::string tag)
void MyCreatePANData (std::string tag)

Protected Member Functions

Bool_t PassAnalysisCuts (Int_t event=0)
Bool_t PassBasicCuts (Int_t event=0)
Float_t HitF (Int_t itr)
Float_t ETrkF (Int_t itr)
Float_t DeDx (Int_t itr)
Float_t EvtLength ()
Float_t Trkphsig (Int_t itrk)
Float_t Trkplanes (Int_t itrk)
Float_t Evtphsig ()
Float_t RecoMuDCosY (Int_t itrk)
Float_t RecoMuZn (Int_t itrk)
Float_t RecoMuAZM (Int_t itrk)

Protected Attributes

TMultiLayerPerceptron * fneural
TFile * fLikeliFile
TH1F * fLikeliHist [6]

Detailed Description

Definition at line 14 of file MadEdAnalysis.h.


Constructor & Destructor Documentation

MadEdAnalysis::MadEdAnalysis ( TChain *  chainSR = 0,
TChain *  chainMC = 0,
TChain *  chainTH = 0,
TChain *  chainEM = 0 
)

Definition at line 15 of file MadEdAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, MadBase::InitChain(), MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.

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 }

MadEdAnalysis::MadEdAnalysis ( JobC j,
string  path,
int  entries 
)

Definition at line 40 of file MadEdAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, MadBase::fChain, MadBase::fJC, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.

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 }

MadEdAnalysis::~MadEdAnalysis (  ) 

Definition at line 60 of file MadEdAnalysis.cxx.

00061 {
00062 }


Member Function Documentation

void MadEdAnalysis::DataHist ( std::string  tag  ) 

Definition at line 2427 of file MadEdAnalysis.cxx.

References MadQuantities::CCNCSepVars(), NtpSRVertex::dcosz, NtpSRTrack::ds, NtpSREvent::end, ETrkF(), MadBase::eventSummary, EvtLength(), Evtphsig(), MadBase::GetEntry(), HitF(), MadQuantities::IsFidAll(), NtpSRTrack::lin, MadBase::LoadEvent(), MyLikeliQE(), NtpSRPlane::n, MadBase::Nentries, NeuNetEval(), NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpTrack, NtpSREvent::ntrack, MadQuantities::PassTrackCuts(), NtpSREvent::ph, NtpSREvent::plane, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), NtpSRPulseHeight::sigcor, NtpSREvent::trk, Trkphsig(), Trkplanes(), NtpSREvent::vtx, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

02427                                          {
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   //#declare lots of histos:
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   //old discrim histos:
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   //read in independent root file for use with likelihood analysis, etc.
02592   TFile *likeliFile = new TFile("PidQE_R1.9n.root","READ");
02593   TH1F *myhist[9];
02594   //myhist = NULL;
02595 
02596   if(likeliFile) { //if file is found
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   //TFile *f = new TFile("NeuNetnear.root");
02622   //TMultiLayerPerceptron *neural=this->NeuNetTrain(f);
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     //if(eventSummary->nevent!=1) continue;
02635 
02636    
02637     
02638     if(eventSummary->nevent==0) continue;
02639 
02640 
02641     for(int i=0;i<eventSummary->nevent;i++){
02642 
02643  
02644        //fiducial cuts on vtx for far detector
02645 
02646        // if(ntpEvent->vtx.z<0.5||ntpEvent->vtx.z>29.4
02647        //||(ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5)
02648        //||fabs(ntpEvent->vtx.x)>3.5||fabs(ntpEvent->vtx.y)>3.5) continue;
02649  
02650        // if(ntpEvent->vtx.z<0.5||ntpEvent->vtx.z>29.4
02651        //  ||(ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5)
02652        //  ||fabs(ntpEvent->vtx.x)>3.5||fabs(ntpEvent->vtx.y)>3.5) fdcut=1;
02653 
02654 
02655       //fiducial cuts on vtx near detector
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; //no event found
02662    
02663       
02664 
02665         //if(ntpEvent->ntrack!=1) continue;
02666 
02667         int *tracks = ntpEvent->trk;  //get array of track indices from this event
02668  
02669         EvtSum=eventSummary->nevent;
02670 
02671         for(int j=0;j<ntpEvent->ntrack;j++){  //loop over the tracks
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             //Use the NN
02736 
02737              nnval=0;
02738              nnval=NeuNetEval(0);
02739 
02740              // cout << nnval << endl;
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         } //loop over tracks
02809 
02810      
02811     }  //loop over eventsummaries
02812 
02813 
02814   } // end of entries
02815   
02816   tree->Write();
02817   save->Write();
02818   save->Close();
02819 
02820 
02821 }

Float_t MadEdAnalysis::DeDx ( Int_t  itr  )  [protected]

Definition at line 1265 of file MadEdAnalysis.cxx.

References NtpSRTrack::ds, MadBase::LoadTrack(), MadBase::ntpTrack, NtpSRTrack::ph, and NtpSRPulseHeight::sigcor.

Referenced by MCHist().

01265                                     {
01266 
01267   float ntpED = ntpTrack->ds;
01268   if(LoadTrack(itr)) return ntpTrack->ph.sigcor/(100.*ntpED);
01269   return 0.;
01270 }

Float_t MadEdAnalysis::ETrkF ( Int_t  itr  )  [protected]

Definition at line 1254 of file MadEdAnalysis.cxx.

References MadBase::LoadTrack(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSRTrack::ph, NtpSREvent::ph, and NtpSRPulseHeight::sigcor.

Referenced by DataHist(), and MCHist().

01254                                      {
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 }

Float_t MadEdAnalysis::EvtLength (  )  [protected]

Definition at line 1219 of file MadEdAnalysis.cxx.

References MadBase::LoadEvent(), NtpSRPlane::n, MadBase::ntpEvent, and NtpSREvent::plane.

Referenced by DataHist(), LikeliPID(), MCHist(), NeuNetEval(), and PIDHisto().

01219                                   {
01220   if(LoadEvent(0)) return ntpEvent->plane.n;
01221   return 0;
01222 }

Float_t MadEdAnalysis::Evtphsig (  )  [protected]

Definition at line 1237 of file MadEdAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::ntpEvent, NtpSREvent::ph, and NtpSRPulseHeight::sigcor.

Referenced by DataHist(), LikeliPID(), MCHist(), NeuNetEval(), and PIDHisto().

01237                                  {
01238   if(LoadEvent(0)) return ntpEvent->ph.sigcor;
01239   return 0;  
01240 }

Float_t MadEdAnalysis::HitF ( Int_t  itr  )  [protected]

Definition at line 1243 of file MadEdAnalysis.cxx.

References MadBase::LoadTrack(), NtpSREvent::nstrip, NtpSRTrack::nstrip, MadBase::ntpEvent, and MadBase::ntpTrack.

Referenced by DataHist(), and MCHist().

01243                                     {
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 }

Float_t * MadEdAnalysis::LikeliPID (  ) 

Definition at line 1411 of file MadEdAnalysis.cxx.

References MuELoss::e, EvtLength(), Evtphsig(), Trkphsig(), and Trkplanes().

Referenced by MakeEff().

01411                                   {
01412 
01413   TFile pfile("Pid_R1.9n.root","READ"); 
01414   
01415 
01416 
01417   //get detector type:
01418   //if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01419   //  pfile("Pid_R1.9f.root","READ");
01420   //}
01421   //else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01422   // pfile("Pid_R1.9n.root","READ");
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;  // ensure no floating errors 
01482   if (P_nc==0) P_nc=1e-04; // ensure no floating errors 
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;  // Pidvalsk (super k)
01490   PD[1] = P_cc/(P_cc+P_nc);  // Pidval2
01491  
01492   return Pid;
01493   
01494 }

void MadEdAnalysis::MakeEff ( std::string  tag  ) 

Definition at line 2141 of file MadEdAnalysis.cxx.

References NtpSREvent::end, MadBase::eventSummary, MadBase::GetEntry(), MadQuantities::IAction(), MadQuantities::INu(), MadQuantities::IResonance(), LikeliPID(), MadBase::LoadEvent(), MadBase::LoadTruthForRecoTH(), MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpEvent, NtpSREvent::ntrack, MadQuantities::PassCuts(), NtpSREvent::trk, MadQuantities::TrueNuEnergy(), NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

02141                                         {
02142 
02143   std::string savename = "EffHist_" + tag + ".root";
02144   TFile save(savename.c_str(),"RECREATE"); 
02145   save.cd();
02146   
02147   //#declare lots of histos:
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   //  TProfile *MuEffC = new TProfile("MuEffC","Events identified as CC",200,0.,20.,0.,1.);
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   // Call Pidhisto
02173   //PIDHisto();
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; //no event found
02195     
02196    
02197 
02198       int mc_entry = 0;
02199       if(!LoadTruthForRecoTH(i,mc_entry)) continue; //no mc event found
02200  
02201     
02202    
02203       //if you get to here, then both event and truth records are present
02204       int inu = INu(mc_entry);         //get some truth info for
02205       int iaction = IAction(mc_entry); //the appropriate MC event
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       //Fill Event && QE Efficiency
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         } //END OF QE  
02238 
02239 
02240       } //end of CC
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;  //get array of track indices from this event
02258      
02259 
02260 
02261 
02262 
02263 
02264 
02265 
02266 
02267      
02268       if(abs(inu)==14&&iaction==1) {
02269    
02270            
02271         //Fill track efficiency
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           //Fill single track effciency
02281           Int_t trks=0;
02282           if (ntpEvent->ntrack==1){ //select single tracks
02283                 
02284             trks=1;
02285             for(int k=0;k<ntpEvent->ntrack;k++){  //loop over the tracks
02286               int index = tracks[k];    
02287               MuEffS->Fill(NuEn,trks);
02288              
02289 
02290                           
02291               // Fill Good track efficiency
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               //fiducial cuts on vtx
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             } //end of for loop
02327           }
02328           else { MuEffS->Fill(NuEn,trks);}
02329 
02330 
02331         } //end of track efficiency
02332         else {
02333           MuEff->Fill(NuEn,trk);}
02334  
02335 
02336 
02337 
02338 
02339       
02340       } //end of true nu
02341 
02342 
02343       
02344     
02345 
02346 
02347       //fiducial cuts on vtx far detector
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       //fiducial cuts on vtx near detector
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++){  //loop over the tracks
02369         int index = tracks[j];             //get the index
02370 
02371 
02372 
02373 
02374 
02375         if(abs(inu)==14&&iaction==1) {
02376                  
02377 
02378         
02379                  
02380 
02381         //Selection Efficiency all cuts
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; //All cuts
02388               
02389           MuEffA->Fill(NuEn,muon_countall);
02390         }
02391       
02392 
02393 else { MuEffA->Fill(NuEn,muon_countall);}
02394 
02395 
02396 
02397 
02398 
02399  
02400           } //true CC
02401 
02402        
02403       } // End of tracks loop
02404      
02405   
02406 
02407     } //End of eventSummary loop
02408 
02409 
02410 
02411 
02412   }
02413   
02414   save.Write();
02415   save.Close();
02416 }

void MadEdAnalysis::MCHist ( std::string  tag  ) 

Definition at line 2824 of file MadEdAnalysis.cxx.

References MadQuantities::CCNCSepVars(), NtpSRVertex::dcosz, DeDx(), NtpSRTrack::ds, NtpSREvent::end, ETrkF(), MadBase::eventSummary, EvtLength(), Evtphsig(), MadBase::GetEntry(), HitF(), MadQuantities::IAction(), MadQuantities::INu(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), NtpSRTrack::lin, MadBase::LoadEvent(), MadBase::LoadTruthForRecoTH(), MyLikeliQE(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpTrack, MadBase::ntpTruth, NtpSREvent::ntrack, MadQuantities::PassTrackCuts(), NtpSREvent::ph, NtpSREvent::plane, MadQuantities::Q2(), MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), NtpSRPulseHeight::sigcor, NtpSREvent::trk, Trkphsig(), Trkplanes(), MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueShwEnergy(), NtpSREvent::vtx, NtpSRTrack::vtx, MadQuantities::W2(), NtpSRVertex::x, NtpMCTruth::y, NtpSRVertex::y, MadQuantities::Y(), and NtpSRVertex::z.

02824                                        {
02825 
02826  
02827 
02828   //True Quantities
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   //Reco Quantities
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 //#declare lots of histos:
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   //old discrim histos:
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    //read in independent root file for use with likelihood analysis, etc.
03216    TFile *likeliFile = new TFile("PidQE_R1.9n.root","READ");
03217    TH1F *myhist[9];
03218    //myhist = NULL;
03219 
03220    if(likeliFile) { //if file is found
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    //TFile *f = new TFile("NeuNetnear.root");
03244    //TMultiLayerPerceptron *neural=this->NeuNetTrain(f);
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      // if(eventSummary->nevent!=1) continue;
03256 
03257   
03258    
03259     
03260      if(eventSummary->nevent==0) continue;
03261 
03262 
03263      for(int i=0;i<eventSummary->nevent;i++){
03264 
03265 
03266        //fiducial cuts on vtx for far detector
03267 
03268        // if(ntpEvent->vtx.z<0.5||ntpEvent->vtx.z>29.4
03269        //||(ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5)
03270        //||fabs(ntpEvent->vtx.x)>3.5||fabs(ntpEvent->vtx.y)>3.5) continue;
03271  
03272        // if(ntpEvent->vtx.z<0.5||ntpEvent->vtx.z>29.4
03273        //  ||(ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5)
03274        //  ||fabs(ntpEvent->vtx.x)>3.5||fabs(ntpEvent->vtx.y)>3.5) fdcut=1;
03275 
03276 
03277       //fiducial cuts on vtx near detector
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; //no event found
03290          int mc_entry = 0;
03291 
03292     
03293          if(!LoadTruthForRecoTH(i,mc_entry)) continue;//no mc event found
03294 
03295          //if you get to here, then both event and truth records are present
03296          int inu = this->INu(mc_entry);         //get some truth info for
03297          int iaction = this->IAction(mc_entry); //the appropriate MC event
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;  //get array of track indices from this event
03317  
03318          EvtSum=eventSummary->nevent;
03319 
03320          for(int j=0;j<ntpEvent->ntrack;j++){  //loop over the tracks
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              //Use the NN
03492 
03493              nnval=0;
03494              // nnval=NeuNetEval(0);
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              //            cout <<mu_angzt<< " " << mu_angz<< " "<< cosz << " " <<ntpTruth->p4mu1[0] << " " << ntpTruth->p4mu1[1] << "  " << ntpTruth->p4mu1[2]<< "  " << ntpTruth->p4mu1[3] << "  " << NuEn << " " << " " << ntpTruth->p4shw[3]<< " " << RecoNuEn << " " << ntpHeader->GetSnarl() << endl;
03586        
03587           
03588 
03589            }
03590 
03591      
03592          } //loop over tracks
03593 
03594       
03595      }  //loop over eventsummaries
03596 
03597 
03598    } // end of entries
03599   
03600    tree->Write();
03601    save->Write();
03602    save->Close();
03603 }

void MadEdAnalysis::MyCreatePAN ( std::string  tag  ) 

Definition at line 246 of file MadEdAnalysis.cxx.

References NtpSRTrack::end, NtpSRMomentum::eqp, MadBase::eventSummary, MadQuantities::Flavour(), VldContext::GetDetector(), MadBase::GetEntry(), GnumiInterface::GetParent(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), RecHeader::GetVldContext(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), MadQuantities::Initial_state(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), MadBase::isST, Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), MadBase::mcrecord, NtpSRTrack::momentum, MyIsFidVtxrz(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, NtpSREvent::ntrack, MadQuantities::Nucleus(), PassAnalysisCuts(), PassBasicCuts(), NtpSREventSummary::ph, NtpSRShower::ph, NtpSRTrack::ph, PID(), NtpSRTrack::plane, NtpSREvent::plane, MadQuantities::Q2(), NtpSRMomentum::qp, NtpSRMomentum::range, RecoMuAZM(), MadQuantities::RecoMuDCosNeu(), RecoMuDCosY(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), MadQuantities::RecoShwEnergySqrt(), run(), NtpSRPulseHeight::sigcor, GnumiInterface::Status(), MadQuantities::Target4Vector(), MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSREvent::vtx, NtpSRTrack::vtx, MadQuantities::W2(), NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::Y(), NtpSRVertex::z, and NuParent::Zero().

00246                                             {
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   //PAN Quantities 
00269   //Truth:
00270   Float_t true_enu;       // true neutrino energy (GeV)
00271   Float_t true_pxnu;      // true neutrino momentum-x (GeV)
00272   Float_t true_pynu;      // true neutrino momentum-y (GeV)
00273   Float_t true_pznu;      // true neutrino momentum-z (GeV)
00274   Float_t true_etgt;       // true target energy (GeV)
00275   Float_t true_pxtgt;      // true target momentum-x (GeV)
00276   Float_t true_pytgt;      // true target momentum-y (GeV)
00277   Float_t true_pztgt;      // true target momentum-z (GeV)
00278   Float_t true_emu;       // true muon energy
00279   Float_t true_eshw;      // true shower energy
00280   Float_t true_x;         // true x
00281   Float_t true_y;         // true y
00282   Float_t true_q2;        // true q2
00283   Float_t true_w2;        // true w2
00284   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
00285   Float_t true_dircosz;   // track z-direction cosine
00286   Float_t true_vtxx;      // true x vtx
00287   Float_t true_vtxy;      // true y vtx
00288   Float_t true_vtxz;      // true z vtx
00289 
00290   //For Neugen:
00291   Int_t flavour;        // true flavour: 1-e 2-mu 3-tau
00292   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
00293   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
00294   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
00295   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
00296   Int_t had_fs;           // hadronic final state, number between 200-300
00297   
00298   //For Beam Reweighting:
00299   NuParent *nuparent = new NuParent();
00300 
00301   //Reco Quantities
00302   Int_t detector;         // Near = 1, Far = 2;
00303   Int_t run;              // run number
00304   Int_t snarl;            // snarl number
00305   Int_t nevent;           // number of events in snarl
00306   Int_t event;            // event index
00307   Int_t mcevent;          // mc event index
00308   Int_t ntrack;           // number of tracks in event
00309   Int_t nshower;          // number of showers in event
00310 
00311   Int_t is_fid;           // pass fid cut. true = 1, false = 0
00312   Int_t is_fidrk;         // pass fid cut for rock muons. true = 1, false = 0
00313                           // is_fidrk=1 (true beam), is_fidrk=0 (rock muon)
00314 
00315   Int_t is_pcev;          // partially contained event
00316   Int_t is_cev;           // fully contained. true = 1, false = 0 
00317   Int_t is_mc;            // is a corresponding mc event found
00318   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
00319   Float_t pid0;           // pid parameter (usual method)
00320   Float_t pid1;           // pid parameter (alternative method 1)
00321   Float_t pid2;           // pid parameter (alternative method 2)
00322   Float_t Pidvalsk;       // Super K parameter
00323 
00324   Int_t pass;             // pass analysis cuts. true = 1, false = 0
00325 
00326   Float_t reco_enu;       // reco neutrino energy
00327   Float_t reco_emu;       // reco muon energy
00328   Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
00329   Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
00330   Float_t reco_qe_enu;    // reco qe neutrino energy
00331   Float_t reco_qe_q2;     // reco qe q2
00332   Float_t reco_y;         // reco y
00333   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
00334   Float_t reco_dircosz;   // reco track vtx z-dircos
00335   Float_t reco_dircosy;   // = -(zenith)
00336   Float_t reco_azimuth;   // = azimuth
00337   Float_t reco_vtxx;      // reco x vtx
00338   Float_t reco_vtxy;      // reco y vtx
00339   Float_t reco_vtxz;      // reco z vtx
00340 
00341   Float_t evtlength;      // reco event length
00342   Float_t trklength;      // reco track length
00343   Float_t trkmomrange;    // reco track momentum from range
00344   Float_t trkqp;          // reco track q/p
00345   Float_t trkeqp;         // reco track q/p error
00346   Float_t trkphfrac;      // reco pulse height fraction in track
00347   Float_t trkphplane;     // reco average track pulse height/plane
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   //Trees
00361   TTree *tree = new TTree("pan","pan");
00362   //Truth
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   //Neugen
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   //NuParent
00390   tree->Branch("nuparent","NuParent",&nuparent,8000,1);
00391   //Reco
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     //fill tree once for each reconstructed event
00465     for(int i=0;i<nevent;i++){ 
00466       if(!LoadEvent(i)) continue; //no event found
00467 
00468       totph = 0;
00469 
00470       //set event index:
00471       event = i;
00472       ntrack = ntpEvent->ntrack;
00473       nshower = ntpEvent->nshower;
00474 
00475       //zero all tree values
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       // Pidvalsk = 0;
00497       //get detector type:
00498       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00499         detector = 2;
00500         totph=eventSummary->ph.sigcor;
00501         //fiducial criteria on vtx for far detector
00502         is_fid = 1;
00503         if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
00504            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00505            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00506                 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00507            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
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         //get total charge in 
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         //fiducial criteria on vtx for near detector
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       //fiducial criteria on vtx for rock muons
00534       //temporary fix
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       // cout << MyIsFidVtx(event) << " " << is_fidrk << endl;
00560 
00561       }
00562       
00563       //check for a corresponding mc event      
00564       if(LoadTruthForRecoTH(i,mcevent)) {
00565         Float_t *vtx = TrueVtx(mcevent);
00566         Float_t *nu3mom = TrueNuMom(mcevent);
00567         Float_t *targ4vec = Target4Vector(mcevent);
00568         //true info for tree:
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) ;//gnumi.GetParent(strecord,mcevent,*nuparent);
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         //reco info for tree:
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         //index will -1 if no track/shower present in event
00624         int track_index   = -1;
00625         LoadLargestTrackFromEvent(i,track_index);
00626         int shower_index  = -1;
00627         LoadLargestShowerFromEvent(i,shower_index);
00628 
00629         //methods should all be protected against index = -1
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){ //check track is present before using ntpTrack
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         // cout << reco_dircosz << " " << RecoMuDCosY(track_index) << " " <<  RecoMuZn(track_index) << " " << RecoMuAZM(track_index) << endl; 
00673 
00674 
00675 
00676 
00677 
00678         if(PassAnalysisCuts(event)) pass = 1;
00679       }//end of PassBasicCuts
00680 
00681       //select neutrino beam events
00682       if(is_fid==1&&is_pcev==1&&pass_basic==1){
00683         ofs << ntpHeader->GetSnarl() << endl;
00684       }
00685 
00686 
00687 
00688       //fill the tree
00689       tree->Fill();
00690     }
00691   }
00692   delete nuparent;
00693   
00694   save.cd();
00695   tree->Write();
00696   save.Write();
00697   save.Close();
00698 }

void MadEdAnalysis::MyCreatePANData ( std::string  tag  ) 

Definition at line 708 of file MadEdAnalysis.cxx.

References NtpSRTrack::end, NtpSRMomentum::eqp, MadBase::eventSummary, VldContext::GetDetector(), MadBase::GetEntry(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), RecHeader::GetVldContext(), MadQuantities::IsFidAll(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadTrack(), NtpSRTrack::momentum, NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, NtpSREvent::ntrack, PassAnalysisCuts(), PassBasicCuts(), NtpSREventSummary::ph, NtpSRShower::ph, NtpSRTrack::ph, PID(), NtpSRTrack::plane, NtpSREvent::plane, NtpSRMomentum::qp, NtpSRMomentum::range, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), run(), NtpSRPulseHeight::sigcor, GnumiInterface::Status(), NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

00708                                                 {
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   //PAN Quantities 
00726 
00727 
00728   //Reco Quantities
00729   Int_t detector;         // Near = 1, Far = 2;
00730   Int_t run;              // run number
00731   Int_t snarl;            // snarl number
00732   Int_t nevent;           // number of events in snarl
00733   Int_t event;            // event index
00734   Int_t mcevent;          // mc event index
00735   Int_t ntrack;           // number of tracks in event
00736   Int_t nshower;          // number of showers in event
00737 
00738   Int_t is_fid;           // pass fid cut. true = 1, false = 0
00739   Int_t is_fidrk;         // pass fid cut for rock muons. true = 1, false = 0
00740                           // is_fidrk=1(true beam), is_fidrk=0(rock muon)
00741   Int_t is_cev;           // fully contained. true = 1, false = 0 
00742   Int_t is_mc;            // is a corresponding mc event found
00743   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
00744   Float_t pid0;           // pid parameter (usual method)
00745   Float_t pid1;           // pid parameter (alternative method 1)
00746   Float_t pid2;           // pid parameter (alternative method 2)
00747   Float_t Pidvalsk;       // Super K parameter
00748 
00749   Int_t pass;             // pass analysis cuts. true = 1, false = 0
00750 
00751   Float_t reco_enu;       // reco neutrino energy
00752   Float_t reco_emu;       // reco muon energy
00753   Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
00754   Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
00755   Float_t reco_qe_enu;    // reco qe neutrino energy
00756   Float_t reco_qe_q2;     // reco qe q2
00757   Float_t reco_y;         // reco y
00758   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
00759   Float_t reco_dircosz;   // reco track vtx z-dircos
00760   Float_t reco_vtxx;      // reco x vtx
00761   Float_t reco_vtxy;      // reco y vtx
00762   Float_t reco_vtxz;      // reco z vtx
00763 
00764   Float_t evtlength;      // reco event length
00765   Float_t trklength;      // reco track length
00766   Float_t trkmomrange;    // reco track momentum from range
00767   Float_t trkqp;          // reco track q/p
00768   Float_t trkeqp;         // reco track q/p error
00769   Float_t trkphfrac;      // reco pulse height fraction in track
00770   Float_t trkphplane;     // reco average track pulse height/plane
00771 
00772   //Trees
00773   TTree *tree = new TTree("pan","pan");
00774  
00775   //Reco
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   //tree->Branch("pid1",&pid1,"pid1/F",32000);
00793   //tree->Branch("pid2",&pid2,"pid2/F",32000);
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     //fill tree once for each reconstructed event
00834     for(int i=0;i<nevent;i++){ 
00835       if(!LoadEvent(i)) continue; //no event found
00836 
00837       totph = 0;
00838 
00839       //set event index:
00840       event = i;
00841       ntrack = ntpEvent->ntrack;
00842       nshower = ntpEvent->nshower;
00843 
00844       //zero all tree values
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       // Pidvalsk = 0;
00856       //get detector type:
00857       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00858         detector = 2;
00859         totph=eventSummary->ph.sigcor;
00860         //fiducial criteria on vtx for far detector
00861         is_fid = 1;
00862         if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
00863            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00864            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00865                 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00866            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
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         //get total charge in 
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         //fiducial criteria on vtx for near detector
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       //fiducial criteria on vtx for rock muons
00893       //temporary fix
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       // cout << MyIsFidVtx(event) << " " << is_fidrk << endl;
00912 
00913       }
00914       
00915    
00916 
00917  
00918 
00919 
00920       if(PassBasicCuts(event)) {
00921         //reco info for tree:
00922 
00923         // cout << flavour << endl;       
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         //index will -1 if no track/shower present in event
00932         int track_index   = -1;
00933         LoadLargestTrackFromEvent(i,track_index);
00934         int shower_index  = -1;
00935         LoadLargestShowerFromEvent(i,shower_index);
00936 
00937         //methods should all be protected against index = -1
00938         reco_emu          = RecoMuEnergy(0,track_index);
00939         reco_eshw         = RecoShwEnergy(shower_index);
00940         //reco_eshw_sqrt    = RecoShwEnergySqrt(shower_index);
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){ //check track is present before using ntpTrack
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       //fill the tree
00979       tree->Fill();
00980     }
00981   }
00982    
00983   save.cd();
00984   tree->Write();
00985   save.Write();
00986   save.Close();
00987 }

Bool_t MadEdAnalysis::MyIsFidVtx ( Int_t  itrk = 0  ) 

Definition at line 67 of file MadEdAnalysis.cxx.

References NtpSRTrack::end, MadBase::LoadTrack(), MadBase::ntpTrack, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by PID().

00067                                           {
00068  if(!LoadTrack(itrk)) return false;
00069   //if(ntpTrack->fidvtx.dr<0.5||ntpTrack->fidvtx.dz<0.5) return false;
00070 
00071 
00072 
00073   //temporary fix
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 }

Bool_t MadEdAnalysis::MyIsFidVtxrz ( Int_t  itrk = 0  ) 

Definition at line 94 of file MadEdAnalysis.cxx.

References NtpSRFiducial::dr, NtpSRFiducial::dz, NtpSRTrack::fidvtx, MadBase::LoadTrack(), and MadBase::ntpTrack.

Referenced by MyCreatePAN().

00094                                             {
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 }

Float_t * MadEdAnalysis::MyLikeliQE ( TH1F **  hist  ) 

Definition at line 1497 of file MadEdAnalysis.cxx.

References MadQuantities::CCNCSepVars(), NtpSRTrack::ds, NtpSREvent::nstrip, NtpSRTrack::nstrip, MadBase::ntpEvent, MadBase::ntpTrack, NtpSRTrack::ph, and NtpSRPulseHeight::sigcor.

Referenced by DataHist(), and MCHist().

01497                                              {
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   //NC:      
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   //std::cout << ProbNC << std::endl;
01527   
01528   //CC:
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   //std::cout << ProbCC << std::endl;
01539 
01540   //QE:
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   //pid pars
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 }

void MadEdAnalysis::MyMakeMyFile ( std::string  tag  ) 

Definition at line 992 of file MadEdAnalysis.cxx.

References MadBase::eventSummary, NtpSRTrack::fit, MadQuantities::Flavour(), VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), MadQuantities::IAction(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadTHTrack(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRFitTrack::pass, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

00992                                              {
00993 
00994   std::string savename = "DPHistos_" + tag + ".root";
00995   TFile save(savename.c_str(),"RECREATE"); 
00996   save.cd();
00997   
00998   //#declare lots of histos:
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     //cout << "\n\n**ENTRY=" << i << endl;
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     //Int_t run= ntpHeader->GetRun();
01046     //Int_t snarl = ntpHeader->GetSnarl();
01047   
01048 
01049     if(nevent==0) continue;
01050     
01051     Int_t trkcount=0;
01052     Int_t shwcount=0;
01053 
01054         
01055     //fill tree once for each reconstructed event
01056     for(int j=0;j<eventSummary->nevent;j++){
01057       //      cout << "*EVENT=" << j << endl;
01058       if(!LoadEvent(j)) continue; //no event found
01059 
01060       //      ntpEvent->Print();
01061       //set event index:
01062       //Int_t event = j;
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       //      cout << "Event ph=" << ntpEvent->ph.sigcor << endl;
01071       //      cout << "TRACKS:" << endl;
01072       for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
01073         LoadTrack(ii);
01074         //      cout << " track #" << ii << " ph=" << ntpTrack->ph.sigcor << endl;
01075         totph+=ntpTrack->ph.sigcor;
01076         LoadTHTrack(ii);
01077         //cout << "track=" << ii << " neumc=" << ntpTHTrack->neumc << endl;
01078       }
01079 
01080       trkcount+=ntrack;
01081 
01082       //       cout << "SHOWERS:" << endl;
01083        for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
01084         LoadShower(ii);
01085         //      cout << " shower #" << ii << " ph=" << ntpShower->ph.sigcor << endl;
01086         totph+=ntpShower->ph.sigcor;
01087         }
01088 
01089        shwcount+=nshower;
01090  
01091        // cout << "totph=" << totph << endl;
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       //get detector type:
01100       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01101         detector = 2;
01102         //fiducial criteria on vtx for far detector
01103         is_fid = 1;
01104         if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
01105            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
01106            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
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         //fiducial criteria on vtx for near detector
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       //check for a corresponding mc event
01119  
01120       Int_t mcevent=-1;
01121 
01122       if(LoadTruthForRecoTH(j,mcevent)) {
01123         //true info for tree:
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         //      cout << "track index=" << track_index << endl;
01136 
01137         Int_t trkpass=ntpTrack->fit.pass;
01138 
01139 
01140         if (is_fid && trkpass) {
01141           
01142           //    cout << "pass basic\n";
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         //      cout << "phfrac=" << phfrac << endl;
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           //cout << "evlength:phfrac:phplane=" << evlength << ":"
01167           //   << phfrac << ":" << phplane << endl;
01168           //cout << "event:mcevent" << j << ":" << mcevent << endl;
01169           //cout << "run:snarl" << run << ":" << snarl << endl;
01170           evlength_nc->Fill(evlength);
01171           phfrac_nc->Fill(phfrac);
01172           phplane_nc->Fill(phplane);     
01173         }
01174 
01175 
01176         }
01177       }
01178     }
01179   }
01180   //#close file  
01181   save.Write();
01182   save.Close();
01183 }

void MadEdAnalysis::MyMakeQEFile ( std::string  tag  ) 

Definition at line 1586 of file MadEdAnalysis.cxx.

References MadQuantities::CCNCSepVars(), NtpSRTrack::ds, NtpSREvent::end, MadBase::eventSummary, VldContext::GetDetector(), MadBase::GetEntry(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), RecHeader::GetVldContext(), MadQuantities::IAction(), MadQuantities::INu(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadTruthForRecoTH(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSREvent::nstrip, NtpSRTrack::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpTrack, NtpSREvent::ntrack, PassBasicCuts(), MadQuantities::PassTrackCuts(), NtpSRTrack::ph, NtpSREvent::ph, NtpSREvent::plane, MadQuantities::Q2(), MadQuantities::RecoMuEnergy(), run(), NtpSRPulseHeight::sigcor, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSREvent::vtx, MadQuantities::W2(), NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::Y(), and NtpSRVertex::z.

01586                                              {
01587 
01588   std::string savename = "PidQE_" + tag + ".root";
01589   TFile save(savename.c_str(),"RECREATE"); 
01590   save.cd();
01591   
01592   //#declare lots of histos:
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   //old discrim histos:
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   //Truth:
01834   Float_t true_enu;       // true neutrino energy (GeV)
01835   Float_t true_emu;       // true muon energy
01836   Float_t true_eshw;      // true shower energy
01837   Float_t true_x;         // true x
01838   Float_t true_y;         // true y
01839   Float_t true_q2;        // true q2
01840   Float_t true_w2;        // true w2
01841   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
01842   Float_t true_dircosz;   // track z-direction cosine
01843   Float_t true_vtxx;      // true x vtx
01844   Float_t true_vtxy;      // true y vtx
01845   Float_t true_vtxz;      // true z vtx
01846   Int_t process = 0;
01847 
01848   //Reco Quantities
01849   Int_t detector;         // Near = 1, Far = 2;
01850   Int_t run;              // run number
01851   Int_t snarl;            // snarl number
01852   Int_t nevent;           // number of events in snarl
01853   Int_t event;            // event index
01854   Int_t mcevent;          // mc event index
01855   Int_t ntrack;           // number of tracks in event
01856   Int_t nshower;          // number of showers in event
01857   Int_t inu = 0;
01858   Int_t iaction = 0;
01859   Int_t is_fid = 0;       // pass fid cut. true = 1, false = 0
01860   //Int_t is_cev;           // fully contained. true = 1, false = 0 
01861   Int_t is_mc;            // is a corresponding mc event found
01862   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
01863 
01864 
01865   //Int_t pass;             // pass analysis cuts. true = 1, false = 0
01866 
01867   //Float_t reco_enu;       // reco neutrino energy
01868   Float_t reco_emu = 0.;    // reco muon energy
01869   Float_t reco_eshw = 0.;   // reco shower energy  (shw.ph.gev/1.23)
01870   //Float_t reco_qe_enu;    // reco qe neutrino energy
01871   //Float_t reco_qe_q2;     // reco qe q2
01872   //Float_t reco_y;         // reco y
01873   //Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
01874   Float_t reco_dircosz = 0.; // reco track vtx z-dircos
01875   //Float_t reco_vtxx;      // reco x vtx
01876   //Float_t reco_vtxy;      // reco y vtx
01877   //Float_t reco_vtxz;      // reco z vtx
01878   //Float_t fitpass;        // track fit pass 
01879   //Float_t evtlength;      // reco event length
01880   //Float_t trklength;      // reco track length
01881   //Float_t trkmomrange;    // reco track momentum from range
01882   //Float_t trkqp;          // reco track q/p
01883   //Float_t trkeqp;         // reco track q/p error
01884   //Float_t trkphfrac;      // reco pulse height fraction in track
01885   //Float_t trkphplane;     // reco average track pulse height/plane
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     //fill tree once for each reconstructed event
01912     for(int i=0;i<eventSummary->nevent;i++){ 
01913       if(!LoadEvent(i)) continue; //no event found
01914 
01915       //set event index:
01916       event = i;
01917       ntrack = ntpEvent->ntrack;
01918       nshower = ntpEvent->nshower;
01919 
01920       //zero all tree values
01921      
01922       detector = -1; mcevent = -1;
01923     
01924 
01925       //get detector type:
01926       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01927         detector = 2;
01928         //fiducial criteria on vtx for far detector
01929         is_fid = 1;
01930         if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
01931            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
01932            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
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         //fiducial criteria on vtx for near detector
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     //check for a corresponding mc event
01946       if(LoadTruthForRecoTH(i,mcevent)) {
01947 
01948         mcindex           = mcevent;
01949         Float_t *vtx = TrueVtx(mcevent);
01950         //true info for tree:
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         //reco info for tree:
01974         recoindex         = event;
01975         pass_basic        = 1;
01976 
01977 
01978        
01979         //index will -1 if no track/shower present in event
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           //if(sqrt(true_w2)>1.0) {
01998           if(process!=1001) {
01999             NEvent_CC->Fill(eventSummary->nevent);
02000 
02001             //if(eventSummary->nevent==1){
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             //TrkMomFrac0_CC->Fill(reco_emu/(reco_emu+reco_eshw));
02016 
02017             if(reco_dircosz>0.) dsdz_CC->Fill(1./reco_dircosz);
02018 
02019             //   NHitTrkFrac_Vs_VtxShwEnergy_CC->Fill(reco_eshw,
02020             //                           float(ntpTrack->nstrip)
02021             //                           /float(ntpEvent->nstrip));
02022             // NHitTrkFrac_Vs_Y_CC->Fill(ntpTruth->y,float(ntpTrack->nstrip)
02023             //                /float(ntpEvent->nstrip));
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        } //process
02040       else {
02041         NEvent_QE->Fill(eventSummary->nevent); 
02042         //if(eventSummary->nevent==1){
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     } //end of true cc events
02077 
02078      
02079       if(iaction==0){
02080       NEvent_NC->Fill(eventSummary->nevent);
02081       //if(eventSummary->nevent==1){
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           //float trkenergy = this->RecoMuEnergy(0,0);
02092           //float shwenergy = this->RecoShwEnergy(1);
02093           //float shwenergy0 = this->RecoShwEnergy(0);
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      }//end of iaction==0
02115 
02116 
02117 
02118 
02119 
02120 
02121 
02122 
02123 
02124 
02125 
02126     
02127 
02128        } //if passbasiccut
02129 
02130     } //eventsummaryloop
02131   
02132   } //for entries
02133   save.Write();
02134   save.Close();
02135   
02136 }

void MadEdAnalysis::MyReadPIDFile ( std::string  tag  ) 

Definition at line 1188 of file MadEdAnalysis.cxx.

References fLikeliFile, and fLikeliHist.

01188                                               {
01189 
01190   fLikeliFile = new TFile(tag.c_str(),"READ");
01191  
01192   if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found
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 }

Float_t MadEdAnalysis::NeuNetEval ( Int_t  event = 0  ) 

Definition at line 1305 of file MadEdAnalysis.cxx.

References EvtLength(), Evtphsig(), fneural, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), Trkphsig(), and Trkplanes().

Referenced by DataHist(), and PID().

01305                                              {
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   //Use the NN
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 }

TMultiLayerPerceptron * MadEdAnalysis::NeuNetTrain ( TFile *  f  ) 

Definition at line 1273 of file MadEdAnalysis.cxx.

01273                                                           {
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   //training sample 500 evts
01283   treeN->Draw(">>elist1","flavor==1","",652,0);
01284   treeN->Draw(">>elist2","flavor==0","",2127,0);
01285   
01286   elist1->Add(elist2);
01287 
01288   //test sample 500 evts
01289   treeN->Draw(">>elist3","flavor==1","",659,652);
01290   treeN->Draw(">>elist4","flavor==0","",1656,2127);
01291   elist3->Add(elist4);
01292   //elist3->Print();
01293 
01294   //Train NN
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 }

Bool_t MadEdAnalysis::PassAnalysisCuts ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 109 of file MadEdAnalysis.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadEvent(), MadBase::ntpHeader, and PID().

Referenced by MyCreatePAN(), and MyCreatePANData().

00110 {
00111   if(!LoadEvent(event)) return false;
00112   //  if(PassBasicCuts(event)) return true;
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 }

Bool_t MadEdAnalysis::PassBasicCuts ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 121 of file MadEdAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, NtpSREvent::ntrack, and MadQuantities::PassTrackCuts().

Referenced by MyCreatePAN(), MyCreatePANData(), and MyMakeQEFile().

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   // if(!MyIsFidVtx()) return false;  
00130   return true;
00131 }

Float_t MadEdAnalysis::PID ( Int_t  event = 0,
Int_t  method = 0 
) [virtual]

Reimplemented from MadAnalysis.

Definition at line 168 of file MadEdAnalysis.cxx.

References MadBase::eventSummary, fLikeliHist, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MyIsFidVtx(), NtpSRPlane::n, NeuNetEval(), MadBase::ntpTrack, NtpSREventSummary::ph, NtpSRTrack::ph, NtpSRTrack::plane, NtpSREventSummary::plane, and NtpSRPulseHeight::sigcor.

Referenced by MyCreatePAN(), MyCreatePANData(), and PassAnalysisCuts().

00168                                                   {
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   // cout << "BP0\n";
00184   Float_t var1=eventSummary->plane.n;
00185   Float_t var2=0;
00186   Float_t var3=0;
00187 
00188   // cout << "BP1\n";
00189 
00190   if (eventSummary->ph.sigcor!=0) var2 = float(ntpTrack->ph.sigcor)
00191                                     /float(eventSummary->ph.sigcor);
00192   //  cout << "BP2\n";
00193 
00194   if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor)
00195                               /float(ntpTrack->plane.n);
00196   //  cout << "BP3\n";
00197   //  fLikeliHist[0]->Print();
00198 
00199   // cout << "var1=" << var1 << " var2=" << var2 << " var3=" << var3 << endl;
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   //  cout << "prob1=" << prob1 << " prob2=" << prob2 << " prob3=" << prob3 << endl;
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   // cout << "prob1=" << prob1 << " prob2=" << prob2 << " prob3=" << prob3 << endl;
00222   //  cout << "PROBCC=" << ProbCC << " PROBNC=" << ProbNC << endl;
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 }

void MadEdAnalysis::PIDHisto (  ) 

Definition at line 1331 of file MadEdAnalysis.cxx.

References MadBase::eventSummary, EvtLength(), Evtphsig(), MadBase::GetEntry(), MadQuantities::IAction(), MadQuantities::INu(), MadBase::LoadEvent(), MadBase::LoadTruthForRecoTH(), MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpEvent, NtpSREvent::ntrack, Trkphsig(), and Trkplanes().

01331                              {
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     // Fill and write out histograms for PID Calculations
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; //no event found
01368       int mc_entry = 0;
01369       if(!LoadTruthForRecoTH(i,mc_entry)) continue; //no mc event found
01370  
01371 
01372       //if you get to here, then both event and truth records are present
01373       int inu = INu(mc_entry);         //get some truth info for
01374       int iaction = IAction(mc_entry); //the appropriate MC event
01375 
01376 
01377       //int *tracks = ntpEvent->trk;  //get array of track indices from this event
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 }

Float_t MadEdAnalysis::RecoMuAZM ( Int_t  itrk  )  [protected]

Definition at line 156 of file MadEdAnalysis.cxx.

References NtpSRCosmicRay::azimuth, NtpSRTrack::cr, MadBase::LoadTrack(), and MadBase::ntpTrack.

Referenced by MyCreatePAN().

00156                                           { //azimuth
00157   if(!LoadTrack(itrk)) return 0.;
00158   return ntpTrack->cr.azimuth;
00159 }

Float_t MadEdAnalysis::RecoMuDCosY ( Int_t  itrk  )  [protected]

Definition at line 138 of file MadEdAnalysis.cxx.

References NtpSRVertex::dcosy, MadBase::LoadTrack(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by MyCreatePAN().

00138                                             { //dy/ds
00139   if(!LoadTrack(itrk)) return 0.;
00140   return ntpTrack->vtx.dcosy;
00141 }

Float_t MadEdAnalysis::RecoMuZn ( Int_t  itrk  )  [protected]

Definition at line 146 of file MadEdAnalysis.cxx.

References NtpSRTrack::cr, MadBase::LoadTrack(), MadBase::ntpTrack, and NtpSRCosmicRay::zenith.

00146                                          { //zenith
00147   if(!LoadTrack(itrk)) return 0.;
00148   return ntpTrack->cr.zenith;
00149 }

void MadEdAnalysis::SetMLP ( TMultiLayerPerceptron *  neural  )  [inline]

Definition at line 55 of file MadEdAnalysis.h.

References fneural.

00055 {fneural=neural;}

Float_t MadEdAnalysis::Trkphsig ( Int_t  itrk  )  [protected]

Definition at line 1225 of file MadEdAnalysis.cxx.

References MadBase::LoadTrack(), MadBase::ntpTrack, NtpSRTrack::ph, and NtpSRPulseHeight::sigcor.

Referenced by DataHist(), LikeliPID(), MCHist(), NeuNetEval(), and PIDHisto().

01225                                           {
01226   if(LoadTrack(itrk)) return ntpTrack->ph.sigcor;
01227   return 0.;
01228 }

Float_t MadEdAnalysis::Trkplanes ( Int_t  itrk  )  [protected]

Definition at line 1231 of file MadEdAnalysis.cxx.

References MadBase::LoadTrack(), NtpSRPlane::n, MadBase::ntpTrack, and NtpSRTrack::plane.

Referenced by DataHist(), LikeliPID(), MCHist(), NeuNetEval(), and PIDHisto().

01231                                            {
01232   if(LoadTrack(itrk)) return ntpTrack->plane.n;
01233   return 0.;
01234 }


Member Data Documentation

TFile* MadEdAnalysis::fLikeliFile [protected]

Definition at line 24 of file MadEdAnalysis.h.

Referenced by MyReadPIDFile().

TH1F* MadEdAnalysis::fLikeliHist[6] [protected]

Definition at line 25 of file MadEdAnalysis.h.

Referenced by MyReadPIDFile(), and PID().

TMultiLayerPerceptron* MadEdAnalysis::fneural [protected]

Definition at line 19 of file MadEdAnalysis.h.

Referenced by NeuNetEval(), and SetMLP().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1