MadAnalysis Class Reference

#include <MadAnalysis.h>

Inheritance diagram for MadAnalysis:

MadQuantities MadBase MadCBSQEAnalysis MadCluAnalysis MadDpAnalysis MadEdAnalysis MadMKAnalysis MadPIDAnalysis MadTestAnalysis MadTVAnalysis MadUserAnalysis List of all members.

Public Member Functions

 MadAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadAnalysis (JobC *, string, int)
virtual ~MadAnalysis ()
void ReInit (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
void ReInit (JobC *, string, int)
void CreatePAN (std::string)
void CreateANtpPAN (std::string)
void SetFileNameTag (std::string)
void SetHistBookInfo (int, float, float)
void RecoMCExperiment (int startnum, double NearPOT, double FarPOT)
void RecoExperiment (int startnum, double NearExpectedNormToPOT, double FarExpectedNormToPOT)
void RecoMC (int startnum, double NearPOT, double FarPOT)
void SetOscillationFunction (Double_t(*fcn)(Double_t *, Double_t *), Float_t, Float_t, int, double *)
void SetOscillationFunction (TF1 *, double *)
void VaryFitParam (int ix, int vx)
Bool_t Do2DFit (float n_a=90., float min_a=0.5, float max_a=1.0, float n_b=100., float min_b=0.0, float max_b=0.005, float n_c=1., float f_c=0.)
Bool_t Do2DFitNearFar (float n_a=90., float min_a=0.5, float max_a=1.0, float n_b=100., float min_b=0.0, float max_b=0.005, float n_c=1., float f_c=0.)
void SetNeugenReweightInfo (Int_t, Int_t *, Float_t *, Float_t *, Float_t *, Float_t *, std::string *)
MadChi2Calc * GetChi2Calc ()
void SetEShiftErr (float err)
int SetDataInfo (TH1F *hist=0, int numev=0, float pot=0, float sigfrac=0, int det=1)
TH1F * GetDataHist (int det)
int GetNumDataEvents (int det=1)
float GetDataPOT (int det=1)
float GetSignalFracInData (int det=1)
void NoOutFile ()
void EndJob ()

Protected Member Functions

virtual Bool_t PassTruthSignal (Int_t mcevent=0)
virtual Bool_t PassBasicCuts (Int_t event=0)
virtual Bool_t PassAnalysisCuts (Int_t event=0)
virtual Float_t PID (Int_t event=0, Int_t method=0)
virtual Float_t RecoAnalysisEnergy (Int_t event=0)
virtual Float_t GetWeight (Int_t event=0)
virtual Bool_t AddUserBranches (TTree *tree)
virtual void CallUserFunctions (Int_t)

Protected Attributes

int NumberOfBins
float RangeLowerLimit
float RangeUpperLimit
TH1F * MCSignalHist [2]
TH1F * MCBkgdHist [2]
TH1F * DataHist [2]
TH1F * DataSignalHist [2]
float DataPOT [2]
int numDataEvts [2]
float signalFracInData [2]
TH1F * DataBkgdHist [2]
vector< mcev_reweight > MCEvents [2]
float MCPOT [2]
int numMCEvts [2]
Int_t NgenBins [3]
Double_t NgenMin [3]
Double_t NgenMax [3]
Double_t NgenMean [3]
Double_t NgenErr [3]
std::string NgenName [3]
TF1 * OscillationFunction
int NumOscPars
double * OscillationParameters
int * varyX
TH1F * BestFit [2]
std::string tag
MadChi2Calc * Chi2Calc
TH2F * Chi2Surf
float ChiMin
float NDOF
float Par1MinVal
float Par2MinVal
float EShiftErr
Bool_t writeOut

Detailed Description

Definition at line 48 of file MadAnalysis.h.

Constructor & Destructor Documentation

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

Definition at line 36 of file MadAnalysis.cxx.

References BestFit, Chi2Calc, Chi2Surf, ChiMin, MadBase::Clear(), DataBkgdHist, DataHist, DataPOT, DataSignalHist, MadBase::emrecord, EShiftErr, MadBase::isEM, MadBase::isMC, MadBase::isTH, MCBkgdHist, MCEvents, MCPOT, MadBase::mcrecord, MCSignalHist, NDOF, MadBase::Nentries, NgenBins, NgenErr, NgenMax, NgenMean, NgenMin, NgenName, NumberOfBins, numDataEvts, numMCEvts, NumOscPars, OscillationFunction, OscillationParameters, Par1MinVal, Par2MinVal, RangeLowerLimit, RangeUpperLimit, MadBase::record, signalFracInData, MadBase::strecord, tag, MadBase::threcord, varyX, MadBase::whichSource, and writeOut.

00038    :MadQuantities(chainSR,chainMC,chainTH,chainEM)
00039 {
00040 
00041   NumberOfBins = 1000;
00042   RangeLowerLimit = 0.;
00043   RangeUpperLimit = 50.;
00044 
00045   MCSignalHist[0] = NULL;   MCSignalHist[1] = NULL;
00046   MCBkgdHist[0] = NULL;     MCBkgdHist[1] = NULL;
00047   DataHist[0] = NULL;       DataHist[1] = NULL;
00048   DataSignalHist[0] = NULL; DataSignalHist[1] = NULL;
00049   DataPOT[0] = 0;           DataPOT[1] = 0;
00050   numDataEvts[0] = 0;       numDataEvts[1] = 0;
00051   signalFracInData[0] = 0;  signalFracInData[1] = 0;
00052   DataBkgdHist[0] = NULL;   DataBkgdHist[1] = NULL;
00053   MCEvents[0].clear();      MCEvents[1].clear();
00054   MCPOT[0] = 0;             MCPOT[1] = 0;
00055   numMCEvts[0] = 0;         numMCEvts[1] = 0;
00056   BestFit[0] = NULL;        BestFit[1] = NULL;
00057 
00058   OscillationFunction = NULL;
00059   OscillationParameters = NULL;
00060   NumOscPars = 0;
00061   varyX = new int[2];
00062   varyX[0] = 1;
00063   varyX[1] = 2;
00064 
00065   NgenBins[0] = 5;   NgenBins[1] = 0;  NgenBins[2] = 0;
00066   NgenMin[0]  = 0.5; NgenMin[1]  = 0.; NgenMin[2] = 0.;
00067   NgenMax[0]  = 1.5; NgenMax[1]  = 0.; NgenMax[2] = 0.;
00068   NgenMean[0] = 1.;  NgenMean[1] = 0.; NgenMean[2] = 0.;
00069   NgenErr[0]  = 0.1; NgenErr[1]  = 0.; NgenErr[2] = 0.;
00070   NgenName[0] = "neugen:ma_qe";        NgenName[1] = "empty"; 
00071   NgenName[2] = "empty";
00072 
00073   tag.append("NoTag");
00074 
00075   Chi2Calc = new MadChi2Calc(0);
00076   Chi2Surf = NULL;
00077   ChiMin=99999.;
00078   NDOF = 0;
00079   Par1MinVal = 0.;
00080   Par2MinVal = 0.;
00081 
00082   EShiftErr = 0.;
00083 
00084   writeOut = true;
00085 
00086   if(!chainSR) {
00087     record = 0;
00088     strecord = 0;
00089     emrecord = 0;
00090     mcrecord = 0;
00091     threcord = 0;
00092     Clear();
00093     whichSource = -1;
00094     isMC = true;
00095     isTH = true;
00096     isEM = true;
00097     Nentries = -1;
00098     return;
00099   }
00100   
00101   //  InitChain(chainSR,chainMC,chainTH,chainEM);
00102   whichSource = 0;
00103 
00104 }

MadAnalysis::MadAnalysis ( JobC *  ,
string  ,
int   
)

Definition at line 107 of file MadAnalysis.cxx.

References BestFit, Chi2Calc, Chi2Surf, ChiMin, MadBase::Clear(), DataBkgdHist, DataHist, DataPOT, DataSignalHist, MadBase::emrecord, EShiftErr, MadBase::fChain, MadBase::fJC, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, MCBkgdHist, MCEvents, MCPOT, MadBase::mcrecord, MCSignalHist, NDOF, MadBase::Nentries, NumberOfBins, numDataEvts, numMCEvts, NumOscPars, OscillationFunction, OscillationParameters, Par1MinVal, Par2MinVal, RangeLowerLimit, RangeUpperLimit, MadBase::record, signalFracInData, MadBase::strecord, tag, MadBase::threcord, varyX, MadBase::whichSource, and writeOut.

00108    : MadQuantities(j,path,entries)
00109 {
00110 
00111   NumberOfBins = 1000;
00112   RangeLowerLimit = 0.;
00113   RangeUpperLimit = 50.;
00114 
00115   MCSignalHist[0] = NULL;   MCSignalHist[1] = NULL;
00116   MCBkgdHist[0] = NULL;     MCBkgdHist[1] = NULL;
00117   DataHist[0] = NULL;       DataHist[1] = NULL;
00118   DataSignalHist[0] = NULL; DataSignalHist[1] = NULL;
00119   DataPOT[0] = 0;           DataPOT[1] = 0;
00120   numDataEvts[0] = 0;       numDataEvts[1] = 0;
00121   signalFracInData[0] = 0;  signalFracInData[1] = 0;
00122   DataBkgdHist[0] = NULL;   DataBkgdHist[1] = NULL;
00123   MCEvents[0].clear();      MCEvents[1].clear();
00124   MCPOT[0] = 0;             MCPOT[1] = 0;
00125   numMCEvts[0] = 0;         numMCEvts[1] = 0;
00126   BestFit[0] = NULL;        BestFit[1] = NULL;
00127 
00128   OscillationFunction = NULL;
00129   OscillationParameters = NULL;
00130   NumOscPars = 0;
00131   varyX = new int[2];
00132   varyX[0] = 1;
00133   varyX[1] = 2;
00134 
00135   tag.append("NoTag");
00136 
00137   Chi2Calc = new MadChi2Calc(0);
00138   Chi2Surf = NULL;
00139   ChiMin=99999.;
00140   NDOF = 0;
00141   Par1MinVal = 0.;
00142   Par2MinVal = 0.;
00143 
00144   EShiftErr = 0.;
00145 
00146   writeOut = true;
00147 
00148   record = 0;
00149   strecord = 0;
00150   emrecord = 0;
00151   mcrecord = 0;
00152   threcord = 0;
00153   Clear();
00154   isMC = true;
00155   isTH = true;
00156   isEM = true;
00157   Nentries = entries;
00158   whichSource = 1;
00159   jcPath = path;
00160   fJC = j;
00161   fChain = NULL;
00162 
00163 }

MadAnalysis::~MadAnalysis (  )  [virtual]

Definition at line 166 of file MadAnalysis.cxx.

References BestFit, Chi2Calc, Chi2Surf, DataBkgdHist, DataHist, DataSignalHist, EndJob(), MCBkgdHist, MCSignalHist, OscillationFunction, varyX, and writeOut.

00167 {
00168 
00169   if(writeOut) {
00170     EndJob();
00171     cout << "Done EndJob()" << endl;
00172   }
00173   delete Chi2Calc;
00174   delete Chi2Surf;
00175   delete MCSignalHist[0];   delete MCSignalHist[1];
00176   delete MCBkgdHist[0];     delete MCBkgdHist[1];
00177   delete DataHist[0];       delete DataHist[1];
00178   delete DataSignalHist[0]; delete DataSignalHist[1];
00179   delete DataBkgdHist[0];   delete DataBkgdHist[1];
00180   delete BestFit[0];        delete BestFit[1];
00181   delete OscillationFunction;
00182   delete [] varyX;
00183   cout << "deleting MadAnalysis object" << endl;
00184 }

Member Function Documentation

virtual Bool_t MadAnalysis::AddUserBranches ( TTree *  tree  )  [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, and MadUserAnalysis.

Definition at line 117 of file MadAnalysis.h.

Referenced by CreatePAN().

00117 {if(tree) return true; return false;}

virtual void MadAnalysis::CallUserFunctions ( Int_t   )  [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, and MadUserAnalysis.

Definition at line 118 of file MadAnalysis.h.

00118 {return;}

void MadAnalysis::CreateANtpPAN ( std::string   ) 

Definition at line 976 of file MadAnalysis.cxx.

References ANtpTrackInfo::begPlane, ANtpEventInfo::begPlane, ANtpHeaderInfo::day, Munits::day, ANtpHeaderInfo::detector, ANtpTrackInfo::endPlane, ANtpEventInfo::endPlane, ANtpRecoInfo::eventLength, ANtpHeaderInfo::events, MadBase::eventSummary, GnumiInterface::FillANtpTruth(), ANtpInfoObjectFillerBeam::FillBeamMCTruthInformation(), ANtpInfoObjectFiller::FillEventInformation(), ANtpInfoObjectFiller::FillMCTruthInformation(), ANtpInfoObjectFiller::FillShowerInformation(), ANtpInfoObjectFiller::FillTrackInformation(), ANtpTrackInfo::fitMomentum, VldTimeStamp::GetDate(), VldContext::GetDetector(), RecDataHeader::GetRun(), VldTimeStamp::GetSec(), RecPhysicsHeader::GetSnarl(), RecDataHeader::GetSubRun(), VldTimeStamp::GetTime(), VldContext::GetTimeStamp(), RecHeader::GetVldContext(), ANtpRecoInfo::hadronicY, ANtpHeaderInfo::hour, Munits::hour, NtpSREvent::index, ANtpEventInfo::index, ANtpRecoInfo::inFiducialVolume, MadQuantities::IsFidAll(), MadQuantities::IsFidVtxEvt(), ANtpRecoInfo::isFullyContained, MadBase::isMC, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadTruthForRecoTH(), MadBase::mcrecord, ANtpHeaderInfo::minute, Munits::minute, ANtpHeaderInfo::month, month, ANtpRecoInfo::muDCosZVtx, ANtpRecoInfo::muEnergy, MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, MadBase::ntpTruth, ANtpRecoInfo::nuDCos, ANtpRecoInfo::nuEnergy, ANtpRecoInfo::pass, PassAnalysisCuts(), PassBasicCuts(), ANtpRecoInfo::passesCuts, PID(), ANtpRecoInfo::qeNuEnergy, ANtpRecoInfo::qeQ2, ANtpTrackInfo::rangeMomentum, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), ANtpShowerInfo::Reset(), ANtpTrackInfo::Reset(), ANtpEventInfo::Reset(), ANtpRecoInfo::Reset(), ANtpAnalysisInfo::Reset(), ANtpTruthInfoBeam::Reset(), ANtpHeaderInfo::Reset(), ANtpHeaderInfo::run, ANtpHeaderInfo::second, Munits::second, ANtpAnalysisInfo::separationParameter, ANtpRecoNtpManipulator::SetRecord(), ANtpRecoInfo::showerEnergy, ANtpTrackInfo::sigmaQoverP, ANtpRecoInfo::sigmaQoverP, ANtpHeaderInfo::snarl, GnumiInterface::Status(), NtpMCRecord::stdhep, NtpStRecord::stdhep, MadBase::strecord, ANtpHeaderInfo::subRun, ANtpRecoInfo::trackLength, ANtpRecoInfo::trackMomentum, ANtpRecoInfo::trackRange, ANtpHeaderInfo::utc, ANtpEventInfo::vtxX, ANtpRecoInfo::vtxX, ANtpEventInfo::vtxY, ANtpRecoInfo::vtxY, ANtpEventInfo::vtxZ, ANtpRecoInfo::vtxZ, ANtpHeaderInfo::year, and Munits::year.

00976                                             {
00977 
00978   std::string savename = "PAN_" + tag + ".root";
00979   TFile save(savename.c_str(),"RECREATE"); 
00980   save.cd();
00981   
00982   GnumiInterface gnumi;
00983   if(!gnumi.Status()) {
00984     cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00985          << " Will not be filling NuParent info"
00986          << endl;
00987     cout << "Set environmental variable $GNUMIAUX to point to the "
00988          << "directory containing the gnumi files to fix this!"
00989          << endl;
00990 
00991   }
00992 
00993   //PAN Quantities: See AnalysisNtuples package for info
00994   ANtpHeaderInfo *antpHeader       = new ANtpHeaderInfo();
00995   ANtpTruthInfoBeam *antpTruth     = new ANtpTruthInfoBeam();
00996   ANtpAnalysisInfo *antpAnalysis   = new ANtpAnalysisInfo();
00997   ANtpRecoInfo *antpReco           = new ANtpRecoInfo();
00998   ANtpEventInfo *antpEvent         = new ANtpEventInfo();
00999   ANtpTrackInfo *antpTrack         = new ANtpTrackInfo();
01000   ANtpShowerInfo *antpShower       = new ANtpShowerInfo();
01001   ANtpRecoNtpManipulator *ntpManip = new ANtpRecoNtpManipulator();
01002   ANtpInfoObjectFillerBeam filla;
01003 
01004   //Quantities missing from ANtp classes:
01005   Int_t mcevent;          // mc event index
01006   Int_t is_mc;            // is a corresponding mc event found
01007   
01008   //Tree:
01009   TTree *tree = new TTree("pan","pan");
01010   tree->Branch("header","ANtpHeaderInfo",&antpHeader,8000,1);
01011   if(isMC) {
01012     tree->Branch("is_mc",&is_mc,"is_mc/I",32000);    
01013     tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
01014     tree->Branch("truth","ANtpTruthInfoBeam",&antpTruth,8000,1);
01015   }
01016   tree->Branch("reco","ANtpRecoInfo",&antpReco,8000,1);
01017   tree->Branch("analysis","ANtpAnalysisInfo",&antpAnalysis,8000,1);
01018   tree->Branch("event","ANtpEventInfo",&antpEvent,8000,1);
01019   tree->Branch("track","ANtpTrackInfo",&antpTrack,8000,1);
01020   tree->Branch("shower","ANtpShowerInfo",&antpShower,8000,1);
01021 
01022   for(int i=0;i<Nentries;i++){
01023     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
01024                             << "% done" << std::endl;
01025     
01026     if(!this->GetEntry(i)) continue;
01027     if(eventSummary->nevent==0) continue;
01028     
01029     //fill tree once for each reconstructed event
01030     for(int j=0;j<eventSummary->nevent;j++){ 
01031       
01032       if(!LoadEvent(j)) continue; //no event found
01033 
01034       //zero all tree values
01035       antpHeader->Reset();
01036       is_mc = 0;
01037       mcevent = -1;
01038       antpTruth->Reset();
01039       antpAnalysis->Reset();
01040       antpReco->Reset();
01041       antpEvent->Reset();
01042       antpTrack->Reset();
01043       antpShower->Reset();
01044       ntpManip->SetRecord(strecord);
01045 
01046       antpHeader->events = eventSummary->nevent;    
01047       antpHeader->run    = ntpHeader->GetRun();
01048       antpHeader->subRun = ntpHeader->GetSubRun();
01049       antpHeader->snarl  = ntpHeader->GetSnarl();
01050 
01051       UInt_t year,month,day,hour,minute,second;
01052       ntpHeader->GetVldContext().GetTimeStamp().GetDate(kFALSE,0,&year,
01053                                                         &month,&day);
01054       ntpHeader->GetVldContext().GetTimeStamp().GetTime(kFALSE,0,&hour,
01055                                                         &minute,&second);
01056       antpHeader->year   = year;
01057       antpHeader->month  = month;
01058       antpHeader->day    = day;
01059       antpHeader->hour   = hour;
01060       antpHeader->minute = minute;
01061       antpHeader->second = second;
01062       antpHeader->utc    = ntpHeader->GetVldContext().GetTimeStamp().GetSec();
01063 
01064       antpEvent->index = j;
01065       filla.FillEventInformation(ntpManip,ntpEvent,antpEvent);
01066       
01067       //get detector type:
01068       antpHeader->detector = ntpHeader->GetVldContext().GetDetector();
01069       
01070       //fiducial criteria on vtx
01071       antpReco->inFiducialVolume = 1;
01072       if(!IsFidVtxEvt(ntpEvent->index)) 
01073         antpReco->inFiducialVolume = 0;    
01074       
01075       //check for a corresponding mc event      
01076       if(LoadTruthForRecoTH(antpEvent->index,mcevent)) {
01077         //true info for tree:
01078         is_mc = 1;
01079         filla.FillMCTruthInformation(ntpTruth,antpTruth);
01080         if(isST) filla.FillBeamMCTruthInformation(ntpTruth,strecord->stdhep,antpTruth);
01081         else filla.FillBeamMCTruthInformation(ntpTruth,mcrecord->stdhep,antpTruth);
01082         if(gnumi.Status()) {
01083           if(isST) gnumi.FillANtpTruth(strecord,mcevent,*antpTruth);
01084           else gnumi.FillANtpTruth(mcrecord,mcevent,*antpTruth);
01085         }
01086       }
01087 
01088       if(PassBasicCuts(antpEvent->index)) {
01089         
01090         int track_index                 = -1;
01091         LoadLargestTrackFromEvent(antpEvent->index,track_index);
01092         if(track_index!=-1)  filla.FillTrackInformation(ntpTrack,antpTrack);
01093         
01094         int shower_index                = -1;
01095         LoadLargestShowerFromEvent(antpEvent->index,shower_index);
01096         if(shower_index!=-1) filla.FillShowerInformation(ntpShower,antpShower);
01097 
01098         antpReco->passesCuts        = 1;
01099         antpReco->eventLength   = TMath::Abs(antpEvent->endPlane - 
01100                                                      antpEvent->begPlane);
01101         antpReco->trackLength   = TMath::Abs(antpTrack->endPlane -
01102                                                      antpTrack->begPlane);
01103         antpReco->trackMomentum = antpTrack->fitMomentum;
01104         antpReco->trackRange    = antpTrack->rangeMomentum;
01105         antpReco->sigmaQoverP   = antpTrack->sigmaQoverP;
01106         
01107         antpReco->muEnergy      = RecoMuEnergy(0,track_index);
01108         antpReco->showerEnergy  = RecoShwEnergy(shower_index);
01109         antpReco->nuEnergy      = ( antpReco->muEnergy + 
01110                                             antpReco->showerEnergy );
01111         
01112         antpReco->qeNuEnergy     = RecoQENuEnergy(track_index);
01113         antpReco->qeQ2           = RecoQEQ2(track_index);
01114         if (antpReco->nuEnergy>0) {
01115           antpReco->hadronicY    = ( antpReco->showerEnergy / 
01116                                              antpReco->nuEnergy );
01117         }
01118         
01119         antpReco->muDCosZVtx    = RecoMuDCosZ(track_index);
01120         antpReco->nuDCos        = RecoMuDCosNeu(track_index);
01121         antpReco->vtxX          = antpEvent->vtxX;
01122         antpReco->vtxY          = antpEvent->vtxY;
01123         antpReco->vtxZ          = antpEvent->vtxZ;
01124         antpReco->isFullyContained  = IsFidAll(track_index);
01125 
01126         antpAnalysis->separationParameter = PID(antpEvent->index,0);
01127         if(PassAnalysisCuts(antpEvent->index)) antpReco->pass = 1;
01128         else antpReco->pass = 0;
01129       }
01130       //fill the tree
01131       tree->Fill();
01132     }
01133   }
01134 
01135   delete antpHeader;
01136   delete antpTruth;
01137   delete antpAnalysis;
01138   delete antpReco;
01139   delete antpEvent;
01140   delete antpTrack;
01141   delete antpShower;
01142 
01143   save.cd();
01144   tree->Write();
01145   save.Write();
01146   save.Close();
01147 }

void MadAnalysis::CreatePAN ( std::string   ) 

Definition at line 244 of file MadAnalysis.cxx.

References AddUserBranches(), MadBase::dmxStatus, MadBase::eventSummary, BeamMonSpill::fBpmInt, MadQuantities::fGeantinoEnergy, MadQuantities::fMaxFSKaonEnergy, MadQuantities::fMaxFSNeutronEnergy, MadQuantities::fMaxFSPiMinusEnergy, MadQuantities::fMaxFSPiPlusEnergy, MadQuantities::fMaxFSPiZeroEnergy, MadQuantities::fMaxFSProtonEnergy, MadQuantities::fNumFSGeantino, MadQuantities::fNumFSKaon, MadQuantities::fNumFSNeutron, MadQuantities::fNumFSPiMinus, MadQuantities::fNumFSPiPlus, MadQuantities::fNumFSPiZero, MadQuantities::fNumFSProton, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTargProfX, BeamMonSpill::fTargProfY, MadQuantities::fTotFSKaonEnergy, MadQuantities::fTotFSNeutronEnergy, MadQuantities::fTotFSPiMinusEnergy, MadQuantities::fTotFSPiPlusEnergy, MadQuantities::fTotFSPiZeroEnergy, MadQuantities::fTotFSProtonEnergy, BDSpillAccessor::Get(), DataUtil::GetDetector(), MadBase::GetEntry(), RecPhysicsHeader::GetRemoteSpillType(), RecDataHeader::GetRun(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), NtpSREvent::index, NtpSRSlice::index, SpillTimeFinder::Instance(), MadQuantities::IsFidVtxEvt(), Detector::kFar, SimFlag::kMC, Detector::kNear, NtpSREventSummary::litime, MadBase::LoadEvent(), MadBase::LoadSlice(), BDSpillAccessor::LoadSpill(), MadBase::LoadStrip(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSRDmxStatus::nonphysicalfail, NtpSREvent::nshower, NtpSREvent::nstrip, NtpSREventSummary::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpSlice, MadBase::ntpStrip, NtpSREvent::ntrack, NtpSRSlice::plane, NtpSRStrip::plane, run(), BMSpillAna::SelectSpill(), BMSpillAna::SetSpill(), MadQuantities::SetStdHepVars(), BMSpillAna::SetTimeDiff(), NtpSREvent::slc, GnumiInterface::Status(), NtpSREvent::stp, NtpSRStrip::strip, NtpSRVertex::t, NtpSRStrip::time0, NtpSRStrip::time1, NtpSREventSummary::trigtime, NtpSRDmxStatus::validplanesfail, NtpSRDmxStatus::vertexplanefail, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, NtpSRVertex::z, and NuParent::Zero().

00244                                         {
00245 
00246   this->GetEntry(0);
00247   const bool file_is_mc
00248     =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
00249 
00250   std::string savename = "PAN_" + tag + ".root";
00251   TFile save(savename.c_str(),"RECREATE");
00252   save.cd();
00253   
00254   GnumiInterface *gnumi = 0;
00255   if(file_is_mc){
00256     gnumi = new GnumiInterface();
00257     if(!gnumi->Status()) {
00258       cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00259            << " Will not be filling NuParent info"
00260            << endl;
00261       cout << "Set environmental variable $GNUMIAUX to point to the "
00262            << "directory containing the gnumi files to fix this!"
00263            << endl;
00264     }
00265     else {
00266       const char* gnumidir= getenv("GNUMIAUX");
00267       cout<<"Read gnumi files from "<<gnumidir<<endl;
00268     }
00269   }
00270 
00271   //make a BMSpillAna so we can tell if it's goodbeam
00272   BMSpillAna bmana;
00273 
00274   //PAN Quantities 
00275   //Truth:
00276   Float_t true_enu;       // true neutrino energy (GeV)
00277   Float_t true_pxnu;      // true neutrino momentum-x (GeV)
00278   Float_t true_pynu;      // true neutrino momentum-y (GeV)
00279   Float_t true_pznu;      // true neutrino momentum-z (GeV)
00280   Float_t true_etgt;       // true target energy (GeV)
00281   Float_t true_pxtgt;      // true target momentum-x (GeV)
00282   Float_t true_pytgt;      // true target momentum-y (GeV)
00283   Float_t true_pztgt;      // true target momentum-z (GeV)
00284   Float_t true_emu;       // true muon energy
00285   Float_t true_elep;       // true muon energy
00286   Float_t true_eshw;      // true shower energy
00287   Float_t true_x;         // true x
00288   Float_t true_y;         // true y
00289   Float_t true_q2;        // true q2
00290   Float_t true_w2;        // true w2
00291   Float_t true_emfrac;    // true emfrac
00292   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
00293   Float_t true_dircosz;   // track z-direction cosine
00294   Float_t true_vtxx;      // true x vtx
00295   Float_t true_vtxy;      // true y vtx
00296   Float_t true_vtxz;      // true z vtx
00297 
00298   //TruthHelper:
00299   Float_t th_evt_pur;     // truth helper event purity
00300   Float_t th_evt_all;     // truth helper event completeness all
00301   Float_t th_evt_slc;     // truth helper event completeness slc
00302   Float_t th_shw_pur;     // truth helper shower purity
00303   Float_t th_shw_all;     // truth helper shower completeness all
00304   Float_t th_shw_slc;     // truth helper shower completeness slc
00305   Float_t th_trk_pur;     // truth helper track purity
00306   Float_t th_trk_all;     // truth helper track completeness all
00307   Float_t th_trk_slc;     // truth helper track completeness slc
00308 
00309   //For Neugen:
00310   Int_t flavour;          // true flavour: 1-e 2-mu 3-tau
00311   Int_t nooscflavour;     // true flavour before osc: 1-e 2-mu 3-tau
00312   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
00313   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
00314   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
00315   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
00316   Int_t had_fs;           // hadronic final state, number between 200-300
00317   
00318   //For Beam Reweighting:
00319   NuParent *nuparent = new NuParent();
00320 
00321   //Reco Quantities
00322   Int_t detector;         // Near = 1, Far = 2;
00323   Int_t run;              // run number
00324   Int_t snarl;            // snarl number
00325   Int_t subrun;           // subrun number
00326   UInt_t trgsrc;          // trigger source
00327   double trgtime;         // trigger time
00328   Int_t spilltype;        // comes from mdSpillTypeEnum
00329                           // 0=none, 1=reported, 2=predicted
00330                           // 3=fake, 4=local
00331   Int_t nevent;           // number of events in snarl
00332   Int_t event;            // event index
00333   Int_t slice;            // slice
00334   Int_t mcevent;          // mc event index
00335   Int_t ntrack;           // number of tracks in event
00336   Int_t nshower;          // number of showers in event
00337   double litime;          //time last li event in snarl, -1 if none
00338   UChar_t dmx_stat;       //demux failures (for fd)
00339 
00340   Int_t is_fid;           // pass fid cut. true = 1, false = 0
00341   Int_t is_cev;           // fully contained. true = 1, false = 0 
00342   Int_t is_mc;            // is a corresponding mc event found
00343   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
00344   Float_t pid0;           // pid parameter (usual method)
00345   Float_t pid1;           // pid parameter (alternative method 1)
00346   Float_t pid2;           // pid parameter (alternative method 2)
00347   Int_t pass;             // pass analysis cuts. true = 1, false = 0
00348 
00349   Float_t reco_enu;       // reco neutrino energy
00350   Float_t reco_emu;       // reco muon energy
00351   Float_t reco_eshw;      // reco shower energy (generic)
00352   Float_t reco_eshw_linCC;// reco shower energy using lin CC method
00353   Float_t reco_eshw_wtCC; // reco shower energy using deweighting for CC 
00354   Float_t reco_eshw_linNC;// reco shower energy using lin NC method
00355   Float_t reco_eshw_wtNC; // reco shower energy using deweighting for NC 
00356   Float_t reco_eshw_em;   // reco shower energy for em showers
00357   Float_t reco_qe_enu;    // reco qe neutrino energy
00358   Float_t reco_qe_q2;     // reco qe q2
00359   Float_t reco_y;         // reco y
00360   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
00361   Float_t reco_dircosz;   // reco track vtx z-dircos
00362   Float_t reco_vtxx;      // reco x vtx
00363   Float_t reco_vtxy;      // reco y vtx
00364   Float_t reco_vtxz;      // reco z vtx
00365   Float_t reco_emfrac;    // reco emfrac
00366   Float_t reco_emfrac_prob;// reco emfrac with fit prob constraint
00367 
00368   Float_t evtlength;      // reco event length
00369   Float_t trklength;      // reco track length
00370   Float_t shwlength;      // reco shower length
00371   Float_t slclength;      // reco slice length
00372   Double_t closestevent_s; // spatial distance to closest event in snarl
00373   Double_t closestevent_t; // time to closest event in snarl
00374   Float_t fracRehitStrip;  // fraction of hits in event hit 
00375 
00376   Float_t trkmomrange;    // reco track momentum from range
00377   Float_t trkqp;          // reco track q/p
00378   Float_t trkeqp;         // reco track q/p error
00379   Float_t trkphfrac;      // reco pulse height fraction in track
00380   Float_t trkphplane;     // reco average track pulse height/plane
00381 
00382   Double_t shwstptimemin; // min strip time in shower
00383   Double_t shwstptimemax; // max strip time in shower
00384   Double_t shwstptimerms; // phw rms of strip times in shower
00385   
00386   // beam monitoring variables
00387   Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill;
00388   //goodbeam==1 if spill meets criteria defined below 
00389   //(see the line where goodbeam gets set)
00390   //beamcnf is an int as defined in BeamMonSpill::StatusBits
00391   UInt_t beamcnf;
00392   //leaving nuTarZ in until we've phased out the summary ntuples
00393   Double_t nuTarZ;
00394   //goodbeam==0 if spill doesnt meet these criteria
00395   Int_t goodbeam;
00396 
00397   //these next variables weren't available in the old beam monitoring ntuples
00398   UInt_t horn_on, target_in, n_batches;
00399   Float_t tor101, tr101d, tortgt, trtgtd;
00400   Float_t hadmon, mumon1, mumon2, mumon3;
00401 
00402   //Trees
00403   TTree *tree = new TTree("pan","pan");
00404   //Truth
00405   tree->Branch("true_enu",&true_enu,"true_enu/F",32000);
00406   tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",32000);
00407   tree->Branch("true_pynu",&true_pynu,"true_pynu/F",32000);
00408   tree->Branch("true_pznu",&true_pznu,"true_pznu/F",32000);
00409   tree->Branch("true_etgt",&true_etgt,"true_etgt/F",32000);
00410   tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",32000);
00411   tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",32000);
00412   tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",32000);
00413   tree->Branch("true_emu",&true_emu,"true_emu/F",32000);
00414   tree->Branch("true_elep",&true_elep,"true_elep/F",32000);
00415   tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000);
00416   tree->Branch("true_x",&true_x,"true_x/F",32000);
00417   tree->Branch("true_y",&true_y,"true_y/F",32000);
00418   tree->Branch("true_q2",&true_q2,"true_q2/F",32000);
00419   tree->Branch("true_w2",&true_w2,"true_w2/F",32000);
00420   tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000);
00421   tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000);
00422   tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000);
00423   tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000);
00424   tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000);
00425   tree->Branch("true_emfrac",&true_emfrac,"true_emfrac/F",32000);
00426   //Truth Helper:
00427   tree->Branch("th_evt_pur",&th_evt_pur,"th_evt_pur/F",32000);
00428   tree->Branch("th_evt_all",&th_evt_all,"th_evt_all/F",32000);
00429   tree->Branch("th_evt_slc",&th_evt_slc,"th_evt_slc/F",32000);
00430   tree->Branch("th_shw_pur",&th_shw_pur,"th_shw_pur/F",32000);
00431   tree->Branch("th_shw_all",&th_shw_all,"th_shw_all/F",32000);
00432   tree->Branch("th_shw_slc",&th_shw_slc,"th_shw_slc/F",32000);
00433   tree->Branch("th_trk_pur",&th_trk_pur,"th_trk_pur/F",32000);
00434   tree->Branch("th_trk_all",&th_trk_all,"th_trk_all/F",32000);
00435   tree->Branch("th_trk_slc",&th_trk_slc,"th_trk_slc/F",32000);
00436   //stdhep
00437   tree->Branch("stdhep_ng",&fNumFSGeantino,"stdhep_ng/I",32000);
00438   tree->Branch("stdhep_eg",&fGeantinoEnergy,"stdhep_eg/F",32000);
00439   tree->Branch("stdhep_np",&fNumFSProton,"stdhep_np/I",32000);
00440   tree->Branch("stdhep_ep",&fTotFSProtonEnergy,"stdhep_ep/F",32000);
00441   tree->Branch("stdhep_mp",&fMaxFSProtonEnergy,"stdhep_mp/F",32000);
00442   tree->Branch("stdhep_nn",&fNumFSNeutron,"stdhep_nn/I",32000);
00443   tree->Branch("stdhep_en",&fTotFSNeutronEnergy,"stdhep_en/F",32000);
00444   tree->Branch("stdhep_mn",&fMaxFSNeutronEnergy,"stdhep_mn/F",32000);
00445   tree->Branch("stdhep_npp",&fNumFSPiPlus,"stdhep_npp/I",32000);
00446   tree->Branch("stdhep_epp",&fTotFSPiPlusEnergy,"stdhep_epp/F",32000);
00447   tree->Branch("stdhep_mpp",&fMaxFSPiPlusEnergy,"stdhep_mpp/F",32000);
00448   tree->Branch("stdhep_npm",&fNumFSPiMinus,"stdhep_npm/I",32000);
00449   tree->Branch("stdhep_epm",&fTotFSPiMinusEnergy,"stdhep_epm/F",32000);
00450   tree->Branch("stdhep_mpm",&fMaxFSPiMinusEnergy,"stdhep_mpm/F",32000);
00451   tree->Branch("stdhep_npz",&fNumFSPiZero,"stdhep_npz/I",32000);
00452   tree->Branch("stdhep_epz",&fTotFSPiZeroEnergy,"stdhep_epz/F",32000);
00453   tree->Branch("stdhep_mpz",&fMaxFSPiZeroEnergy,"stdhep_mpz/F",32000);
00454   tree->Branch("stdhep_nk",&fNumFSKaon,"stdhep_nk/I",32000);
00455   tree->Branch("stdhep_ek",&fTotFSKaonEnergy,"stdhep_ek/F",32000);
00456   tree->Branch("stdhep_mk",&fMaxFSKaonEnergy,"stdhep_mk/F",32000);
00457   //Neugen
00458   tree->Branch("flavour",&flavour,"flavour/I",32000);
00459   tree->Branch("nooscflavour",&nooscflavour,"nooscflavour/I",32000);
00460   tree->Branch("process",&process,"process/I",32000);
00461   tree->Branch("nucleus",&nucleus,"nucleus/I",32000);
00462   tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000);
00463   tree->Branch("initial_state",&initial_state,"initial_state/I",32000);
00464   tree->Branch("had_fs",&had_fs,"had_fs/I",32000);
00465   //NuParent
00466   tree->Branch("nuparent","NuParent",&nuparent,8000,1);
00467   //Reco
00468   tree->Branch("detector",&detector,"detector/I",32000);
00469   tree->Branch("run",&run,"run/I",32000);
00470   tree->Branch("snarl",&snarl,"snarl/I",32000);
00471   tree->Branch("subrun",&subrun,"subrun/I",32000);
00472   tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000);
00473   tree->Branch("trgtime",&trgtime,"trgtime/D",32000);
00474   tree->Branch("spilltype",&spilltype,"spiltype/I",32000);
00475   tree->Branch("event",&event,"event/I",32000);
00476   tree->Branch("slice",&slice,"slice/I",32000);
00477   tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
00478   tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
00479   tree->Branch("nshower",&nshower,"nshower/I",32000);
00480   tree->Branch("litime",&litime,"litime/D");
00481   tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b");
00482 
00483   tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
00484   tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
00485   tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
00486   tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
00487   tree->Branch("pid0",&pid0,"pid0/F",32000);
00488   tree->Branch("pid1",&pid1,"pid1/F",32000);
00489   tree->Branch("pid2",&pid2,"pid2/F",32000);
00490   tree->Branch("pass",&pass,"pass/I",32000);
00491 
00492   tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
00493   tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
00494   tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
00495   tree->Branch("reco_eshw_linCC",&reco_eshw_linCC,"reco_eshw_linCC/F",32000);
00496   tree->Branch("reco_eshw_wtCC",&reco_eshw_wtCC,"reco_eshw_wtCC/F",32000);
00497   tree->Branch("reco_eshw_linNC",&reco_eshw_linNC,"reco_eshw_linNC/F",32000);
00498   tree->Branch("reco_eshw_wtNC",&reco_eshw_wtNC,"reco_eshw_wtNC/F",32000);
00499   tree->Branch("reco_eshw_em",&reco_eshw_em,"reco_eshw_em/F",32000);
00500   tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
00501   tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);
00502   tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
00503   tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
00504   tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);
00505   tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
00506   tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
00507   tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
00508   tree->Branch("reco_emfrac",&reco_emfrac,"reco_emfrac/F",32000);
00509   tree->Branch("reco_emfrac_prob",&reco_emfrac_prob,
00510                "reco_emfrac_prob/F",32000);
00511   
00512   tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
00513   tree->Branch("trklength",&trklength,"trklength/F",32000);
00514   tree->Branch("shwlength",&shwlength,"shwlength/F",32000);
00515   tree->Branch("slclength",&slclength,"slclength/F",32000);
00516   tree->Branch("closestevent_s",&closestevent_s,"closestevent_s/D",32000);
00517   tree->Branch("closestevent_t",&closestevent_t,"closestevent_t/D",32000);
00518   tree->Branch("fracRehitStrip",&fracRehitStrip,"fracRehitStrip/F",32000);
00519   
00520   tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
00521   tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
00522   tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
00523   tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
00524   tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
00525 
00526   tree->Branch("shwstptimemin",&shwstptimemin,"shwstptimemin/D",32000);
00527   tree->Branch("shwstptimemax",&shwstptimemax,"shwstptimemax/D",32000);
00528   tree->Branch("shwstptimerms",&shwstptimerms,"shwstptimerms/D",32000);
00529 
00530   tree->Branch("hbw",&hbw,"hbw/D");
00531   tree->Branch("vbw",&vbw,"vbw/D");
00532   tree->Branch("hpos1",&hpos1,"hpos1/D");
00533   tree->Branch("vpos1",&vpos1,"vpos1/D");
00534   tree->Branch("hpos2",&hpos2,"hpos2/D");
00535   tree->Branch("vpos2",&vpos2,"vpos2/D");
00536   tree->Branch("hornI",&hornI,"hornI/D");
00537   tree->Branch("closestspill",&closestspill,"closestspill/D");
00538   tree->Branch("beamcnf",&beamcnf,"beamcnf/i");
00539   tree->Branch("nuTarZ",&nuTarZ,"nuTarZ/D");
00540   tree->Branch("goodbeam",&goodbeam,"goodbeam/I");
00541   tree->Branch("horn_on",&horn_on,"horn_on/i");
00542   tree->Branch("target_in",&target_in,"target_in/i");
00543   tree->Branch("n_batches",&n_batches,"n_batches/i");
00544   tree->Branch("tor101",&tor101,"tor101/F");
00545   tree->Branch("tr101d",&tr101d,"tr101d/F");
00546   tree->Branch("tortgt",&tortgt,"tortgt/F");
00547   tree->Branch("trtgtd",&trtgtd,"trtgtd/F");
00548   tree->Branch("hadmon",&hadmon,"hadmon/F");
00549   tree->Branch("mumon1",&mumon1,"mumon1/F");
00550   tree->Branch("mumon2",&mumon2,"mumon2/F");
00551   tree->Branch("mumon3",&mumon3,"mumon3/F");
00552   
00553   this->AddUserBranches(tree);
00554 
00555   //pot counting
00556   float pottor101=0.;
00557   float pottortgt=0.;
00558   float pot=0.;
00559   int NRUNS=0;
00560   int NSNARLS=0;
00561   TTree *pottree = new TTree("pottree","pottree");
00562   pottree->Branch("pot",&pot,"pot/F");
00563   pottree->Branch("pottor101",&pottor101,"pottor101/F");
00564   pottree->Branch("pottortgt",&pottortgt,"pottortgt/F");
00565   pottree->Branch("nruns",&NRUNS,"nruns/I");
00566   pottree->Branch("nsnarls",&NSNARLS,"nsnarls/I");
00567 
00568   int lastrun=-1;
00569   for(int i=0;i<Nentries;i++){
00570     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00571                             << "% done" << std::endl;
00572 
00573     if(!this->GetEntry(i)) continue;
00574     NSNARLS++;
00575     nevent = eventSummary->nevent;
00576     run = ntpHeader->GetRun();
00577     if(run!=lastrun){
00578       NRUNS++;
00579       lastrun = run;
00580     }
00581     snarl = ntpHeader->GetSnarl();
00582     trgsrc = ntpHeader->GetTrigSrc();
00583     spilltype = ntpHeader->GetRemoteSpillType();
00584     trgtime = eventSummary->trigtime;
00585     litime = eventSummary->litime;    
00586     dmx_stat=0;
00587     dmx_stat+=(dmxStatus->nonphysicalfail);
00588     dmx_stat+=((dmxStatus->validplanesfail<<1));
00589     dmx_stat+=((dmxStatus->vertexplanefail<<2));
00590 
00592 
00593     hbw = vbw = hpos1 = vpos1 = hpos2 = vpos2 = 0.0;
00594     nuTarZ = hornI = closestspill = 0.0;
00595     beamcnf = horn_on = target_in = n_batches = 0;
00596     goodbeam=0;
00597     tor101 = tr101d = tortgt = trtgtd = 0.0;
00598     hadmon = mumon1 = mumon2 = mumon3 = 0.0;
00599     
00600     if(!file_is_mc){ //file is mc
00601       VldContext vc = ntpHeader->GetVldContext();
00602       VldTimeStamp vts = vc.GetTimeStamp();
00603       BDSpillAccessor &sa = BDSpillAccessor::Get();
00604       const BeamMonSpill *bs = sa.LoadSpill(vts);
00605       hbw=bs->fProfWidX*1000.;//make it in mm
00606       vbw=bs->fProfWidY*1000.;//make it in mm
00607       hpos1=bs->fTargProfX*1000.;//make it in mm
00608       vpos1=bs->fTargProfY*1000.;//make it in mm
00609       n_batches=0;
00610       float btint=0.;
00611       hpos2=0.;
00612       vpos2=0.;
00613       for(int bt=0;bt<6;bt++){
00614         hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
00615         vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);
00616         if(bs->fBpmInt[bt]>0.001){
00617           n_batches++;
00618         }
00619         btint+=bs->fBpmInt[bt];
00620       }
00621       if(btint!=0){
00622         hpos2=hpos2*1000./btint;//make it in mm
00623         vpos2=vpos2*1000./btint;//make it in mm
00624       }
00625       else{
00626         hpos2=-999.;
00627         vpos2=-999.;
00628       }
00629       
00630       hornI=bs->fHornCur;       
00631       SpillTimeFinder &sfind= SpillTimeFinder::Instance();
00632       closestspill = sfind.GetTimeToNearestSpill(vc);
00633       
00634       beamcnf =bs->GetStatusBits().beam_type;
00635       nuTarZ=0.;
00636       horn_on = bs->GetStatusBits().horn_on;
00637       target_in = bs->GetStatusBits().target_in;
00638       
00639       tor101=bs->fTor101;
00640       tr101d=bs->fTr101d;
00641       tortgt=bs->fTortgt;
00642       trtgtd=bs->fTrtgtd;
00643       hadmon=bs->fHadInt;
00644       mumon1=bs->fMuInt1;
00645       mumon2=bs->fMuInt2;
00646       mumon3=bs->fMuInt3;
00647       
00648       goodbeam=0;
00649       
00650       //use Mark D. BMSpillAna to set good beam
00651       bmana.SetSpill(*bs);
00652       bmana.SetTimeDiff(closestspill);
00653       if(bmana.SelectSpill()){
00654         goodbeam=1;
00655       }
00656       else{
00657         goodbeam=0;
00658       }
00659       
00660       if(goodbeam==1&&spilltype!=3){//don't count fake spills
00661         pot+=tortgt;
00662         pottor101+=tor101;
00663         pottortgt+=tortgt;
00664       }
00665     }    
00666     
00667     std::vector<Double_t> evtx(nevent,0);
00668     std::vector<Double_t> evty(nevent,0);
00669     std::vector<Double_t> evtz(nevent,0);
00670     std::vector<Double_t> evtt(nevent,0);
00671     //get vertex for each event in slice
00672     for(int j=0;j<nevent;j++){ 
00673       if(!LoadEvent(j)) {       
00674         evtx[j] = -100.;
00675         evty[j] = -100.;
00676         evtz[j] = -100.;
00677         evtt[j] = -100.;
00678       }
00679       evtx[j] = ntpEvent->vtx.x;
00680       evty[j] = ntpEvent->vtx.y;
00681       evtz[j] = ntpEvent->vtx.z;
00682       evtt[j] = ntpEvent->vtx.t;
00683     }
00684 
00685     //find and store earliest time of hit strip
00686     Double_t stripArrayTime[500][192][2];
00687     for(int ii=0;ii<500;ii++) { 
00688       for(int jj=0;jj<192;jj++) {
00689         stripArrayTime[ii][jj][0] = stripArrayTime[ii][jj][1] = -1.0;
00690       }
00691     }
00692     for(int j=0;j<int(eventSummary->nstrip);j++){
00693       if(!LoadStrip(j)) continue;
00694       Int_t pl = ntpStrip->plane;
00695       Int_t st = ntpStrip->strip;
00696       if(ntpStrip->time1<stripArrayTime[pl][st][1] || 
00697          stripArrayTime[pl][st][1]<=0) {
00698         stripArrayTime[pl][st][1] = ntpStrip->time1;
00699       }
00700       if(ntpStrip->time0<stripArrayTime[pl][st][0] || 
00701          stripArrayTime[pl][st][0]<=0) {
00702         stripArrayTime[pl][st][0] = ntpStrip->time0;
00703       }
00704     }
00705 
00706     if(nevent==0) continue;
00707     //fill tree once for each reconstructed event
00708     for(int j=0;j<nevent;j++){ 
00709       if(!LoadEvent(j)) continue; //no event found 
00710       //set event index:
00711       event = j;
00712       if(LoadSlice(ntpEvent->slc)) {
00713         slice = ntpSlice->index;
00714         slclength = ntpSlice->plane.n;
00715       }
00716       ntrack = ntpEvent->ntrack;
00717       nshower = ntpEvent->nshower;
00718       
00719       //zero all tree values
00720       true_enu = 0; true_emu = 0; true_elep = 0; true_eshw = 0; 
00721       true_pxnu = 0; true_pynu = 0; true_pznu = 0;
00722       true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
00723       true_dircosneu = 0; true_dircosz = 0; true_emfrac = 0;
00724       true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
00725       th_evt_pur = 0; th_evt_all = 0; th_evt_slc = 0;
00726       th_shw_pur = 0; th_shw_all = 0; th_shw_slc = 0;
00727       th_trk_pur = 0; th_trk_all = 0; th_trk_slc = 0;
00728       flavour = 0; nooscflavour = 0; process = 0; nucleus = 0; cc_nc = 0;
00729       initial_state = 0; had_fs = 0;
00730       true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;
00731       nuparent->Zero();
00732       detector = -1; mcevent = -1;
00733       is_fid = is_cev = is_mc = pass_basic = pass = 0;
00734       pid0 = pid1 = pid2 = 0.; 
00735       reco_enu = reco_emu = reco_eshw = 0; 
00736       reco_eshw_linCC = reco_eshw_wtCC = reco_eshw_linNC = 0; 
00737       reco_eshw_wtNC = reco_eshw_em = 0;
00738       reco_qe_enu = reco_qe_q2 = reco_y = reco_dircosneu = 0; 
00739       reco_dircosz = reco_emfrac = reco_emfrac_prob = 0;
00740       reco_vtxx = reco_vtxy = reco_vtxz = 0;
00741       evtlength = trklength = shwlength = 0;
00742       closestevent_s = closestevent_t = -1.;
00743       fracRehitStrip = 0;
00744       trkphfrac = trkphplane = trkmomrange = trkqp = trkeqp = 0;
00745       shwstptimemin = shwstptimemax = shwstptimerms = 0.0;
00746       this->SetStdHepVars(-1); //zero values
00747 
00748       //get detector type:
00749       if(ntpHeader->GetVldContext().
00750          GetDetector()==Detector::kFar) detector = 2;
00751       else if(ntpHeader->GetVldContext().
00752               GetDetector()==Detector::kNear) detector = 1;
00753 
00754       //fiducial vertex
00755       is_fid = 1;
00756       if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00757 
00758       //data cleaning      
00759       for(int k=0; k<ntpEvent->nstrip; k++){
00760         if(!LoadStrip(ntpEvent->stp[k])) continue;
00761         Int_t pl = ntpStrip->plane;
00762         Int_t st = ntpStrip->strip;
00763         if(ntpStrip->time0>stripArrayTime[pl][st][0]) {
00764           fracRehitStrip+=1;
00765         }
00766         else if(ntpStrip->time1>stripArrayTime[pl][st][1]) {
00767           fracRehitStrip+=1;
00768         }
00769       }
00770       fracRehitStrip/=ntpEvent->nstrip;
00771 
00772       Bool_t good_truth = false;
00773       //check for a corresponding mc event      
00774       if(LoadTruthForRecoTH(j,mcevent)) {
00775         good_truth = true;
00776         Float_t *vtx = TrueVtx(mcevent);
00777         Float_t *nu3mom = TrueNuMom(mcevent);
00778         Float_t *targ4vec = Target4Vector(mcevent);
00779         //true info for tree:
00780         is_mc             = 1;
00781         nucleus           = Nucleus(mcevent);
00782         flavour           = Flavour(mcevent);
00783         nooscflavour      = ntpTruth->inunoosc;
00784         initial_state     = Initial_state(mcevent);
00785         had_fs            = HadronicFinalState(mcevent);
00786         process           = IResonance(mcevent); 
00787         cc_nc             = IAction(mcevent);
00788         true_enu          = TrueNuEnergy(mcevent); 
00789         true_pxnu         = nu3mom[0];
00790         true_pynu         = nu3mom[1];
00791         true_pznu         = nu3mom[2];
00792         true_emu          = TrueMuEnergy(mcevent); 
00793         true_elep         = TrueLeptonEnergy(mcevent); 
00794         true_eshw         = TrueShwEnergy(mcevent); 
00795         true_x            = X(mcevent);
00796         true_y            = Y(mcevent);
00797         true_q2           = Q2(mcevent);
00798         true_w2           = W2(mcevent);
00799         true_dircosz      = TrueMuDCosZ(mcevent);
00800         true_dircosneu    = TrueMuDCosNeu(mcevent);
00801         true_vtxx         = vtx[0];
00802         true_vtxy         = vtx[1];
00803         true_vtxz         = vtx[2];
00804         true_etgt         = targ4vec[3];
00805         true_pxtgt        = targ4vec[0];
00806         true_pytgt        = targ4vec[1];
00807         true_pztgt        = targ4vec[2];
00808         true_emfrac       = ntpTruth->emfrac;
00809 
00810         if(LoadTHEvent(j)) {
00811           th_evt_pur      = ntpTHEvent->purity;
00812           th_evt_all      = ntpTHEvent->completeall;
00813           th_evt_slc      = ntpTHEvent->completeslc;
00814         }
00815 
00816         this->SetStdHepVars(mcevent);
00817 
00818         if(gnumi->Status()) {
00819           if(isST) gnumi->GetParent(strecord,mcevent,*nuparent);
00820           else gnumi->GetParent(mcrecord,mcevent,*nuparent);
00821         }
00822 
00823         delete [] vtx;
00824         delete [] nu3mom;
00825         delete [] targ4vec;
00826       }
00827     
00828       //call user functions
00829       this->CallUserFunctions(event);
00830 
00831       //reco info for tree:
00832       reco_vtxx         = ntpEvent->vtx.x;
00833       reco_vtxy         = ntpEvent->vtx.y;
00834       reco_vtxz         = ntpEvent->vtx.z;
00835       evtlength         = ntpEvent->plane.n;
00836       if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00837         evtlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00838       }
00839 
00840       if(nevent>1) {
00841         for(int k=0;k<nevent;k++){
00842           if(k!=j){
00843             Double_t t_sep = TMath::Abs(evtt[k]-evtt[j]);
00844             if(closestevent_t<0 || t_sep<closestevent_t){
00845               closestevent_t = t_sep;
00846               closestevent_s = TMath::Sqrt(TMath::Power(evtx[k]-evtx[j],2)+
00847                                            TMath::Power(evty[k]-evty[j],2)+
00848                                            TMath::Power(evtz[k]-evtz[j],2));
00849               
00850             }
00851           }
00852         }
00853       }
00854 
00855       //index will be -1 if no track/shower present in event
00856       int track_index   = -1;
00857       int shower_index  = -1;
00858       if(LoadLargestTrackFromEvent(j,track_index)) {
00859         if(!LoadShowerAtTrackVertex(j,track_index,shower_index)) {
00860           LoadLargestShowerFromEvent(j,shower_index);
00861         }
00862       }
00863       else LoadLargestShowerFromEvent(j,shower_index);  
00864 
00865       if(good_truth) {
00866         if(LoadTHShower(shower_index)) {
00867           th_shw_pur = ntpTHShower->purity;
00868           th_shw_all = ntpTHShower->completeall;
00869           th_shw_slc = ntpTHShower->completeslc;
00870         }
00871         if(LoadTHTrack(track_index)) {
00872           th_trk_pur = ntpTHTrack->purity;
00873           th_trk_all = ntpTHTrack->completeall;
00874           th_trk_slc = ntpTHTrack->completeslc;
00875         }
00876       }
00877       
00878       //methods should all be protected against index = -1
00879       reco_emu          = RecoMuEnergy(0,track_index);
00880       reco_eshw_linCC   = RecoShwEnergy(shower_index,0);
00881       reco_eshw_wtCC    = RecoShwEnergy(shower_index,1);
00882       reco_eshw_linNC   = RecoShwEnergy(shower_index,2);
00883       reco_eshw_wtNC    = RecoShwEnergy(shower_index,3);
00884       reco_eshw_em      = RecoShwEnergy(shower_index,4);
00885       reco_eshw         = reco_eshw_linCC;
00886       reco_enu          = reco_emu + reco_eshw;
00887       reco_qe_enu       = RecoQENuEnergy(track_index);
00888       reco_qe_q2        = RecoQEQ2(track_index);
00889       if (reco_enu>0) reco_y = reco_eshw/reco_enu;
00890       reco_dircosz      = RecoMuDCosZ(track_index);
00891       reco_dircosneu    = RecoMuDCosNeu(track_index);
00892       is_cev            = IsFidAll(track_index);
00893 
00894       Double_t emfrac0 = 0;
00895       Double_t emfrac1 = 0;
00896       reco_emfrac       = RecoEMFrac(shower_index,0,emfrac0,emfrac1);
00897       reco_emfrac_prob  = RecoEMFrac(shower_index,1,emfrac0,emfrac1);
00898       
00899       //check track is present before using ntpTrack
00900       if(PassTrackCuts(track_index)){ 
00901         trklength         = ntpTrack->plane.n;
00902         trkmomrange       = ntpTrack->momentum.range;
00903         trkqp             = ntpTrack->momentum.qp;
00904         trkeqp            = ntpTrack->momentum.eqp;
00905         Float_t phtrack=ntpTrack->ph.sigcor;
00906         Float_t phevent=ntpEvent->ph.sigcor;
00907         if (phevent>0) {trkphfrac=phtrack/phevent;}
00908         if(ntpTrack->plane.n>0) {
00909           trkphplane      = ntpTrack->ph.sigcor/ntpTrack->plane.n;
00910         }
00911       }
00912       else {
00913         trklength         = 0;
00914         trkmomrange       = 0;
00915         trkqp             = 0;
00916         trkeqp            = 0;
00917         trkphfrac         = 0;
00918         trkphplane        = 0;
00919       }
00920       
00921       if(shower_index!=-1){
00922         shwlength = ntpShower->plane.n; 
00923         Double_t sum_x = 0;
00924         Double_t sum_x2 = 0;
00925         for(int k=0;k<ntpShower->nstrip;k++){
00926           if(!LoadStrip(ntpShower->stp[k])) continue;
00927           Double_t t0 = ntpStrip->time0;
00928           Double_t t1 = ntpStrip->time1;
00929           Double_t avet = 0;
00930           if(t0>0&&t1>0) avet = (t0+t1)/2.;
00931           else if(t0>0) avet = t0;
00932           else if(t1>0) avet = t1;
00933           if(k==0) shwstptimemin = shwstptimemax = avet;
00934           else {
00935             if(avet>shwstptimemax) shwstptimemax = avet;
00936             if(avet<shwstptimemin) shwstptimemin = avet;
00937           }
00938           sum_x +=  (ntpStrip->ph0.pe + ntpStrip->ph1.pe)*avet;
00939           sum_x2 += (ntpStrip->ph0.pe + ntpStrip->ph1.pe)*avet*avet;
00940         }
00941         shwstptimemin -= eventSummary->trigtime;
00942         shwstptimemax -= eventSummary->trigtime;
00943         if(ntpShower->ph.pe>0){
00944           sum_x2 /= ntpShower->ph.pe;
00945           sum_x  /= ntpShower->ph.pe;
00946           shwstptimerms = sum_x2 - sum_x*sum_x;
00947           if(shwstptimerms>0) shwstptimerms = TMath::Sqrt(shwstptimerms);
00948           else shwstptimerms = 0;
00949         }
00950       }
00951 
00952       if(PassBasicCuts(event)) {
00953         pass_basic        = 1;      
00954         pid0              = PID(event,0);
00955         pid1              = PID(event,1);
00956         pid2              = PID(event,2);
00957         if(PassAnalysisCuts(event)) pass = 1;
00958       }
00959 
00960       //fill the tree
00961       tree->Fill();
00962     }
00963   }
00964   delete nuparent;
00965 
00966   save.cd();
00967   pottree->Fill();
00968   pottree->Write();
00969   tree->Write();
00970   //save.Write();
00971   save.Close();
00972 
00973 }

Bool_t MadAnalysis::Do2DFit ( float  n_a = 90.,
float  min_a = 0.5,
float  max_a = 1.0,
float  n_b = 100.,
float  min_b = 0.0,
float  max_b = 0.005,
float  n_c = 1.,
float  f_c = 0. 
)

Definition at line 1524 of file MadAnalysis.cxx.

References BestFit, Chi2Calc, Chi2Surf, ChiMin, DataHist, DataPOT, EShiftErr, MadChi2Calc::GetChi2(), MCEvents, MCPOT, MCSignalHist, NDOF, NumberOfBins, OscillationFunction, Par1MinVal, Par2MinVal, RangeLowerLimit, RangeUpperLimit, and varyX.

01527 {
01528 
01529   //This method is used only for fitting Far Detector spectrum
01530 
01531   if(!OscillationFunction||!DataHist[1]||!MCSignalHist[1]) {
01532     cout << "**ERROR** Did not find one of the following: " << endl;
01533     cout << "Oscillation Function," 
01534          << " Far Detector Reconstructed Data Histogram," 
01535          << " Far Detector Reconstructed MC Histogram" << endl;
01536     cout << "Ensure these are set and try again" << endl;
01537     return false;
01538   }
01539 
01540   //center of last bin of RecoHist
01541   Float_t midupbin = DataHist[1]->GetBinCenter(NumberOfBins);
01542 
01543   //Set up chi2 surface
01544   Float_t Par1Min=min_a;
01545   Float_t Par1Max=max_a;
01546   Int_t Par1Bins=int(n_a);
01547   Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
01548   Float_t Par2Min=min_b;
01549   Float_t Par2Max=max_b;
01550   Int_t Par2Bins=int(n_b);
01551   Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
01552 
01553   Int_t EShiftBins=int(n_c);
01554   Float_t fracEShift=f_c;
01555   Float_t EShiftWidth=2.*f_c/n_c;
01556   
01557   NDOF = NumberOfBins - 2; //two fit params
01558   
01559   Chi2Surf = new TH2F("Chi2Surf",
01560                       "ChiSquare Surface",
01561                       Par1Bins,Par1Min,Par1Max,
01562                       Par2Bins,Par2Min,Par2Max);
01563   
01564   cout << "Starting fit:" << endl;
01565   
01566   for(int i=0;i<Par1Bins;i++){
01567     float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
01568     if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
01569     for(int j=0;j<Par2Bins;j++){
01570       float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
01571       OscillationFunction->SetParameter(varyX[0],par1);
01572       OscillationFunction->SetParameter(varyX[1],par2);
01573       
01574       float ChiSq = 999999; //best ChiSq value in the E shift loop
01575       TH1F *TempBestFit = NULL; //to hold best fit histo for each EShift loop
01576       
01577       //Include a loop to test systematic shifts in reco energy
01578       for(int k=0;k<EShiftBins;k++){
01579         float EShift = 1. - fracEShift + k*EShiftWidth;
01580         
01581         TH1F *temp_sig = new TH1F("temp_sig","temp_sig",
01582                                   NumberOfBins,RangeLowerLimit,
01583                                   RangeUpperLimit);
01584         TH1F *temp_bkg = new TH1F("temp_bkg","temp_bkg",
01585                                   NumberOfBins,RangeLowerLimit,
01586                                   RangeUpperLimit);
01587         
01588         vector<mcev_reweight>::iterator mcbeg = MCEvents[1].begin();
01589         vector<mcev_reweight>::iterator mcend = MCEvents[1].end();
01590         while(mcbeg!=mcend){
01591           Float_t Energy = mcbeg->TrueNuEnergy;
01592           Float_t Ereco  = mcbeg->RecoNuEnergy*EShift;
01593           Float_t weight = mcbeg->Weight;
01594           if(mcbeg->PassTruth){
01595             Double_t oscprob = OscillationFunction->Eval(Energy);
01596             if(oscprob>1.0) oscprob=1.0;
01597             if(Ereco>RangeUpperLimit) Ereco = midupbin;
01598             temp_sig->Fill(Ereco,(1.-oscprob)*weight);
01599           }
01600           else {
01601             if(Ereco>RangeUpperLimit) Ereco = midupbin;
01602             temp_bkg->Fill(Ereco,weight);
01603           }
01604           mcbeg++;
01605         }
01606         
01607         //Chi2 Calculation
01608         float chisq = Chi2Calc->GetChi2(DataHist[1],temp_sig,temp_bkg,
01609                                         MCPOT[1]/DataPOT[1]);
01610         //need to add penalty to chisq depending on sys error on energy scale
01611         if(EShiftErr>0) chisq += TMath::Power((1.-EShift)/EShiftErr,2);
01612         if(chisq<ChiSq) {
01613           delete TempBestFit;
01614           TempBestFit = temp_sig;
01615           TempBestFit->SetName("TempFarBestFit");
01616           ChiSq = chisq;
01617         }
01618         else delete temp_sig;
01619         delete temp_bkg;
01620       }
01621       
01622       if(ChiSq<ChiMin){
01623         delete BestFit[1];
01624         ChiMin=ChiSq;
01625         Par1MinVal=par1;
01626         Par2MinVal=par2;
01627         TempBestFit->Scale(DataPOT[1]/MCPOT[1]);
01628         TempBestFit->SetName("FarBestFit");
01629         BestFit[1] = TempBestFit;
01630       }
01631       else delete TempBestFit;
01632       Chi2Surf->Fill(par1,par2,ChiSq);
01633     }
01634   }
01635   std::cout << "Minimum value of Chi2 is " << ChiMin << std::endl;
01636   std::cout << "at Par1 = " << Par1MinVal << " and Par2 = "
01637             << Par2MinVal << std::endl;
01638   
01639   return true;
01640 }

Bool_t MadAnalysis::Do2DFitNearFar ( float  n_a = 90.,
float  min_a = 0.5,
float  max_a = 1.0,
float  n_b = 100.,
float  min_b = 0.0,
float  max_b = 0.005,
float  n_c = 1.,
float  f_c = 0. 
)

Definition at line 1643 of file MadAnalysis.cxx.

References BestFit, Chi2Calc, Chi2Surf, ChiMin, MCReweight::ComputeWeight(), DataHist, DataPOT, EShiftErr, MadChi2Calc::GetChi2(), MCReweight::Instance(), Munits::m, MCEvents, MCPOT, MCSignalHist, NDOF, NgenBins, NgenErr, NgenMax, NgenMean, NgenMin, NgenName, NumberOfBins, OscillationFunction, Par1MinVal, Par2MinVal, RangeLowerLimit, RangeUpperLimit, Registry::Set(), and varyX.

01646 {
01647 
01648   //Get instance of MCReweight object:
01649   MCReweight &fReweight = MCReweight::Instance();
01650 
01651   //This method is used for simultaneously fitting Near & Far Detector 
01652   //spectra taking into account nuisance parameters
01653 
01654   if(!OscillationFunction||!DataHist[0]||!MCSignalHist[0]) {
01655     cout << "**ERROR** Did not find one of the following: " << endl;
01656     cout << "Oscillation Function," 
01657          << " Reconstructed Data Histograms," 
01658          << " Reconstructed MC Histograms" << endl;
01659     cout << "Ensure these are set and try again" << endl;
01660     return false;
01661   }
01662 
01663   //center of last bin of RecoHist
01664   Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01665 
01666   //Set up chi2 surface
01667   Float_t Par1Min=min_a;
01668   Float_t Par1Max=max_a;
01669   Int_t Par1Bins=int(n_a);
01670   Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
01671   Float_t Par2Min=min_b;
01672   Float_t Par2Max=max_b;
01673   Int_t Par2Bins=int(n_b);
01674   Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
01675 
01676   Int_t EShiftBins=int(n_c);
01677   Float_t fracEShift=f_c;
01678   Float_t EShiftWidth=2.*f_c/n_c;
01679   
01680   NDOF = 2*(NumberOfBins -1) - 2; //Near,Far hists, shape fits, 2 fit params
01681   
01682   Chi2Surf = new TH2F("Chi2Surf",
01683                       "ChiSquare Surface",
01684                       Par1Bins,Par1Min,Par1Max,
01685                       Par2Bins,Par2Min,Par2Max);
01686   
01687   cout << "Starting fit:" << endl;
01688   
01689   for(int i=0;i<Par1Bins;i++){
01690     float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
01691     if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
01692     for(int j=0;j<Par2Bins;j++){
01693       float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
01694       OscillationFunction->SetParameter(varyX[0],par1);
01695       OscillationFunction->SetParameter(varyX[1],par2);
01696       
01697       float ChiSq = 999999; //best ChiSq value in the E shift loop
01698       TH1F *TempBestFit[2]; //to hold best fit histo for each EShift loop
01699       TempBestFit[0] = NULL; TempBestFit[1] = NULL;
01700 
01701       //Include a loop to test systematic shifts in reco energy
01702       for(int k=0;k<EShiftBins;k++){
01703         float EShift = 1. - fracEShift + k*EShiftWidth;
01704         
01705         //Include another few loops to deal with systematics:
01706         //don't need to use all of them, and can add more if necessary
01707         for(int l=0;l<=NgenBins[0];l++){
01708           double NgenPar0 = NgenMin[0] + l*(NgenMax[0] -
01709                                            NgenMin[0])/NgenBins[0];
01710           for(int m=0;m<=NgenBins[1];m++){
01711             double NgenPar1 = NgenMin[1] + m*(NgenMax[1] -
01712                                              NgenMin[1])/NgenBins[1];       
01713             for(int n=0;n<=NgenBins[2];n++){
01714               double NgenPar2 = NgenMin[2] + n*(NgenMax[2] -
01715                                                NgenMin[2])/NgenBins[2];
01716           
01717               Registry rwtconfig;
01718               rwtconfig.Set(NgenName[0].data(),NgenPar0);
01719               rwtconfig.Set(NgenName[1].data(),NgenPar1);
01720               rwtconfig.Set(NgenName[2].data(),NgenPar2);
01721             
01722               TH1F *temp_sig_ND = new TH1F("temp_sig_ND","temp_sig_ND",
01723                                            NumberOfBins,RangeLowerLimit,
01724                                            RangeUpperLimit);
01725               TH1F *temp_bkg_ND = new TH1F("temp_bkg_ND","temp_bkg_ND",
01726                                            NumberOfBins,RangeLowerLimit,
01727                                            RangeUpperLimit);
01728               
01729               vector<mcev_reweight>::iterator mcbeg = MCEvents[0].begin();
01730               vector<mcev_reweight>::iterator mcend = MCEvents[0].end();
01731               while(mcbeg!=mcend){
01732                 Float_t Ereco  = mcbeg->RecoNuEnergy*EShift;
01733                 Float_t weight = mcbeg->Weight;
01734                 Registry event;
01735                 event.Set("event:energy",mcbeg->TrueNuEnergy);
01736                 event.Set("event:iaction",mcbeg->IAction);
01737                 event.Set("event:inu",mcbeg->INu);
01738                 event.Set("event:process",mcbeg->Process);
01739                 event.Set("event:initial_state",mcbeg->InitState);
01740                 event.Set("event:nucleus",mcbeg->Nucleus);
01741                 event.Set("event:q2",mcbeg->KinQ2);
01742                 Float_t mcreweight = fReweight.ComputeWeight(&event,
01743                                                              &rwtconfig);
01744                 
01745                 if(mcbeg->PassTruth){
01746                   if(Ereco>RangeUpperLimit) Ereco = midupbin;
01747                   temp_sig_ND->Fill(Ereco,weight*mcreweight);
01748                 }
01749                 else {
01750                   if(Ereco>RangeUpperLimit) Ereco = midupbin;
01751                   temp_bkg_ND->Fill(Ereco,weight*mcreweight);
01752                 }
01753                 mcbeg++;
01754               }
01755               
01756               TH1F *temp_sig_FD = new TH1F("temp_sig_FD","temp_sig_FD",
01757                                            NumberOfBins,RangeLowerLimit,
01758                                            RangeUpperLimit);
01759               TH1F *temp_bkg_FD = new TH1F("temp_bkg_FD","temp_bkg_FD",
01760                                            NumberOfBins,RangeLowerLimit,
01761                                            RangeUpperLimit);
01762               
01763               mcbeg = MCEvents[1].begin();
01764               mcend = MCEvents[1].end();
01765               while(mcbeg!=mcend){
01766                 Float_t Energy = mcbeg->TrueNuEnergy;
01767                 Float_t Ereco  = mcbeg->RecoNuEnergy*EShift;
01768                 Float_t weight = mcbeg->Weight;
01769                 Registry event;
01770                 event.Set("event:energy",mcbeg->TrueNuEnergy);
01771                 event.Set("event:iaction",mcbeg->IAction);
01772                 event.Set("event:inu",mcbeg->INu);
01773                 event.Set("event:process",mcbeg->Process);
01774                 event.Set("event:initial_state",mcbeg->InitState);
01775                 event.Set("event:nucleus",mcbeg->Nucleus);
01776                 event.Set("event:q2",mcbeg->KinQ2);
01777                 Float_t mcreweight = fReweight.ComputeWeight(&event,
01778                                                              &rwtconfig);
01779                 if(mcbeg->PassTruth){
01780                   Double_t oscprob = OscillationFunction->Eval(Energy);
01781                   if(oscprob>1.0) oscprob=1.0;
01782                   if(Ereco>RangeUpperLimit) Ereco = midupbin;
01783                   temp_sig_FD->Fill(Ereco,(1.-oscprob)*weight*mcreweight);
01784                 }
01785                 else {
01786                   if(Ereco>RangeUpperLimit) Ereco = midupbin;
01787                   temp_bkg_FD->Fill(Ereco,weight*mcreweight);
01788                 }
01789                 mcbeg++;
01790               }
01791               
01792               //Chi2 Calculation
01793               float chisq = Chi2Calc->GetChi2(DataHist[0],
01794                                               temp_sig_ND,temp_bkg_ND,
01795                                               MCPOT[0]/DataPOT[0]);
01796               chisq      += Chi2Calc->GetChi2(DataHist[1],
01797                                               temp_sig_FD,temp_bkg_FD,
01798                                               MCPOT[1]/DataPOT[1]);
01799 
01800               //need to add penalty to chisq depending on sys error
01801               if(EShiftErr>0) chisq += TMath::Power((1.-EShift)/EShiftErr,2);
01802               if(NgenMax[0]!=NgenMin[0])
01803                 chisq += TMath::Power((NgenMean[0]-NgenPar0)/NgenErr[0],2);
01804               if(NgenMax[1]!=NgenMin[1])
01805                 chisq += TMath::Power((NgenMean[1]-NgenPar1)/NgenErr[1],2);
01806               if(NgenMax[2]!=NgenMin[2])
01807                 chisq += TMath::Power((NgenMean[2]-NgenPar2)/NgenErr[2],2);
01808               
01809               if(chisq<ChiSq) {
01810                 delete TempBestFit[0];
01811                 TempBestFit[0] = temp_sig_ND;
01812                 TempBestFit[0]->SetName("TempNearBestFit");
01813                 delete TempBestFit[1];
01814                 TempBestFit[1] = temp_sig_FD;
01815                 TempBestFit[1]->SetName("TempFarBestFit");
01816                 ChiSq = chisq;
01817               }
01818               else {
01819                 delete temp_sig_ND;
01820                 delete temp_sig_FD;
01821               }
01822               delete temp_bkg_ND;
01823               delete temp_bkg_FD;
01824             }
01825           }
01826         }
01827       }
01828       if(ChiSq<ChiMin){
01829         delete BestFit[0];
01830         delete BestFit[1];
01831         ChiMin=ChiSq;
01832         Par1MinVal=par1;
01833         Par2MinVal=par2;
01834         TempBestFit[0]->Scale(DataPOT[0]/MCPOT[0]);
01835         TempBestFit[0]->SetName("NearBestFit");
01836         BestFit[0] = TempBestFit[0];
01837         TempBestFit[1]->Scale(DataPOT[1]/MCPOT[1]);
01838         TempBestFit[1]->SetName("FarBestFit");
01839         BestFit[1] = TempBestFit[1];
01840       }
01841       else {
01842         delete TempBestFit[0];
01843         delete TempBestFit[1];
01844       }
01845       Chi2Surf->Fill(par1,par2,ChiSq);
01846     }
01847   }
01848   std::cout << "Minimum value of Chi2 is " << ChiMin << std::endl;
01849   std::cout << "at Par1 = " << Par1MinVal << " and Par2 = "
01850             << Par2MinVal << std::endl;
01851   
01852   return true;
01853 }

void MadAnalysis::EndJob (  ) 

Definition at line 1856 of file MadAnalysis.cxx.

References BestFit, Chi2Surf, ChiMin, DataBkgdHist, DataHist, DataPOT, DataSignalHist, EShiftErr, MadChi2Calc::GetBinToBinError(), MadChi2Calc::GetBkdXSecError(), MadChi2Calc::GetCalcType(), GetChi2Calc(), MadChi2Calc::GetNormalisationError(), MCBkgdHist, MCPOT, MCSignalHist, NDOF, numDataEvts, numMCEvts, NumOscPars, OscillationParameters, Par1MinVal, Par2MinVal, signalFracInData, and tag.

Referenced by ~MadAnalysis().

01857 {
01858   
01859   std::string outputFileName = "MadFile_" + tag + ".root";
01860   TFile *outputFile = new TFile(outputFileName.c_str(),"RECREATE");
01861   outputFile->cd();
01862   for(int i=0;i<2;i++){
01863     if(MCSignalHist[i]) MCSignalHist[i]->Write();
01864     if(MCBkgdHist[i]) MCBkgdHist[i]->Write();
01865     if(DataHist[i]) DataHist[i]->Write();
01866     if(DataSignalHist[i]) DataSignalHist[i]->Write();
01867     if(DataBkgdHist[i]) DataBkgdHist[i]->Write();
01868     if(BestFit[i]) BestFit[i]->Write();
01869   }
01870   if(Chi2Surf) Chi2Surf->Write();
01871   TTree *varTree = new TTree("Params","Params");
01872   varTree->Branch("DataPOTNear",&DataPOT[0],"DataPOTNear/F",32000);
01873   varTree->Branch("DataPOTFar",&DataPOT[1],"DataPOTFar/F",32000);
01874   varTree->Branch("numDataEvtsNear",&numDataEvts[0],"numDataEvtsNear/I",32000);
01875   varTree->Branch("numDataEvtsFar",&numDataEvts[1],"numDataEvtsFar/I",32000);
01876   varTree->Branch("signalFracInDataNear",&signalFracInData[0],
01877                   "signalFracInDataNear/F",32000);
01878   varTree->Branch("signalFracInDataFar",&signalFracInData[1],
01879                   "signalFracInDataFar/F",32000);
01880   
01881   varTree->Branch("MCPOTNear",&MCPOT[0],"MCPOTNear/F",32000);
01882   varTree->Branch("MCPOTFar",&MCPOT[1],"MCPOTFar/F",32000);
01883   varTree->Branch("numMCEvtsNear",&numMCEvts[0],"numMCEvtsNear/I",32000);
01884   varTree->Branch("numMCEvtsFar",&numMCEvts[1],"numMCEvtsFar/I",32000);
01885   
01886   varTree->Branch("NumOscPars",&NumOscPars,"NumOscPars/I",32000);
01887   varTree->Branch("OscillationParameters",OscillationParameters,
01888                   "OscillationParameters[10]/D",32000);
01889   
01890   varTree->Branch("ChiMin",&ChiMin,"ChiMin/F",32000);
01891   varTree->Branch("NDOF",&NDOF,"NDOF/F",32000);
01892   varTree->Branch("Par1MinVal",&Par1MinVal,"Par1MinVal/F",32000);
01893   varTree->Branch("Par2MinVal",&Par2MinVal,"Par2MinVal/F",32000);
01894 
01895   varTree->Branch("EShiftErr",&EShiftErr,"EShiftErr/F",32000);
01896 
01897   double NormErr = this->GetChi2Calc()->GetNormalisationError();
01898   varTree->Branch("NormErr",&NormErr,"NormErr/D",32000);
01899   double BinErr = this->GetChi2Calc()->GetBinToBinError();
01900   varTree->Branch("BinErr",&BinErr,"BinErr/D",32000);
01901   double CalcType = this->GetChi2Calc()->GetCalcType();
01902   varTree->Branch("CalcType",&CalcType,"CalcType/D",32000);
01903   double BkdXSecErr = this->GetChi2Calc()->GetBkdXSecError();
01904   varTree->Branch("BkdXSecErr",&BkdXSecErr,"BkdXSecErr/D",32000);
01905   
01906   varTree->Fill();
01907   varTree->Write();
01908   outputFile->Close();
01909 
01910 }

MadChi2Calc* MadAnalysis::GetChi2Calc (  )  [inline]

Definition at line 160 of file MadAnalysis.h.

References Chi2Calc.

Referenced by EndJob().

00160 {return Chi2Calc;} //so user can set options

TH1F* MadAnalysis::GetDataHist ( int  det  )  [inline]

Definition at line 165 of file MadAnalysis.h.

References DataHist.

00165                              { if(det!=0&&det!=1) return 0; 
00166                                return DataHist[det]; }

float MadAnalysis::GetDataPOT ( int  det = 1  )  [inline]

Definition at line 168 of file MadAnalysis.h.

References DataPOT, and det.

00168 { return DataPOT[det]; }

int MadAnalysis::GetNumDataEvents ( int  det = 1  )  [inline]

Definition at line 167 of file MadAnalysis.h.

References det, and numDataEvts.

00167 { return numDataEvts[det]; }

float MadAnalysis::GetSignalFracInData ( int  det = 1  )  [inline]

Definition at line 169 of file MadAnalysis.h.

References det, and signalFracInData.

00169 { return signalFracInData[det]; }

virtual Float_t MadAnalysis::GetWeight ( Int_t  event = 0  )  [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, and MadUserAnalysis.

Definition at line 115 of file MadAnalysis.h.

Referenced by RecoExperiment(), RecoMC(), and RecoMCExperiment().

00115 {if(event>=0) return 1.0; return 0;}

void MadAnalysis::NoOutFile (  )  [inline]

Definition at line 171 of file MadAnalysis.h.

References writeOut.

00171 { writeOut = false;}  //don't write output file

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

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, MadDpAnalysis, MadEdAnalysis, MadPIDAnalysis, MadTestAnalysis, and MadUserAnalysis.

Definition at line 104 of file MadAnalysis.h.

Referenced by CreateANtpPAN(), RecoExperiment(), RecoMC(), and RecoMCExperiment().

00104                                                  { if(event>=0) return true; 
00105                                                    return false; }

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

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, MadDpAnalysis, MadEdAnalysis, MadPIDAnalysis, MadTestAnalysis, and MadUserAnalysis.

Definition at line 101 of file MadAnalysis.h.

Referenced by CreateANtpPAN(), and MadMKAnalysis::CreatePAN().

00101                                               { if(event>=0) return true; 
00102                                                 return false; }

virtual Bool_t MadAnalysis::PassTruthSignal ( Int_t  mcevent = 0  )  [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, MadDpAnalysis, MadPIDAnalysis, MadTestAnalysis, and MadUserAnalysis.

Definition at line 97 of file MadAnalysis.h.

References MadBase::LoadTruth().

Referenced by RecoMC(), and RecoMCExperiment().

00097                                                   { 
00098                                        if(!LoadTruth(mcevent)) return false;
00099                                        return true; }

virtual Float_t MadAnalysis::PID ( Int_t  event = 0,
Int_t  method = 0 
) [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, MadDpAnalysis, MadEdAnalysis, MadPIDAnalysis, MadTestAnalysis, and MadUserAnalysis.

Definition at line 107 of file MadAnalysis.h.

Referenced by CreateANtpPAN().

00107                                                     { 
00108                                        if(method>=0 && event>=0) return 1.0;
00109                                        return 0.0; }

virtual Float_t MadAnalysis::RecoAnalysisEnergy ( Int_t  event = 0  )  [inline, protected, virtual]

Reimplemented in MadCBSQEAnalysis, MadCluAnalysis, MadDpAnalysis, MadPIDAnalysis, MadTestAnalysis, and MadUserAnalysis.

Definition at line 111 of file MadAnalysis.h.

References NtpSRStripPulseHeight::gev, MadBase::LoadEvent(), MadBase::ntpEvent, and NtpSREvent::ph.

Referenced by RecoExperiment(), RecoMC(), and RecoMCExperiment().

00111                                                     { 
00112                                        if(!LoadEvent(event)) return 0;
00113                                        return ntpEvent->ph.gev; }

void MadAnalysis::RecoExperiment ( int  startnum,
double  NearExpectedNormToPOT,
double  FarExpectedNormToPOT 
)

Definition at line 1150 of file MadAnalysis.cxx.

References DataHist, DataPOT, MadBase::eventSummary, VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), GetWeight(), MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpHeader, NumberOfBins, numDataEvts, PassAnalysisCuts(), RangeLowerLimit, RangeUpperLimit, RecoAnalysisEnergy(), and signalFracInData.

01153 {
01154 
01155   cout << "Reconstructing Data Energy Spectrum" << endl;
01156   if(!DataHist[0]){
01157     DataHist[0] = 
01158       new TH1F("NearDataHist",
01159                "Near Detector Reconstructed Energy Spectrum",
01160                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01161     DataHist[0]->SetXTitle("E_{reco} (GeV)");
01162     DataHist[0]->SetYTitle("Event Rate");
01163     DataHist[0]->Sumw2();
01164   }
01165   else {
01166     DataHist[0]->Reset();
01167     numDataEvts[0] = 0;
01168     signalFracInData[0] = 0;
01169   }
01170 
01171   if(!DataHist[1]){
01172     DataHist[1] = 
01173       new TH1F("FarDataHist",
01174                "Far Detector Oscillated Reconstructed Energy Spectrum",
01175                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01176     DataHist[1]->SetXTitle("E_{reco} (GeV)");
01177     DataHist[1]->SetYTitle("Event Rate");
01178     DataHist[1]->Sumw2();
01179   }
01180   else {
01181     DataHist[1]->Reset();
01182     numDataEvts[1] = 0;
01183     signalFracInData[1] = 0;
01184   }
01185 
01186   Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01187 
01188   int i = startnum;
01189   while(i<=Nentries){ //while there are more entries in the TChains...
01190     if(!GetEntry(i++)) continue;
01191 
01192     if((i-startnum)%500==0) 
01193       std::cout << 100*float(i-startnum)/float(Nentries-startnum) 
01194                 << "% completed" << std::endl;
01195     
01196     //get which detector this snarl is from
01197     int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01198 
01199     //loop over all events in snarl
01200     for(int j=0;j<eventSummary->nevent;j++){      
01201       numDataEvts[theDetector]+=1; //increment for each event reconstructed
01202       if(PassAnalysisCuts(j)){ //if event passes cuts
01203         signalFracInData[theDetector]+=1; //increment "signal" event counter
01204         float ereco = RecoAnalysisEnergy(j); //get neutrino energy
01205         if(ereco<RangeUpperLimit) //check it is in histo range
01206           DataHist[theDetector]->Fill(ereco,GetWeight()); //fill
01207         else DataHist[theDetector]->Fill(midupbin,GetWeight());
01208       }
01209     }
01210   }
01211   
01212   signalFracInData[0]/=float(numDataEvts[0]);
01213   signalFracInData[1]/=float(numDataEvts[1]);
01214   DataPOT[0]=numDataEvts[0]*NearExpectedNormToPOT;
01215   DataPOT[1]=numDataEvts[1]*FarExpectedNormToPOT;
01216 
01217 }

void MadAnalysis::RecoMC ( int  startnum,
double  NearPOT,
double  FarPOT 
)

Definition at line 1365 of file MadAnalysis.cxx.

References MadBase::eventSummary, FARPOTTONUMEVT, VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), GetWeight(), MadQuantities::IAction(), MadQuantities::Initial_state(), MadQuantities::INu(), MadQuantities::IResonance(), MadBase::LoadTruthForRecoTH(), MCBkgdHist, MCEvents, MCPOT, MCSignalHist, NEARPOTTONUMEVT, MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpHeader, MadQuantities::Nucleus(), NumberOfBins, numDataEvts, numMCEvts, PassAnalysisCuts(), PassTruthSignal(), MadQuantities::Q2(), RangeLowerLimit, RangeUpperLimit, RecoAnalysisEnergy(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueShwEnergy(), MadQuantities::W2(), MadQuantities::X(), and MadQuantities::Y().

01366 {
01367   cout << "Reconstructing MC Event Sample" << endl;
01368 
01369   if(!MCSignalHist[0]){
01370     MCSignalHist[0] = new TH1F("NearMCSignalHist",
01371       "Near Detector Reconstructed Unoscillated Energy Spectrum (Signal only)",
01372                                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01373     MCSignalHist[0]->SetXTitle("E_{reco} (GeV)");
01374     MCSignalHist[0]->SetYTitle("Event Rate");
01375     MCSignalHist[0]->Sumw2();
01376     
01377     MCBkgdHist[0] = new TH1F("NearMCBkgdHist",              
01378     "Near Detector Reconstructed Unoscillated Energy Spectrum (Bkd only)",
01379                              NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01380     MCBkgdHist[0]->SetXTitle("E_{reco} (GeV)");
01381     MCBkgdHist[0]->SetYTitle("Event Rate");
01382     MCBkgdHist[0]->Sumw2();
01383   }
01384   else {
01385     MCSignalHist[0]->Reset();
01386     MCBkgdHist[0]->Reset();
01387     MCPOT[0] = 0;
01388     numMCEvts[0] = 0;
01389   }
01390 
01391   if(!MCSignalHist[1]){
01392     MCSignalHist[1] = new TH1F("FarMCSignalHist",
01393       "Far Detector Reconstructed Unoscillated Energy Spectrum (Signal only)",
01394                                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01395     MCSignalHist[1]->SetXTitle("E_{reco} (GeV)");
01396     MCSignalHist[1]->SetYTitle("Event Rate");
01397     MCSignalHist[1]->Sumw2();
01398     
01399     MCBkgdHist[1] = new TH1F("FarMCBkgdHist",              
01400     "Far Detector Reconstructed Unoscillated Energy Spectrum (Bkd only)",
01401                              NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01402     MCBkgdHist[1]->SetXTitle("E_{reco} (GeV)");
01403     MCBkgdHist[1]->SetYTitle("Event Rate");
01404     MCBkgdHist[1]->Sumw2();
01405   }
01406   else {
01407     MCSignalHist[1]->Reset();
01408     MCBkgdHist[1]->Reset();
01409     MCPOT[1] = 0;
01410     numMCEvts[1] = 0;
01411   }
01412 
01413 
01414   Float_t midupbin = MCSignalHist[0]->GetBinCenter(NumberOfBins);
01415   
01416   MCPOT[0] = NearPOT; MCPOT[1] = FarPOT;
01417   numMCEvts[0] = int(NEARPOTTONUMEVT*NearPOT);
01418   numMCEvts[1] = int(FARPOTTONUMEVT*FarPOT);
01419   int counter[2] = {0};
01420   int i = startnum;
01421   
01422   while((counter[0]<numMCEvts[0]||counter[1]<numMCEvts[1])&&i<=Nentries){
01423     if(!GetEntry(i++)) continue;
01424 
01425     //get which detector this snarl is from
01426     int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01427     
01428     //loop over all events in snarl
01429     for(int j=0;j<eventSummary->nevent;j++){
01430     
01431       //if no good truth candidate, we loose event
01432       int mcevent = -1;
01433       if(LoadTruthForRecoTH(j,mcevent)){ 
01434         
01435         if(PassTruthSignal(mcevent)) counter[theDetector]+=1;
01436         
01437         if(counter[theDetector]%1000==0) 
01438           std::cout << 100*(float(counter[theDetector])
01439                             /float(numDataEvts[theDetector])) 
01440                     << "% completed for Detector " << theDetector 
01441                     << std::endl;
01442           
01443         if(PassAnalysisCuts(j)){
01444           float emu = 0;
01445           float eshw = 0;
01446           float mudcosz = 0;
01447           float ereco = RecoAnalysisEnergy(j);
01448           if(PassTruthSignal(mcevent)) { //expect oscillations from these
01449             if(ereco<RangeUpperLimit) //according to OscillationFunction
01450               MCSignalHist[theDetector]->Fill(ereco,GetWeight());
01451             else MCSignalHist[theDetector]->Fill(midupbin,GetWeight());
01452           }
01453           else {
01454             if(ereco<RangeUpperLimit) 
01455               MCBkgdHist[theDetector]->Fill(ereco,GetWeight());
01456             else MCBkgdHist[theDetector]->Fill(midupbin,GetWeight());
01457           }
01458           //fill a struct with all the needed info for event reweighting
01459           mcev_reweight mmei = {TrueNuEnergy(mcevent),TrueMuEnergy(mcevent),
01460                                 TrueShwEnergy(mcevent),TrueMuDCosZ(mcevent),
01461                                 ereco,emu,eshw,mudcosz,GetWeight(),
01462                                 PassTruthSignal(mcevent),INu(mcevent),
01463                                 IAction(mcevent),IResonance(mcevent),
01464                                 Initial_state(mcevent),Nucleus(mcevent),
01465                                 X(mcevent),Y(mcevent),Q2(mcevent),W2(mcevent)};
01466           MCEvents[theDetector].push_back(mmei);
01467         }
01468       }         
01469     }
01470   }
01471 
01472   if(counter[0]<numMCEvts[0]) numMCEvts[0]=counter[0];
01473   MCPOT[0]=numMCEvts[0]/NEARPOTTONUMEVT;
01474   if(counter[1]<numMCEvts[1]) numMCEvts[1]=counter[1];
01475   MCPOT[1]=numMCEvts[1]/FARPOTTONUMEVT;
01476 
01477 }

void MadAnalysis::RecoMCExperiment ( int  startnum,
double  NearPOT,
double  FarPOT 
)

Definition at line 1220 of file MadAnalysis.cxx.

References DataBkgdHist, DataHist, DataPOT, DataSignalHist, MadBase::eventSummary, FARPOTTONUMEVT, VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), GetWeight(), MadBase::LoadTruthForRecoTH(), NEARPOTTONUMEVT, MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpHeader, NumberOfBins, numDataEvts, OscillationFunction, PassAnalysisCuts(), PassTruthSignal(), RangeLowerLimit, RangeUpperLimit, RecoAnalysisEnergy(), signalFracInData, and MadQuantities::TrueNuEnergy().

01221 {
01222 
01223   cout << "Reconstructing MC Data Energy Spectrum" << endl;
01224   
01225   if(!DataHist[0]){
01226     DataHist[0] = 
01227       new TH1F("NearDataHist",
01228                "Near Detector Reconstructed Energy Spectrum",
01229                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01230     DataHist[0]->SetXTitle("E_{reco} (GeV)");
01231     DataHist[0]->SetYTitle("Event Rate");
01232     DataHist[0]->Sumw2();
01233     
01234     DataSignalHist[0] = 
01235       new TH1F("NearDataSignalHist",
01236                "Near Detector Reconstructed Energy Spectrum (Signal only)",
01237                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01238     DataSignalHist[0]->SetXTitle("E_{reco} (GeV)");
01239     DataSignalHist[0]->SetYTitle("Event Rate");
01240     DataSignalHist[0]->Sumw2();
01241     
01242     DataBkgdHist[0] = 
01243       new TH1F("NearDataBkgdHist",
01244                "Near Detector Reconstructed Energy Spectrum (Background only)",
01245                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01246     DataBkgdHist[0]->SetXTitle("E_{reco} (GeV)");
01247     DataBkgdHist[0]->SetYTitle("Event Rate");
01248     DataBkgdHist[0]->Sumw2();
01249   }
01250   else {
01251     DataHist[0]->Reset();
01252     DataSignalHist[0]->Reset();
01253     DataBkgdHist[0]->Reset();
01254     numDataEvts[0] = 0;
01255     signalFracInData[0] = 0;
01256   }
01257   
01258   if(!DataHist[1]){
01259     DataHist[1] = 
01260       new TH1F("FarDataHist",
01261                "Far Detector Oscillated Reconstructed Energy Spectrum",
01262                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01263     DataHist[1]->SetXTitle("E_{reco} (GeV)");
01264     DataHist[1]->SetYTitle("Event Rate");
01265     DataHist[1]->Sumw2();
01266     
01267     DataSignalHist[1] = new TH1F("FarDataSignalHist",
01268     "Far Detector Oscillated Reconstructed Energy Spectrum (Signal only)",
01269                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01270     DataSignalHist[1]->SetXTitle("E_{reco} (GeV)");
01271     DataSignalHist[1]->SetYTitle("Event Rate");
01272     DataSignalHist[1]->Sumw2();
01273     
01274     DataBkgdHist[1] = new TH1F("FarDataBkgdHist",
01275     "Far Detector Oscillated Reconstructed Energy Spectrum (Background only)",
01276                                NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01277     DataBkgdHist[1]->SetXTitle("E_{reco} (GeV)");
01278     DataBkgdHist[1]->SetYTitle("Event Rate");
01279     DataBkgdHist[1]->Sumw2();
01280   }
01281   else {
01282     DataHist[1]->Reset();
01283     DataSignalHist[1]->Reset();
01284     DataBkgdHist[1]->Reset();
01285     numDataEvts[1] = 0;
01286     signalFracInData[1] = 0;
01287   }
01288 
01289   //center of last bin of RecoHist
01290   Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01291 
01292   DataPOT[0] = NearPOT;
01293   DataPOT[1] = FarPOT;
01294   numDataEvts[0] = int(NEARPOTTONUMEVT*NearPOT);
01295   numDataEvts[1] = int(FARPOTTONUMEVT*FarPOT);
01296 
01297   int counter[2] = {0};
01298   
01299   int i = startnum; 
01300   while((counter[0]<numDataEvts[0]||counter[1]<numDataEvts[1])&&i<=Nentries){
01301     if(!GetEntry(i++)) continue;
01302 
01303     //get which detector this snarl is from
01304     int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01305     
01306     //loop over all events in snarl
01307     for(int j=0;j<eventSummary->nevent;j++){
01308       //if no good truth candidate, we loose event
01309       int mcevent = -1;
01310       if(LoadTruthForRecoTH(j,mcevent)){ 
01311 
01312         if(PassTruthSignal(mcevent)) counter[theDetector]+=1;
01313           
01314         if(counter[theDetector]%1000==0) 
01315           std::cout << 100*(float(counter[theDetector])
01316                             /float(numDataEvts[theDetector])) 
01317                     << "% completed for Detector " << theDetector 
01318                     << std::endl;
01319         
01320         if(PassAnalysisCuts(j)){
01321           float ereco = RecoAnalysisEnergy(j);
01322           if(PassTruthSignal(mcevent)) {      //expect oscillations from these
01323             signalFracInData[theDetector]+=1; //events according to OscFunc
01324             
01325             double oscprob = 0;
01326             if(theDetector==1) 
01327               OscillationFunction->Eval(TrueNuEnergy(mcevent));
01328             
01329             if(ereco<RangeUpperLimit) {
01330               DataHist[theDetector]->Fill(ereco,(1.-oscprob)*GetWeight());
01331               DataSignalHist[theDetector]->Fill(ereco,
01332                                                 (1.-oscprob)*GetWeight());
01333             }
01334             else {
01335               DataHist[theDetector]->Fill(midupbin,(1.-oscprob)*GetWeight());
01336               DataSignalHist[theDetector]->Fill(midupbin,
01337                                                 (1.-oscprob)*GetWeight());
01338             }
01339           }
01340           else {  //don't expect oscillations
01341             if(ereco<RangeUpperLimit) {
01342               DataHist[theDetector]->Fill(ereco,GetWeight());
01343               DataBkgdHist[theDetector]->Fill(ereco,GetWeight());
01344             }
01345             else {
01346               DataHist[theDetector]->Fill(midupbin,GetWeight());
01347               DataBkgdHist[theDetector]->Fill(midupbin,GetWeight());
01348             }
01349           }
01350         }
01351       }
01352     }
01353   }
01354   
01355   if(counter[0]<numDataEvts[0]) numDataEvts[0]=counter[0];
01356   DataPOT[0]=numDataEvts[0]/NEARPOTTONUMEVT;
01357   signalFracInData[0]/=float(numDataEvts[0]);
01358   if(counter[1]<numDataEvts[1]) numDataEvts[1]=counter[1];
01359   DataPOT[1]=numDataEvts[1]/FARPOTTONUMEVT;
01360   signalFracInData[1]/=float(numDataEvts[1]);
01361 
01362 }

void MadAnalysis::ReInit ( JobC *  ,
string  ,
int   
)

Definition at line 195 of file MadAnalysis.cxx.

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

00195                                                        {
00196   
00197   record = 0;
00198   strecord = 0;
00199   emrecord = 0;
00200   mcrecord = 0;
00201   threcord = 0;
00202   Clear();
00203   isMC = true;
00204   isTH = true;
00205   isEM = true;
00206   Nentries = entries;
00207   whichSource = 1;
00208   jcPath = path;
00209   fJC = j;
00210 
00211 }

void MadAnalysis::ReInit ( TChain *  chainSR = 0,
TChain *  chainMC = 0,
TChain *  chainTH = 0,
TChain *  chainEM = 0 
)

Definition at line 187 of file MadAnalysis.cxx.

References MadBase::InitChain().

00188                                                          {
00189   
00190   InitChain(chainSR,chainMC,chainTH,chainEM);
00191 
00192 }

int MadAnalysis::SetDataInfo ( TH1F *  hist = 0,
int  numev = 0,
float  pot = 0,
float  sigfrac = 0,
int  det = 1 
)

Definition at line 229 of file MadAnalysis.cxx.

References DataHist, DataPOT, numDataEvts, and signalFracInData.

00230                                      {
00231 
00232   if(hist) {
00233     DataHist[det] = new TH1F(*hist);  
00234     numDataEvts[det] = numev;
00235     DataPOT[det] = pot;
00236     signalFracInData[det] = sigfrac;    
00237     return 1;
00238   }
00239   return 0;
00240 
00241 }

void MadAnalysis::SetEShiftErr ( float  err  )  [inline]

Definition at line 161 of file MadAnalysis.h.

References EShiftErr.

00161 {EShiftErr = err;}

void MadAnalysis::SetFileNameTag ( std::string   ) 

Definition at line 214 of file MadAnalysis.cxx.

References tag.

00215 {
00216   tag = t; 
00217 }

void MadAnalysis::SetHistBookInfo ( int  ,
float  ,
float   
)

Definition at line 220 of file MadAnalysis.cxx.

References NumberOfBins, RangeLowerLimit, and RangeUpperLimit.

00220                                                             {
00221 
00222   NumberOfBins = bins;
00223   RangeLowerLimit = low;
00224   RangeUpperLimit = up;
00225   
00226 }

void MadAnalysis::SetNeugenReweightInfo ( Int_t  ,
Int_t *  ,
Float_t *  ,
Float_t *  ,
Float_t *  ,
Float_t *  ,
std::string *   
)

Definition at line 1511 of file MadAnalysis.cxx.

References NgenBins, NgenErr, NgenMax, NgenMean, NgenMin, and NgenName.

01514 {
01515   if(n>3) n = 3; //can only handle 3 right now
01516   for(int i=0;i<n;i++){
01517     NgenBins[i] = bins[i]; NgenMin[i]  = min[i];
01518     NgenMax[i]  = max[i];  NgenMean[i] = mean[i];
01519     NgenErr[i]  = err[i];  NgenName[i] = name[i];
01520   }
01521 }

void MadAnalysis::SetOscillationFunction ( TF1 *  ,
double *   
)

Definition at line 1497 of file MadAnalysis.cxx.

References NumOscPars, OscillationFunction, and OscillationParameters.

01498 {
01499 
01500   OscillationParameters = params;
01501   OscillationFunction = new TF1(*form);
01502   NumOscPars = OscillationFunction->GetNpar();
01503   
01504   for(int i=0;i<NumOscPars;i++){
01505     OscillationFunction->SetParameter(i,OscillationParameters[i]);
01506   }
01507 
01508 }

void MadAnalysis::SetOscillationFunction ( Double_t(*)(Double_t *, Double_t *)  fcn,
Float_t  ,
Float_t  ,
int  ,
double *   
)

Definition at line 1480 of file MadAnalysis.cxx.

References NumOscPars, OscillationFunction, and OscillationParameters.

01483 {
01484   
01485   OscillationParameters = params;
01486   OscillationFunction = new TF1("OscFunc",fcn,lowend,
01487                                 highend,numpars);  
01488   NumOscPars = numpars;
01489 
01490   for(int i=0;i<NumOscPars;i++){
01491     OscillationFunction->SetParameter(i,OscillationParameters[i]);
01492   }
01493   
01494 }

void MadAnalysis::VaryFitParam ( int  ix,
int  vx 
) [inline]

Definition at line 147 of file MadAnalysis.h.

References varyX.

00147 {varyX[ix] = vx;}

Member Data Documentation

TH1F* MadAnalysis::BestFit[2] [protected]

Definition at line 82 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), MadAnalysis(), and ~MadAnalysis().

MadChi2Calc* MadAnalysis::Chi2Calc [protected]

Definition at line 86 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), GetChi2Calc(), MadAnalysis(), and ~MadAnalysis().

TH2F* MadAnalysis::Chi2Surf [protected]

Definition at line 87 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), MadAnalysis(), and ~MadAnalysis().

float MadAnalysis::ChiMin [protected]

Definition at line 88 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), and MadAnalysis().

TH1F* MadAnalysis::DataBkgdHist[2] [protected]

Definition at line 65 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), RecoMCExperiment(), and ~MadAnalysis().

TH1F* MadAnalysis::DataHist[2] [protected]

Definition at line 60 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), GetDataHist(), MadAnalysis(), RecoExperiment(), RecoMCExperiment(), SetDataInfo(), and ~MadAnalysis().

float MadAnalysis::DataPOT[2] [protected]

Definition at line 62 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), GetDataPOT(), MadAnalysis(), RecoExperiment(), RecoMCExperiment(), and SetDataInfo().

TH1F* MadAnalysis::DataSignalHist[2] [protected]

Definition at line 61 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), RecoMCExperiment(), and ~MadAnalysis().

float MadAnalysis::EShiftErr [protected]

Definition at line 94 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), MadAnalysis(), and SetEShiftErr().

TH1F* MadAnalysis::MCBkgdHist[2] [protected]

Definition at line 59 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), RecoMC(), and ~MadAnalysis().

vector<mcev_reweight> MadAnalysis::MCEvents[2] [protected]

Definition at line 67 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), and RecoMC().

float MadAnalysis::MCPOT[2] [protected]

Definition at line 68 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), MadAnalysis(), and RecoMC().

TH1F* MadAnalysis::MCSignalHist[2] [protected]

Definition at line 58 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), MadAnalysis(), RecoMC(), and ~MadAnalysis().

float MadAnalysis::NDOF [protected]

Definition at line 89 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), and MadAnalysis().

Int_t MadAnalysis::NgenBins[3] [protected]

Definition at line 71 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

Double_t MadAnalysis::NgenErr[3] [protected]

Definition at line 75 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

Double_t MadAnalysis::NgenMax[3] [protected]

Definition at line 73 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

Double_t MadAnalysis::NgenMean[3] [protected]

Definition at line 74 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

Double_t MadAnalysis::NgenMin[3] [protected]

Definition at line 72 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

std::string MadAnalysis::NgenName[3] [protected]

Definition at line 76 of file MadAnalysis.h.

Referenced by Do2DFitNearFar(), MadAnalysis(), and SetNeugenReweightInfo().

int MadAnalysis::NumberOfBins [protected]

Definition at line 53 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), RecoExperiment(), RecoMC(), RecoMCExperiment(), and SetHistBookInfo().

int MadAnalysis::numDataEvts[2] [protected]

Definition at line 63 of file MadAnalysis.h.

Referenced by EndJob(), GetNumDataEvents(), MadAnalysis(), RecoExperiment(), RecoMC(), RecoMCExperiment(), and SetDataInfo().

int MadAnalysis::numMCEvts[2] [protected]

Definition at line 69 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), and RecoMC().

int MadAnalysis::NumOscPars [protected]

Definition at line 79 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), and SetOscillationFunction().

TF1* MadAnalysis::OscillationFunction [protected]

Definition at line 78 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), RecoMCExperiment(), SetOscillationFunction(), and ~MadAnalysis().

double* MadAnalysis::OscillationParameters [protected]

Definition at line 80 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), and SetOscillationFunction().

float MadAnalysis::Par1MinVal [protected]

Definition at line 90 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), and MadAnalysis().

float MadAnalysis::Par2MinVal [protected]

Definition at line 91 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), EndJob(), and MadAnalysis().

float MadAnalysis::RangeLowerLimit [protected]

Definition at line 54 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), RecoExperiment(), RecoMC(), RecoMCExperiment(), and SetHistBookInfo().

float MadAnalysis::RangeUpperLimit [protected]

Definition at line 55 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), RecoExperiment(), RecoMC(), RecoMCExperiment(), and SetHistBookInfo().

float MadAnalysis::signalFracInData[2] [protected]

Definition at line 64 of file MadAnalysis.h.

Referenced by EndJob(), GetSignalFracInData(), MadAnalysis(), RecoExperiment(), RecoMCExperiment(), and SetDataInfo().

std::string MadAnalysis::tag [protected]

Definition at line 84 of file MadAnalysis.h.

Referenced by EndJob(), MadAnalysis(), and SetFileNameTag().

int* MadAnalysis::varyX [protected]

Definition at line 81 of file MadAnalysis.h.

Referenced by Do2DFit(), Do2DFitNearFar(), MadAnalysis(), VaryFitParam(), and ~MadAnalysis().

Bool_t MadAnalysis::writeOut [protected]

Definition at line 120 of file MadAnalysis.h.

Referenced by MadAnalysis(), NoOutFile(), and ~MadAnalysis().

The documentation for this class was generated from the following files:
Generated on Tue Jun 18 22:01:15 2013 for loon by  doxygen 1.4.7