MadTestAnalysis Class Reference

#include <MadTestAnalysis.h>

Inheritance diagram for MadTestAnalysis:
MadAnalysis MadQuantities MadBase

List of all members.

Public Member Functions

 MadTestAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadTestAnalysis (JobC *, string, int)
 ~MadTestAnalysis ()
void MakeMyFile (std::string, int)
void ReadPIDFile (std::string)
void CreatePAN (std::string, Int_t)

Static Public Member Functions

static bool InFidVol (const Detector::Detector_t &det, const Float_t &x, const Float_t &y, const Float_t &z)

Public Attributes

MadNsID nsid
MadDpID dpid

Protected Member Functions

Bool_t PassTruthSignal (Int_t event=0)
Bool_t PassAnalysisCuts (Int_t event=0)
Float_t PID (Int_t event=0, Int_t method=0)
Bool_t PassBasicCuts (Int_t event=0)
Bool_t PassFid (Int_t event=0)
Float_t RecoAnalysisEnergy (Int_t event=0)
AnnInputBlock AnnVar (Int_t event=0)
Double_t GetAnnPid (AnnInputBlock anninput, Int_t det, Int_t tar, Int_t fid, Int_t prior)
Double_t Sigmoid (Double_t x)
Float_t SingleTimeFrame (Int_t snarlentry, Int_t nentries)

Protected Attributes

TFile * fLikeliFile
TH1F * fLikeliHist [6]

Detailed Description

Definition at line 10 of file MadTestAnalysis.h.


Constructor & Destructor Documentation

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

Definition at line 25 of file MadTestAnalysis.cxx.

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

00027 {
00028 
00029   if(!chainSR) {
00030     record = 0;
00031     strecord = 0;
00032     emrecord = 0;
00033     mcrecord = 0;
00034     threcord = 0;
00035     Clear();
00036     whichSource = -1;
00037     isMC = true;
00038     isTH = true;
00039     isEM = true;
00040     Nentries = -1;
00041     return;
00042   }
00043   
00044   InitChain(chainSR,chainMC,chainTH,chainEM);
00045   whichSource = 0;
00046   fLikeliFile = NULL;
00047   for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00048 
00049 }

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

Definition at line 52 of file MadTestAnalysis.cxx.

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

00053 {
00054   record = 0;
00055   strecord = 0;
00056   emrecord = 0;
00057   mcrecord = 0;
00058   threcord = 0;
00059   Clear();
00060   isMC = true;
00061   isTH = true;
00062   isEM = true;
00063   Nentries = entries;
00064   whichSource = 1;
00065   jcPath = path;
00066   fJC = j;
00067   fLikeliFile = NULL;
00068   for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00069 }

MadTestAnalysis::~MadTestAnalysis (  ) 

Definition at line 72 of file MadTestAnalysis.cxx.

References fLikeliFile.

00073 {
00074   if(fLikeliFile) fLikeliFile->Close();
00075 }


Member Function Documentation

AnnInputBlock MadTestAnalysis::AnnVar ( Int_t  event = 0  )  [protected]

Definition at line 118 of file MadTestAnalysis.cxx.

References AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhcommon, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwdig, AnnInputBlock::aShwph, AnnInputBlock::aShwphper, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwphperpl, AnnInputBlock::aShwphperstp, AnnInputBlock::aShwpl, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aShwstp, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrklen, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkph, AnnInputBlock::aTrkphper, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, AnnInputBlock::aTrkstp, NtpSRPlane::beg, NtpSRPlane::end, NtpSRTrack::fit, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadStrip(), MadBase::LoadTrack(), NtpSRPlane::n, NtpSRShower::ndigit, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, MadBase::ntpEvent, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRFitTrack::pass, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRShower::plane, NtpSRStrip::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRShower::stp, NtpSRTrack::stp, and NtpSRStrip::strip.

Referenced by CreatePAN().

00119 {
00120 AnnInputBlock anninput;
00121   // Initialize
00122 anninput.aTotrk       =0.;
00123 anninput.aTotstp      =0.;
00124 anninput.aTotph       =0.;
00125 anninput.aTotlen      =0.;
00126 anninput.aPhperpl     =0.;
00127 anninput.aPhperstp    =0.;
00128 
00129 
00130 anninput.aTrkpass     =0.;
00131 anninput.aTrkph       =0.;
00132 anninput.aTrklen      =0.;
00133 anninput.aTrkphperpl  =0.;
00134 anninput.aTrkphper    =0.;
00135 anninput.aTrkplu      =0.;
00136 anninput.aTrkplv      =0.;
00137 anninput.aTrkstp      =0.;
00138 
00139 anninput.aShwph       =0.;
00140 anninput.aShwstp      =0.;
00141 anninput.aShwdig      =0.;
00142 anninput.aShwpl       =0.;
00143 anninput.aShwphper    =0.;
00144 anninput.aShwphperpl  =0.;
00145 anninput.aShwphperdig =0.;
00146 anninput.aShwphperstp =0.;
00147 anninput.aShwplu      =0.;
00148 anninput.aShwplv      =0.;
00149 anninput.aPh3         =0.;
00150 anninput.aPh6         =0.;
00151 anninput.aPhlast      =0.;
00152 anninput.aPhcommon    =0.;
00153 //  
00154   
00155 // Calculate variables needed for ANN (and a few more)
00156   
00157   LoadEvent(event) ;  
00158   int track_index   = -1;
00159   LoadLargestTrackFromEvent(event,track_index);
00160   int shower_index  = -1;
00161   LoadLargestShowerFromEvent(event,shower_index);
00162   LoadTrack(track_index);
00163   LoadShower(shower_index);
00164  
00165 anninput.aTotrk       =ntpEvent->ntrack;
00166 anninput.aTotstp      =ntpEvent->nstrip;
00167 anninput.aTotph       =ntpEvent->ph.sigcor;
00168 anninput.aTotlen      =ntpEvent->plane.end-ntpEvent->plane.beg+1; // prepei na to alla3w
00169 anninput.aPhperpl     =ntpEvent->ph.sigcor/ntpEvent->plane.n;
00170 anninput.aPhperstp    =ntpEvent->ph.sigcor/ntpEvent->nstrip;
00171 
00172 if(track_index!=-1) {
00173 anninput.aTrkpass     =ntpTrack->fit.pass;
00174 anninput.aTrkph       =ntpTrack->ph.sigcor;
00175 anninput.aTrklen      =ntpTrack->plane.ntrklike;
00176 if(ntpTrack->plane.ntrklike>0)anninput.aTrkphperpl  =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
00177 if(ntpEvent->ph.sigcor>0)     anninput.aTrkphper    =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
00178 anninput.aTrkplu      =ntpTrack->plane.nu;
00179 anninput.aTrkplv      =ntpTrack->plane.nv;
00180 anninput.aTrkstp      =ntpTrack->nstrip;
00181 }
00182 if(shower_index!=-1) {
00183 anninput.aShwph       =ntpShower->ph.sigcor;
00184 anninput.aShwstp      =ntpShower->nstrip;
00185 anninput.aShwdig      =ntpShower->ndigit;
00186 anninput.aShwpl       =ntpShower->plane.n;
00187 anninput.aShwphper    =ntpShower->ph.sigcor/ntpEvent->ph.sigcor;
00188 anninput.aShwphperpl  =ntpShower->ph.sigcor/ntpShower->plane.n;
00189 anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit;
00190 anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip;
00191 anninput.aShwplu      =ntpShower->plane.nu;
00192 anninput.aShwplv      =ntpShower->plane.nv; 
00193 }
00194 
00195       Int_t ssind,ssplind;
00196       Double_t ssphtot;
00197       Bool_t foundshw,foundtrk;
00198       Int_t planes=ntpEvent->plane.beg;
00199      
00200       Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
00201       ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
00202 // loop over strips to compute ph3 ph6 phlast phcommon
00203      for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
00204        Int_t stp_index=((ntpEvent->stp)[evss]);
00205        LoadStrip(stp_index);    
00206        ssind=ntpStrip->strip;
00207        ssplind=ntpStrip->plane;
00208        ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;       
00209          
00210        foundshw=false;
00211        foundtrk=false;
00212    
00213        Double_t phstrips=0;
00214        Double_t phstript=0;       
00215        // shower strips
00216        Int_t sshwind,sshwplind;
00217        Double_t sshwphtot;
00218        
00219        if(shower_index!=-1) {
00220        for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00221         Int_t stp_indexs=((ntpShower->stp)[jj]);
00222         LoadStrip(stp_indexs);  
00223         
00224         if(!foundshw){
00225           sshwind=ntpStrip->strip;
00226           sshwplind=ntpStrip->plane;
00227           sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;             
00228             if(sshwind==ssind && sshwplind==ssplind) {
00229              foundshw=true;
00230              phstrips=sshwphtot;
00231            }
00232          }
00233        }
00234       }      
00235        // tracks strips
00236        Int_t strkind,strkplind;
00237        Double_t strkphtot;
00238       if(track_index!=-1) {
00239        for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
00240         Int_t stp_indext=((ntpTrack->stp)[jj]);
00241         LoadStrip(stp_indext);  
00242          if(!foundtrk){
00243           strkind=ntpStrip->strip;
00244           strkplind=ntpStrip->plane;
00245           strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;   
00246            if(strkind==ssind && strkplind==ssplind) {
00247             foundtrk=true;
00248             phstript=strkphtot;
00249            }
00250          }
00251        }
00252       }
00253        
00254        if(foundshw && foundtrk) {
00255         hitcommon=hitcommon+1;
00256         phcommon=phcommon+phstrips+phstript;
00257        }  
00258        if(!foundshw && ! foundtrk) {
00259          hitnowere=hitnowere+1; 
00260          phnowere=phnowere+ssphtot; 
00261         }                
00262        if(ssplind>=planes && ssplind<=(planes+3)){
00263         ph3=ph3+ssphtot;
00264        }
00265        else if(ssplind>(planes+3) && ssplind<=(planes+6)){
00266         ph6=ph6+ssphtot; 
00267        }
00268        else {
00269         phlast=phlast+ssphtot;
00270        }     
00271      } // end of looping over event strips....    
00272 
00273 //
00274 anninput.aPh3       =ph3;
00275 anninput.aPh6       =ph6;
00276 anninput.aPhlast    =phlast;
00277 anninput.aPhcommon  =phcommon; 
00278 return anninput;
00279 }

void MadTestAnalysis::CreatePAN ( std::string  tag,
Int_t  scanfilter 
)

Definition at line 724 of file MadTestAnalysis.cxx.

References AnnVar(), NtpSRPlane::beg, MadDpID::CalcPID(), NtpSREventSummary::date, det, dpid, NtpSRTrack::ds, NtpSRPlane::end, NtpSRTrack::end, NtpSRMomentum::eqp, MadBase::eventSummary, BeamMonSpill::fHornCur, MadQuantities::Flavour(), BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTrtgtd, BDSpillAccessor::Get(), GetAnnPid(), ScanList::GetDecision(), VldContext::GetDetector(), MadBase::GetEntry(), SpillInfoBlock::GetHorn(), SpillInfoBlock::GetHpos(), SpillInfoBlock::GetMagnet(), VldTimeStamp::GetNanoSec(), SpillTimeFinder::GetNearestSpill(), GnumiInterface::GetParent(), MadNsID::GetPID(), SpillInfoBlock::GetPot(), RecDataHeader::GetRun(), VldTimeStamp::GetSec(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), SpillInfo::GetSpillInfo(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), SpillInfoBlock::GetTgtpos(), SpillTimeFinder::GetTimeOfNearestSpill(), SpillTimeND::GetTimeStamp(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), SpillInfoBlock::GetVpos(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSREvent::index, MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), MadQuantities::IsFidVtxEvt(), MadBase::isST, Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), BDSpillAccessor::LoadSpill(), MadBase::LoadStrip(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), MadBase::mcrecord, NtpSRTrack::momentum, month, NtpSRDate::month, NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, nsid, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTrack, NtpSREvent::ntrack, MadQuantities::Nucleus(), PassAnalysisCuts(), PassBasicCuts(), PassFid(), NtpSREventSummary::ph, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, PID(), NtpSRTrack::plane, NtpSREvent::plane, MadQuantities::Q2(), NtpSRMomentum::qp, NtpSRMomentum::range, NtpSRPulseHeight::raw, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), MadQuantities::RecoShwEnergySqrt(), run(), NtpSRPulseHeight::sigcor, SingleTimeFrame(), BeamMonSpill::SpillTime(), GnumiInterface::Status(), NtpSREvent::stp, MadBase::strecord, MadQuantities::Target4Vector(), NtpSRStrip::time0, NtpSRStrip::time1, NtpSREventSummary::trigtime, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSRVertex::u, NtpSRVertex::v, NtpSREvent::vtx, NtpSRTrack::vtx, MadQuantities::W2(), MadQuantities::X(), NtpSRVertex::x, MadQuantities::Y(), NtpSRVertex::y, NtpSRDate::year, Munits::year, NtpSRVertex::z, NuParent::Zero(), and SpillInfoBlock::Zero().

00724                                                               {
00725 
00726   std::string savename = "PAN_" + tag + ".root";
00727   TFile save(savename.c_str(),"RECREATE"); 
00728   save.cd();
00729   
00730   GnumiInterface gnumi;
00731   if(!gnumi.Status()) {
00732     cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00733          << " Will not be filling NuParent info"
00734          << endl;
00735     cout << "Set environmental variable $GNUMIAUX to point to the "
00736          << "directory containing the gnumi files to fix this!"
00737          << endl;
00738 
00739   }
00740 
00741  SpillInfo pppinfo;
00742   
00743   //PAN Quantities 
00744   //Truth:
00745   Float_t true_enu;       // true neutrino energy (GeV)
00746   Float_t true_pxnu;      // true neutrino momentum-x (GeV)
00747   Float_t true_pynu;      // true neutrino momentum-y (GeV)
00748   Float_t true_pznu;      // true neutrino momentum-z (GeV)
00749   Float_t true_etgt;       // true target energy (GeV)
00750   Float_t true_pxtgt;      // true target momentum-x (GeV)
00751   Float_t true_pytgt;      // true target momentum-y (GeV)
00752   Float_t true_pztgt;      // true target momentum-z (GeV)
00753   Float_t true_emu;       // true muon energy
00754   Float_t true_eshw;      // true shower energy
00755   Float_t true_x;         // true x
00756   Float_t true_y;         // true y
00757   Float_t true_q2;        // true q2
00758   Float_t true_w2;        // true w2
00759   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
00760   Float_t true_dircosz;   // track z-direction cosine
00761   Float_t true_vtxx;      // true x vtx
00762   Float_t true_vtxy;      // true y vtx
00763   Float_t true_vtxz;      // true z vtx
00764 
00765   //For Neugen:
00766   Int_t flavour;          // true flavour: 1-e 2-mu 3-tau
00767   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
00768   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
00769   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
00770   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
00771   Int_t had_fs;           // hadronic final state, number between 200-300
00772   
00773   //For Beam Reweighting:
00774   NuParent *nuparent = new NuParent();
00775 
00776   //Reco Quantities
00777   Int_t detector;         // Near = 1, Far = 2;
00778   Int_t run;              // run number
00779   Int_t snarl;            // snarl number
00780   Int_t nevent;           // number of events in snarl
00781   Int_t event;            // event index
00782   Int_t subrun;           // subrun number     
00783   Int_t trigsrc;          // trigger source    
00784   Int_t mcevent;          // mc event index
00785   Int_t ntrack;           // number of tracks in event
00786   Int_t nshower;          // number of showers in event
00787 
00788   Int_t is_fid;           // pass fid cut. true = 1, false = 0
00789   Int_t is_cev;           // fully contained. true = 1, false = 0 
00790   Int_t is_mc;            // is a corresponding mc event found
00791   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
00792   Float_t pid0;           // pid parameter (usual method)
00793   Float_t pid1;           // pid parameter (alternative method 1)
00794   Float_t pid2;           // pid parameter (alternative method 2)
00795   Float_t annpid;        // ANN PID
00796   Int_t pass;             // pass analysis cuts. true = 1, false = 0
00797 
00798 
00799   Float_t npid0;           // pid parameter (usual method)
00800   Float_t npid1;           // pid parameter (alternative method 1)
00801   Float_t npid2;           // pid parameter (alternative method 2)
00802   Double_t nannpid;        // ANN PID
00803 
00804 
00805   Float_t reco_enu;       // reco neutrino energy
00806   Float_t reco_emu;       // reco muon energy
00807   Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
00808   Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
00809   Float_t reco_qe_enu;    // reco qe neutrino energy
00810   Float_t reco_qe_q2;     // reco qe q2
00811   Float_t reco_y;         // reco y
00812   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
00813   Float_t reco_dircosz;   // reco track vtx z-dircos
00814   Float_t reco_vtxx;      // reco x vtx
00815   Float_t reco_vtxy;      // reco y vtx
00816   Float_t reco_vtxz;      // reco z vtx
00817 
00818   Float_t evtlength;      // reco event length
00819   Float_t evtph;          // reco event ph
00820   Float_t trklength;      // reco track length
00821   Float_t trkmomrange;    // reco track momentum from range
00822   Float_t trkqp;          // reco track q/p
00823   Float_t trkeqp;         // reco track q/p error
00824   Float_t trkphfrac;      // reco pulse height fraction in track
00825   Float_t trkphplane;     // reco average track pulse height/plane
00826   Float_t trkds;          // track DS 
00827   Float_t rawph;          // Raw ph
00828 // Additional track info
00829   Float_t trkendz;  // track end Z position
00830   Float_t trkendx;  // track end X position
00831   Float_t trkendy;  // track end Y position
00832   Float_t trkendu;  // track end U position
00833   Float_t trkendv;  // track end V position
00834   Float_t trkvtxz;  // track begin Z position
00835   Float_t trkvtxx;  // track begin X position
00836   Float_t trkvtxy;  // track begin Y position
00837   Float_t trkvtxu;  // track begin U position
00838   Float_t trkvtxv;  // track begin V position
00839 
00840   // SCAN DECISION   // <- define an instance of ScanList class here
00841   Int_t scandecision;
00842   ScanList scan;
00843   
00844 // EVENT TIMING
00845    Double_t evttimemin;  //MIN STRIP TIME OF EVENT
00846    Double_t evttimemax;  //MAX STRIP TIME OF EVENT
00847    Double_t snarlt;
00848 
00849 // WEIGHT FOR THE PME BEAM
00850    Double_t w_le_me;  
00851 // BEAM INFO
00852    Float_t  pot,horn,tar,vvpos,hhpos,magn;    
00853 // Beam INFO DATABASE
00854    Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
00855    Double_t timedb;
00856 // TIME INFO DATABASE
00857    Double_t neartdb,fartdb,neardifftdb;
00858 
00859 // IS SNARL CUT IN PIECES 
00860    Float_t singletimeframe;
00861 // Beam Info Block    
00862    SpillInfoBlock *spinfo = new SpillInfoBlock();  
00863    spinfo->Zero(); 
00864    singletimeframe=-999;
00865     
00866   //Trees
00867   TTree *tree = new TTree("pan","pan");
00868   //Truth
00869   tree->Branch("true_enu",&true_enu,"true_enu/F",512000);
00870   tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000);
00871   tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000);
00872   tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000);
00873   tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000);
00874   tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000);
00875   tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000);
00876   tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000);
00877   tree->Branch("true_emu",&true_emu,"true_emu/F",512000);
00878   tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000);
00879   tree->Branch("true_x",&true_x,"true_x/F",512000);
00880   tree->Branch("true_y",&true_y,"true_y/F",512000);
00881   tree->Branch("true_q2",&true_q2,"true_q2/F",512000);
00882   tree->Branch("true_w2",&true_w2,"true_w2/F",512000);
00883   tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000);
00884   tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000);
00885   tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000);
00886   tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000);
00887   tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000);
00888   //Neugen
00889   tree->Branch("flavour",&flavour,"flavour/I",512000);
00890   tree->Branch("process",&process,"process/I",512000);
00891   tree->Branch("nucleus",&nucleus,"nucleus/I",512000);
00892   tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000);
00893   tree->Branch("initial_state",&initial_state,"initial_state/I",512000);
00894   tree->Branch("had_fs",&had_fs,"had_fs/I",512000);
00895   //NuParent
00896   tree->Branch("nuparent","NuParent",&nuparent,8000,1);
00897   //Reco
00898   tree->Branch("detector",&detector,"detector/I",512000);
00899   tree->Branch("run",&run,"run/I",512000);
00900   tree->Branch("snarl",&snarl,"snarl/I",512000);
00901   tree->Branch("event",&event,"event/I",512000);
00902   tree->Branch("mcevent",&mcevent,"mcevent/I",512000);
00903   tree->Branch("ntrack",&ntrack,"ntrack/I",512000);
00904   tree->Branch("nshower",&nshower,"nshower/I",512000);
00905   tree->Branch("subrun",&subrun,"subrun/I",512000);   
00906   tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000); 
00907   
00908   tree->Branch("is_fid",&is_fid,"is_fid/I",512000);
00909   tree->Branch("is_cev",&is_cev,"is_cev/I",512000);
00910   tree->Branch("is_mc",&is_mc,"is_mc/I",512000);
00911   tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000);
00912   tree->Branch("pid0",&pid0,"pid0/F",512000);
00913   tree->Branch("pid1",&pid1,"pid1/F",512000);
00914   tree->Branch("pid2",&pid2,"pid2/F",512000);
00915   tree->Branch("annpid",&annpid,"annpid/F",512000);  
00916   tree->Branch("pass",&pass,"pass/I",512000);
00917 
00918   tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000);
00919   tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000);
00920   tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000);
00921   tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000);
00922   tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000);
00923   tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000);
00924   tree->Branch("reco_y",&reco_y,"reco_y/F",512000);
00925   tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000);
00926   tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000);
00927   tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000);
00928   tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000);
00929   tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000);
00930 
00931   tree->Branch("evtlength",&evtlength,"evtlength/F",512000);
00932   tree->Branch("evtph",&evtph,"evtph/F",512000);
00933   tree->Branch("trklength",&trklength,"trklength/F",512000);
00934   tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000);
00935   tree->Branch("trkqp",&trkqp,"trkqp/F",512000);
00936   tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000);
00937   tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000);
00938   tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000);
00939   tree->Branch("rawph",&rawph,"rawph/F",512000);    // <- !! added this
00940   
00941   tree->Branch("w_le_me",&w_le_me,"w_le_me/D",512000); 
00942   tree->Branch("trkendz",&trkendz,"trkendz/F",512000);
00943   tree->Branch("trkendu",&trkendu,"trkendu/F",512000);
00944   tree->Branch("trkendv",&trkendv,"trkendv/F",512000);
00945   tree->Branch("trkendx",&trkendx,"trkendx/F",512000);
00946   tree->Branch("trkendy",&trkendy,"trkendy/F",512000);
00947 
00948   tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000);
00949   tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000);
00950   tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000);
00951   tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000);
00952   tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000);
00953   tree->Branch("trkds",  &trkds,"trkds/F",512000);
00954   
00955   tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
00956   tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
00957   
00958   // SCAN DECISIONS 
00959   tree->Branch("scandecision",&scandecision,"scandecision/I");
00960 
00961   //BEAM INFO 
00962   tree->Branch("pot",    &pot,    "pot/F",512000);
00963   tree->Branch("horn",   &horn,   "horn/F",512000);
00964   tree->Branch("tar",    &tar,    "tar/F",512000);
00965   tree->Branch("vvpos",  &vvpos,  "vvpos/F",512000);
00966   tree->Branch("hhpos",  &hhpos,  "hhpos/F",512000);
00967   tree->Branch("magn",   &magn,   "magn/F",512000);
00968   tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000);
00969 
00970   // BEAM INFO DB
00971   tree->Branch("potdb",     &potdb,      "potdb/F",512000);
00972   tree->Branch("horndb",    &horndb,     "horndb/F",512000);
00973   tree->Branch("tardb",     &tardb,      "tardb/F",512000);
00974   tree->Branch("vvposdb",   &vvposdb,    "vvposdb/F",512000);
00975   tree->Branch("hhposdb",   &hhposdb,    "hhposdb/F",512000);
00976   tree->Branch("vvwidthdb", &vvwidthdb,  "vvwidthdb/F",512000);
00977   tree->Branch("hhwidthdb", &hhwidthdb,  "hhwidthdb/F",512000);
00978   tree->Branch("timedb",    &timedb,     "timedb/D");
00979 
00980   // TIME DB  
00981   tree->Branch("neartdb",      &neartdb,    "neartdb/D");
00982   tree->Branch("fartdb",       &fartdb,     "fartdb/D");
00983   tree->Branch("neardifftdb",  &neardifftdb,"neardifftdb/D"); 
00984   tree->Branch("snarlt",       &snarlt,     "snarlt/D");
00985   
00986   // new pid variables
00987   tree->Branch("npid0",&npid0,"npid0/F",512000);
00988   tree->Branch("npid1",&npid1,"npid1/F",512000);
00989   tree->Branch("npid2",&npid2,"npid2/F",512000);
00990   tree->Branch("nannpid",&nannpid,"nannpid/D",512000);  
00991   
00992 
00993   
00994 
00995   // MAIN LOOP (over tree entries)
00996 
00997   for(int i=0;i<Nentries;i++){
00998  
00999     if(i%100==0) std::cout << float(i)*100./float(Nentries) 
01000                             << "% done" << std::endl;
01001 
01002     if(!this->GetEntry(i)) continue;
01003 
01004     nevent = eventSummary->nevent;
01005     if(nevent==0) continue;
01006     
01007     run     = ntpHeader->GetRun();
01008     snarl   = ntpHeader->GetSnarl();
01009     subrun  = ntpHeader->GetSubRun();   
01010     trigsrc = ntpHeader->GetTrigSrc(); 
01011     
01012     Double_t snarltime=ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);       
01013     Int_t    month    =eventSummary->date.month; 
01014     Int_t    year     =eventSummary->date.year;
01015 
01016     // SCAN INFO  
01017     
01018     scandecision=scan.GetDecision(run,snarl);
01019 
01020     // only pass through scanned snarls if scanfilter=1
01021     if (scanfilter==0 || scandecision>0) {  
01022 
01023       // Get Beam Monitorint Info if NOT MC: 
01024       
01025       if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01026         pppinfo.GetSpillInfo(month,year,snarltime,*spinfo); 
01027         pot   =spinfo->GetPot();   
01028         horn  =spinfo->GetHorn();
01029         hhpos =spinfo->GetHpos();
01030         vvpos =spinfo->GetVpos();
01031         tar   =spinfo->GetTgtpos();
01032         magn  =spinfo->GetMagnet();   
01033         
01034         potdb       = -999;
01035         horndb      = -999; 
01036         tardb       = -999; 
01037         vvposdb     = -999; 
01038         hhposdb     = -999; 
01039         vvwidthdb   = -999; 
01040         hhwidthdb   = -999; 
01041         timedb      = -999; 
01042         neartdb     = -999; 
01043         fartdb      = -999; 
01044         neardifftdb = -999; 
01045         
01046         // BEAM INFO DB
01047         Detector::Detector_t detectort;
01048         SimFlag::SimFlag_t simflag;
01049         VldTimeStamp  timestamp;
01050         VldContext    cx;
01051         cx        = ntpHeader->GetVldContext(); 
01052         detectort = ntpHeader->GetVldContext().GetDetector();
01053         simflag   = ntpHeader->GetVldContext().GetSimFlag();
01054         timestamp = ntpHeader->GetVldContext().GetTimeStamp();
01055         
01056         BDSpillAccessor &spillAccess=BDSpillAccessor::Get();
01057         const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp);
01058         if(spillInfoDB) { 
01059           if(spillInfoDB->fTortgt !=0.)      potdb     = spillInfoDB->fTortgt;
01060           else if(spillInfoDB->fTrtgtd !=0.) potdb     = spillInfoDB->fTrtgtd;
01061           else if(spillInfoDB->fTor101 !=0.) potdb     = spillInfoDB->fTor101;
01062           else potdb = -999;
01063           
01064           horndb    = spillInfoDB->fHornCur;
01065           tardb     = (int)spillInfoDB->GetStatusBits().target_in;
01066           vvposdb   = spillInfoDB->fTargBpmY[0];
01067           hhposdb   = spillInfoDB->fTargBpmX[0];
01068           vvwidthdb = spillInfoDB->fProfWidY;
01069           hhwidthdb = spillInfoDB->fProfWidX;
01070           timedb    = spillInfoDB->SpillTime();
01071         }
01072         
01073         // TIME DB
01074         SpillTimeFinder& spillfinder = SpillTimeFinder::Instance();
01075         fartdb                       = spillfinder.GetTimeOfNearestSpill(cx);
01076         const SpillTimeND& spnd      = spillfinder.GetNearestSpill(cx);
01077         neartdb                      = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9);
01078         neardifftdb                  = spillfinder.GetTimeToNearestSpill(cx);
01079         snarlt                       = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
01080         
01081       }
01082       singletimeframe=SingleTimeFrame(i,Nentries);
01083       
01084       Int_t trkcount=0;
01085       Int_t shwcount=0;
01086       Float_t totph=0;
01087       
01088       // EVENT LOOP - fill tree once for each reconstructed event
01089 
01090       for(int j=0;j<nevent;j++){ 
01091 
01092         if(!LoadEvent(j)) continue; //no event found
01093 
01094         totph = 0;
01095 
01096         //set event index:
01097         event = j;
01098         ntrack = ntpEvent->ntrack;
01099         nshower = ntpEvent->nshower;
01100         
01101         // Initialize everything
01102         nuparent->Zero();     
01103         true_enu = 0; true_emu = 0; true_eshw = 0; 
01104         true_pxnu = 0; true_pynu = 0; true_pznu = 0;
01105         true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
01106         true_dircosneu = 0; true_dircosz = 0;
01107         true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
01108         flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
01109         initial_state = 0; had_fs = 0;
01110         true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;     
01111         detector = -1; mcevent = -1;
01112         is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; 
01113         pid0 = -999; pid1 = -999; pid2 = -999; pass = 0;
01114         reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0; 
01115         reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0; 
01116         reco_dircosz = 0;
01117         reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999;
01118         evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1;
01119         trkmomrange = -999; trkqp = -999; trkeqp = -999;
01120         trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.;   
01121         w_le_me=-1; trkendz=999; trkendu=999; trkendv=999; trkendx=999;
01122         trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999;
01123         trkds=-1; evttimemin=-1; evttimemax=-1;
01124         evtph=0.;
01125         annpid=-999.;
01126         npid0 = -999; npid1 = -999; npid2 = -999;
01127         nannpid=-999;
01128 
01129         rawph=ntpEvent->ph.raw;
01130  
01131         //get detector type:
01132         if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01133           detector = 2;
01134           totph=eventSummary->ph.sigcor;
01135           //fiducial criteria on vtx for far detector
01136           is_fid = 1;
01137           if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01138         }
01139         else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01140           detector = 1;
01141         
01142           //get total charge in trk/shw
01143           for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
01144             LoadTrack(ii);
01145             totph+=ntpTrack->ph.sigcor;
01146           }
01147           trkcount+=ntrack;
01148           
01149           for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
01150             LoadShower(ii);
01151             totph+=ntpShower->ph.sigcor;
01152           }
01153           shwcount+=nshower;
01154           
01155           //fiducial criteria on vtx for near detector
01156           is_fid = 1;
01157           if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01158         }
01159         
01160         //check for a corresponding mc event      
01161         if(LoadTruthForRecoTH(j,mcevent)) {
01162           Float_t *vtx = TrueVtx(mcevent);
01163           Float_t *nu3mom = TrueNuMom(mcevent);
01164           Float_t *targ4vec = Target4Vector(mcevent);
01165           //true info for tree:
01166           is_mc             = 1;
01167           nucleus           = Nucleus(mcevent);
01168           flavour           = Flavour(mcevent);
01169           initial_state     = Initial_state(mcevent);
01170           had_fs            = HadronicFinalState(mcevent);
01171           process           = IResonance(mcevent); 
01172           cc_nc             = IAction(mcevent);
01173           true_enu          = TrueNuEnergy(mcevent); 
01174           true_pxnu         = nu3mom[0];
01175           true_pynu         = nu3mom[1];
01176           true_pznu         = nu3mom[2];
01177           true_emu          = TrueMuEnergy(mcevent); 
01178           true_eshw         = TrueShwEnergy(mcevent); 
01179           true_x            = X(mcevent);
01180           true_y            = Y(mcevent);
01181           true_q2           = Q2(mcevent);
01182           true_w2           = W2(mcevent);
01183           true_dircosz      = TrueMuDCosZ(mcevent);
01184           true_dircosneu    = TrueMuDCosNeu(mcevent);
01185           true_vtxx         = vtx[0];
01186           true_vtxy         = vtx[1];
01187           true_vtxz         = vtx[2];
01188           true_etgt         = targ4vec[3];
01189           true_pxtgt        = targ4vec[0];
01190           true_pytgt        = targ4vec[1];
01191           true_pztgt        = targ4vec[2];
01192           
01193           if(gnumi.Status()) {
01194             if(isST) gnumi.GetParent(strecord,mcevent,*nuparent);
01195             else     gnumi.GetParent(mcrecord,mcevent,*nuparent);
01196           }
01197           
01198           delete [] vtx;
01199           delete [] nu3mom;
01200           delete [] targ4vec;
01201         }
01202         
01203         
01204         if(PassBasicCuts(event)) {pass_basic=1;} // make this reflect pass/fail flag
01205         
01206         //reco info for tree:
01207         reco_vtxx         = ntpEvent->vtx.x;
01208         reco_vtxy         = ntpEvent->vtx.y;
01209         reco_vtxz         = ntpEvent->vtx.z;
01210         evtlength         = ntpEvent->plane.end-ntpEvent->plane.beg+1;
01211         evtph             = ntpEvent->ph.sigcor; 
01212         
01213         // different definition for neardet - accounts for 
01214         // uninstrumented planes
01215    
01216         // ANN PID // 
01217         
01218         AnnInputBlock anninput=AnnVar(event);
01219         //      AnnInputBlock anninput2=anninput;
01220         //      anninput2.aTotrk=1000;
01221         //      if(!nsid.CompareAnnBlocks(this,event,anninput2)) assert(false);
01222 
01223         Int_t target=0;
01224         if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01225           if(tar<900) target=1;
01226           if(tar>900 && tar<40000.) target=2;
01227           if(tar>40000 && tar <100000) target=3;
01228         }
01229         if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){
01230           Int_t mod=(int)(run/1e6);
01231           if(mod==14) target=1; 
01232           if(mod==16) target=2;
01233           if(mod==18) target=3; 
01234         }
01235         if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){
01236           target=0; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC 
01237         }                       
01238         // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING
01239         
01240         Int_t fid    =2;
01241         Int_t prior = 1;
01242         
01243         if(!PassFid(event)) { 
01244           annpid=-999;
01245           nannpid=-999;
01246         }
01247         else {
01248           Detector::Detector_t det=
01249             ntpHeader->GetVldContext().GetDetector();
01250           // ANN PID //           
01251           AnnInputBlock anninput=AnnVar(event);  
01252           if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01253             annpid=GetAnnPid(anninput,detector,target,fid,prior);
01254             if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999;
01255             
01256             if(!((reco_vtxz>1. &&reco_vtxz<5.)&&sqrt((reco_vtxx-1.4885)*(reco_vtxx-1.4885)+(reco_vtxy-0.1397)*(reco_vtxy-0.1397))<1)) {
01257               annpid=-999.; nannpid=-999;
01258             }
01259           } 
01260           
01261           if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01262             annpid=GetAnnPid(anninput,detector,target,fid,prior);  
01263             if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999;
01264             
01265             if(fid==1) if(!(((reco_vtxz>1. &&reco_vtxz<14.)||(reco_vtxz>17. &&reco_vtxz<27.))&& sqrt(reco_vtxx*reco_vtxx+reco_vtxy*reco_vtxy)<3.5)) {
01266               annpid=-999.; nannpid=-999;
01267             }
01268           }
01269         }
01270         
01271         Double_t timemin=9.*10e30; 
01272         Double_t timemax=-9999999; 
01273         Double_t trgtime=eventSummary->trigtime;
01274 
01275         //  Find max min time of events by loading strips for ND 
01276         if(detector==1){
01277           for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
01278             Int_t stp_index=((ntpEvent->stp)[evss]);
01279             LoadStrip(stp_index); 
01280             Double_t striptime=ntpStrip->time1;
01281             if(striptime<=timemin) timemin=striptime;
01282             if(striptime>=timemax) timemax=striptime;
01283           }
01284           evttimemax=timemax-trgtime;
01285           evttimemin=timemin-trgtime;  
01286         }        
01287         
01288         if(detector==2){  
01289           for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
01290             Int_t stp_index=((ntpEvent->stp)[evss]);
01291             LoadStrip(stp_index);
01292             Double_t striptime1=ntpStrip->time1;
01293             Double_t striptime0=ntpStrip->time0;
01294             Double_t striptime = 0;
01295             if(striptime1>0 && striptime0<0)  striptime=striptime1;
01296             if(striptime0>0 && striptime1<0)  striptime=striptime0;
01297             if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
01298             if(striptime<=timemin) timemin=striptime;
01299             if(striptime>=timemax) timemax=striptime;
01300           }
01301           evttimemax=timemax-trgtime;
01302           evttimemin=timemin-trgtime;
01303         }
01304   
01305         //index will -1 if no track/shower present in event
01306         int track_index   = -1;
01307         LoadLargestTrackFromEvent(j,track_index);
01308         int shower_index  = -1;
01309         LoadLargestShowerFromEvent(j,shower_index);
01310         
01311         //methods should all be protected against index = -1
01312         reco_emu          = RecoMuEnergy(0,track_index);
01313         reco_eshw         = RecoShwEnergy(shower_index);
01314         reco_eshw_sqrt    = RecoShwEnergySqrt(event);
01315         reco_enu          = reco_emu + reco_eshw;
01316         reco_qe_enu       = RecoQENuEnergy(track_index);
01317         if (reco_enu>0) {
01318           reco_y          = reco_eshw/reco_enu;
01319         }
01320         reco_qe_q2        = RecoQEQ2(track_index);
01321         reco_dircosz      = RecoMuDCosZ(track_index);
01322         reco_dircosneu    = RecoMuDCosNeu(track_index);
01323         is_cev            = IsFidAll(track_index);
01324         
01325         if(track_index!=-1){ //check track is present before using ntpTrack
01326           trklength         = ntpTrack->plane.end-ntpTrack->plane.beg+1; 
01327           trkmomrange       = ntpTrack->momentum.range;
01328           trkqp             = ntpTrack->momentum.qp;
01329           trkeqp            = ntpTrack->momentum.eqp;
01330           trkendz           = ntpTrack->end.z;
01331           trkendu           = ntpTrack->end.u;
01332           trkendv           = ntpTrack->end.v;
01333           trkendx           = ntpTrack->end.x;
01334           trkendy           = ntpTrack->end.y;
01335           
01336           trkvtxz           = ntpTrack->vtx.z;
01337           trkvtxu           = ntpTrack->vtx.u;
01338           trkvtxv           = ntpTrack->vtx.v;
01339           trkvtxx           = ntpTrack->vtx.x;
01340           trkvtxy           = ntpTrack->vtx.y;
01341           trkds             = ntpTrack->ds;               
01342           Float_t phtrack=ntpTrack->ph.sigcor;
01343           Float_t phevent=ntpEvent->ph.sigcor;
01344           if (phevent>0) {trkphfrac=phtrack/phevent;}
01345           if(ntpTrack->plane.n>0) {
01346             trkphplane      = ntpTrack->ph.sigcor/ntpTrack->plane.n;
01347           }
01348         }
01349         else {
01350           trklength         = 0;
01351           trkmomrange       = 0;
01352           trkqp             = 0;
01353           trkeqp            = 0;
01354           trkphfrac         = 0;
01355           trkphplane        = 0;
01356         }
01357 
01358         // only calculate likelihood PID for events that pass event 
01359         // fiducial cuts and track quality cuts 
01360 
01361         if(PassBasicCuts(event) && PassFid(event)){
01362           pid0              = PID(event,0);
01363           pid1              = PID(event,1);
01364           pid2              = PID(event,2);
01365           npid0 = dpid.CalcPID(this,event,0);
01366           npid1 = dpid.CalcPID(this,event,1);
01367           npid2 = dpid.CalcPID(this,event,2);
01368           
01369           if(PassAnalysisCuts(event)) pass = 1;
01370         }
01371         else{ 
01372           pid0            = -999.;
01373           pid1            = -999.;
01374           pid2            = -999.;
01375           pass              =  0;
01376           npid0           = -999.;
01377           npid1           = -999.;
01378           npid2           = -999.;
01379 
01380         }
01381         tree->Fill();
01382 
01383       } // end of EVENT LOOP
01384     } //  end of SCANINFO loop 
01385   } // end of NTUPLE ENTRIES loop
01386 
01387   delete nuparent;
01388   delete spinfo;
01389   
01390   save.cd();
01391   tree->Write();
01392   save.Write();
01393   save.Close();
01394 }

Double_t MadTestAnalysis::GetAnnPid ( AnnInputBlock  anninput,
Int_t  det,
Int_t  tar,
Int_t  fid,
Int_t  prior 
) [protected]

Definition at line 297 of file MadTestAnalysis.cxx.

References AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, and Sigmoid().

Referenced by CreatePAN().

00298 {
00299    // DET IS THE DETECTOR 1 = NEAR 2 = FAR
00300    // TAR IS THE ENERGY   1 = LE 2 = PME 3 = PHE
00301    // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS  1 = MINE (MORE STRICT)
00302    // PRIOR IS THE FAR TRAINED METHOD    , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 
00303    //                                      3 = DM2 0.002 SINTHETA2 = 0.95
00304       
00305       int    inneuron  = 13;
00306       int    hidneuron = 15;
00307       double var;
00308      
00309       double out[15]; // input neurons
00310       double rin[15]; // first hidden layer
00311      
00312       double weight1[13][15];// weights first layer
00313       double constant1[15];   // constant first layer
00314       
00315       double weighto[15];// weights second  (output) layer
00316       double constanto[1];   // constant second (output) layer  
00317       double prob;
00318       
00319 // Initialize    
00320 
00321        for(Int_t i=0;i<hidneuron;i++) {
00322          out[i]=0.;
00323          rin[i]=0.; 
00324          constant1[i]=0.;        
00325        }
00326             
00327        for(Int_t i=0;i<inneuron;i++) {
00328         for(Int_t j=0;j<hidneuron;j++){
00329           weight1[i][j]=0.;
00330         }
00331        }
00332 
00333         constanto[0]=0.;
00334         for(Int_t i=0;i<hidneuron;i++){
00335           weighto[i]=0.;
00336         }
00337       Int_t all  = inneuron*hidneuron+hidneuron+hidneuron+1;
00338       Int_t all1 = inneuron*hidneuron+hidneuron;
00339       Int_t n1  =0;
00340       Int_t nw1 =0;
00341       Int_t no  =0;
00342       Int_t nwo =0;
00343       ifstream weightfile;
00344       
00345 //   INPUT VARIABLES
00346       if(det==2){  
00347        if(prior==1){ 
00348         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnew12.dat");    
00349         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnewfid12.dat");        
00350        }
00351        if(prior==2){ 
00352         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1new12.dat");    
00353         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1newfid12.dat");
00354        }
00355        if(prior==3){ 
00356         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2new12.dat");    
00357         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2newfid12.dat");
00358        }
00359       out[0]  = anninput.aTrkphperpl/2000.;
00360       out[1]  = anninput.aTotrk*anninput.aTrkpass;
00361       out[2]  = (anninput.aTotstp/100.);
00362       out[3]  = (anninput.aTotph/50000.);   
00363       out[4]  = (anninput.aTotlen/40);
00364       out[5]  = (anninput.aPhperpl/5000.);      
00365       out[6]  = (anninput.aPhperstp/1000.);
00366       out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.;
00367       out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
00368       }
00369       else if(det==1){
00370        if(tar==1)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle16.dat");    
00371        if(tar==2)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme16.dat");    
00372        if(tar==3)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe16.dat");    
00373        out[0]  = anninput.aTrkphperpl/5000.;
00374        out[1]  = anninput.aTotrk*anninput.aTrkpass;
00375        out[2]  = (anninput.aTotstp/100.);
00376        out[3]  = (anninput.aTotph/100000.);   
00377        out[4]  = (anninput.aTotlen/40);
00378        out[5]  = (anninput.aPhperpl/10000.);    
00379        out[6]  = (anninput.aPhperstp/2000.);  
00380        out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.;
00381        out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00382       }
00383       out[7]  = (anninput.aPh3/(anninput.aTotph+1.));
00384       out[8]  = (anninput.aPh6/(anninput.aTotph+1.));
00385       out[9]  = (anninput.aPhlast/(anninput.aTotph+1.));
00386       out[12] = (anninput.aShwphperdig/800.);
00387       
00388 // INPUT FILE
00389         
00390      for(Int_t k=0;k<all;k++){  
00391        weightfile >> var ;
00392        if(k<all1){
00393          if(k%(inneuron+1)==0){
00394           nw1=0;
00395           constant1[n1]=var;      
00396           n1=n1+1;         
00397          }
00398          else {
00399            weight1[nw1][n1-1]=var;         
00400            nw1=nw1+1;      
00401          }
00402         }
00403         else if(k>=all1){
00404           if((k-all1)%(hidneuron+1)==0){
00405           nwo=0;
00406           constanto[no]=var;      
00407           no=no+1;        
00408          }
00409          else {
00410            weighto[nwo]=var;       
00411            nwo=nwo+1;      
00412          }
00413         }          
00414        }
00415        
00416      weightfile.close();  
00417 // end of reading the weights
00418    
00419        // first layer         
00420        for(Int_t i=0;i<hidneuron;i++){
00421         rin[i]=constant1[i];
00422         for(Int_t j=0;j<inneuron;j++){
00423          rin[i]=rin[i]+weight1[j][i]*out[j];
00424         }
00425        }
00426        // output                        
00427         for(Int_t i=0;i<hidneuron;i++){  
00428          out[i]=Sigmoid(rin[i]);         
00429         }               
00430         
00431         prob=constanto[0];      
00432         for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
00433         
00434                         
00435         if(anninput.aTotlen>40) prob=0.;
00436                 
00437  return prob;
00438 }

bool MadTestAnalysis::InFidVol ( const Detector::Detector_t det,
const Float_t &  x,
const Float_t &  y,
const Float_t &  z 
) [static]

Definition at line 1422 of file MadTestAnalysis.cxx.

References Detector::kFar, and Detector::kNear.

01424                                               {
01425   // Mike Kordosky: 11 August 2005
01426   // Defines a fiducial volume used to pass events which fill the PDFs
01427   // Transcribed from MakeMyFile() above
01428   
01429   bool is_fid=false;
01430   
01431   if(det==Detector::kFar){
01432     
01433     //fiducial criteria on vtx for far detector
01434     /*
01435       if (evlength_cc->GetNbinsX()==60) {
01436       evlength_cc->SetBins(50,0,500);
01437       evlength_nc->SetBins(50,0,500);
01438       }
01439     */
01440     is_fid = true;
01441     if(z<0.5 || z>29.4 ||   //ends
01442        (z<16.5&&z>14.5) ||  //between SMs
01443        sqrt((x*x)           //radial cut
01444             +(y*y))>3.5) is_fid = false;
01445   }
01446   else if(det==Detector::kNear){
01447     
01448     /*
01449       if (evlength_cc->GetBinLowEdge(50)>200) {
01450       evlength_cc->SetBins(60,0,300);
01451       evlength_nc->SetBins(60,0,300);
01452       }
01453     */
01454     //fiducial criteria on vtx for near detector
01455     is_fid = true;
01456     if(z<1 || z>5 ||
01457        sqrt(((x-1.4885)*(x-1.4885)) +
01458             ((y-0.1397)*(y-0.1397)))>1) is_fid = false;
01459   }
01460 
01461   return is_fid;
01462   
01463 }

void MadTestAnalysis::MakeMyFile ( std::string  tag,
int  tarpos 
)

Definition at line 510 of file MadTestAnalysis.cxx.

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

00510                                                         {
00511 
00512   std::string savename = "DPHistos_" + tag + ".root";
00513   TFile save(savename.c_str(),"RECREATE"); 
00514   save.cd();
00515   
00516   //#declare lots of histos:
00517 
00518   TH1F *evlength_nc = new TH1F("Event length - nc",
00519                                "Event length - nc",
00520                                60,0,600);
00521   evlength_nc->SetXTitle("Event length (planes)");
00522   TH1F *evlength_cc = new TH1F("Event length - cc",
00523                                "Event length - cc",
00524                                60,0,600);
00525   evlength_cc->SetXTitle("Event length (planes)");
00526 
00527 
00528   TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
00529                                "Track ph frac - nc",
00530                                50,0,1);
00531   phfrac_nc->SetXTitle("pulse height fraction");
00532 
00533 
00534   TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
00535                                "Track ph frac - cc",
00536                                50,0,1);
00537   phfrac_cc->SetXTitle("pulse height fraction");
00538 
00539 
00540   TH1F *phplane_nc = new TH1F("ph per plane (nc)",
00541                               "ph per plane (nc)",
00542                                50,0,5000);
00543   phplane_nc->SetXTitle("pulse height per plane");
00544   TH1F *phplane_cc = new TH1F("ph per plane (cc)",
00545                               "ph per plane (cc)",
00546                                50,0,5000);
00547   phplane_cc->SetXTitle("pulse height per plane");
00548 
00549 
00550 
00551   for(int i=0;i<Nentries;i++){
00552 
00553     
00554     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00555                             << "% done" << std::endl;
00556    if(!this->GetEntry(i)) continue;
00557 
00558     Int_t nevent = eventSummary->nevent;
00559     if(nevent==0) continue;
00560     
00561     Int_t trkcount=0;
00562     Int_t shwcount=0;
00563 
00564         
00565     //fill tree once for each reconstructed event
00566     for(int j=0;j<eventSummary->nevent;j++){
00567 
00568       if(!LoadEvent(j)) continue; //no event found
00569 
00570       Int_t ntrack = 0; 
00571       ntrack=ntpEvent->ntrack;
00572       Int_t nshower = 0;
00573       nshower=ntpEvent->nshower;
00574 
00575       Float_t totph=0;
00576 
00577 
00578       for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00579         LoadTrack(ii);
00580 
00581         totph+=ntpTrack->ph.sigcor;
00582         LoadTHTrack(ii);
00583 
00584       }
00585 
00586       trkcount+=ntrack;
00587 
00588 
00589        for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00590         LoadShower(ii);
00591 
00592         totph+=ntpShower->ph.sigcor;
00593         }
00594 
00595        shwcount+=nshower;
00596  
00597 
00598 
00599       Int_t cc_nc = 0;
00600       Int_t flavour = 0;
00601       Int_t detector = -1;
00602       Int_t is_fid = 0;
00603       Double_t neu_e  = -1; 
00604       Double_t weight = -1; 
00605       
00606       //get detector type:
00607       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00608         detector = 2;
00609         //fiducial criteria on vtx for far detector
00610         if (evlength_cc->GetNbinsX()==60) {
00611           evlength_cc->SetBins(50,0,500);
00612           evlength_nc->SetBins(50,0,500);
00613         }
00614         is_fid = 1;
00615         if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
00616            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00617            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00618                 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0;
00619       }
00620       else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00621         detector = 1;   
00622 
00623         if (evlength_cc->GetBinLowEdge(50)>200) {
00624           evlength_cc->SetBins(60,0,300);
00625           evlength_nc->SetBins(60,0,300);
00626         }
00627 
00628         //fiducial criteria on vtx for near detector
00629         is_fid = 1;
00630         if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
00631            sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
00632                 ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) is_fid = 0;
00633       }
00634       
00635       //check for a corresponding mc event
00636  
00637       Int_t mcevent=-1;
00638 
00639       if(LoadTruthForRecoTH(j,mcevent)) {
00640         //true info for tree:
00641 
00642         cc_nc             = IAction(mcevent);
00643         flavour           = Flavour(mcevent);
00644         neu_e             = TrueNuEnergy(mcevent);  
00645                           
00646    // tarpos   1 = LE 2 = ME 3 = HE     
00647           if(tarpos==2)  weight=1;
00648           if(tarpos==1)  weight=1;      
00649           if(tarpos==3)  weight=1; 
00650       }
00651 
00652       int track_index   = -1;
00653       LoadLargestTrackFromEvent(j,track_index);
00654 
00655       if (track_index!=-1) {
00656 
00657         Int_t trkpass=ntpTrack->fit.pass;
00658 
00659         if (is_fid && trkpass) {
00660 
00661         Float_t evlength=0;
00662         Float_t phfrac=0;
00663         Float_t phplane=0;
00664 
00665         evlength=ntpEvent->plane.n;
00666         if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00667           evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00668         }
00669 
00670         Float_t phtrack=ntpTrack->ph.sigcor;
00671         Float_t phevent=ntpEvent->ph.sigcor;
00672         if (phevent>0) {phfrac=phtrack/phevent;}
00673 
00674         Float_t trkplane=ntpTrack->plane.n;
00675         if (trkplane>0) {phplane=phtrack/trkplane;}
00676        
00677         if(flavour==2 && cc_nc==1) {
00678           evlength_cc->Fill(evlength,weight);
00679           phfrac_cc->Fill(phfrac,weight);
00680           phplane_cc->Fill(phplane,weight);      
00681         }
00682 
00683         else if(cc_nc==0) {
00684           evlength_nc->Fill(evlength,weight);
00685           phfrac_nc->Fill(phfrac,weight);
00686           phplane_nc->Fill(phplane,weight);      
00687         }
00688 
00689 
00690         }
00691       }
00692     }
00693   }
00694   //#close file  
00695   save.Write();
00696   save.Close();
00697 }

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

Reimplemented from MadAnalysis.

Definition at line 86 of file MadTestAnalysis.cxx.

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

Referenced by CreatePAN().

00087 {
00088   if(!LoadEvent(event)) return false;
00089   //  if(PassBasicCuts(event)) return true;
00090   Float_t pidcut=-0.4;
00091   if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
00092   
00093   if(!(PID(event,0)>pidcut)) return false; 
00094   return true;
00095 }

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

Reimplemented from MadAnalysis.

Definition at line 98 of file MadTestAnalysis.cxx.

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

Referenced by CreatePAN().

00099 { 
00100   if(!LoadEvent(event)) return false;
00101    if(ntpEvent->ntrack==0) return false;
00102    Int_t track = -1;
00103    LoadLargestTrackFromEvent(event,track);  
00104    if(!PassTrackCuts(track)) return false;  
00105    if(!IsFidVtxEvt(ntpEvent->index)) return false;  
00106 //  if(!IsFidVtx(track)) return false;
00107   return true;
00108 }

Bool_t MadTestAnalysis::PassFid ( Int_t  event = 0  )  [protected]

Definition at line 109 of file MadTestAnalysis.cxx.

References NtpSREvent::index, MadQuantities::IsFidVtxEvt(), MadBase::LoadEvent(), and MadBase::ntpEvent.

Referenced by CreatePAN().

00110 { 
00111   if(!LoadEvent(event)) return false;
00112   if(!IsFidVtxEvt(ntpEvent->index)) return false;  
00113   return true;
00114 }

Bool_t MadTestAnalysis::PassTruthSignal ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 78 of file MadTestAnalysis.cxx.

References NtpMCTruth::iaction, NtpMCTruth::inu, MadBase::LoadTruth(), and MadBase::ntpTruth.

00079 {
00080   if(!LoadTruth(mcevent)) return false;
00081   if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
00082   return false;
00083 }

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

Reimplemented from MadAnalysis.

Definition at line 441 of file MadTestAnalysis.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, MadBase::eventSummary, fLikeliHist, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRPlane::n, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpTrack, NtpSRTrack::ph, NtpSREvent::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, and NtpSRPulseHeight::sigcor.

Referenced by CreatePAN(), and PassAnalysisCuts().

00442 {
00443 
00444 
00445   if(!LoadEvent(event)) return -10;
00446   Int_t track = -1;
00447   if(!LoadLargestTrackFromEvent(event,track)) return -10;
00448   if(method!=0) return -10;
00449   
00450   Float_t ProbNC = 1.;
00451   Float_t ProbCC = 1.;
00452   Float_t PidCC = 0.;
00453   
00454 
00455   Float_t var1=eventSummary->plane.n;
00456   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00457     var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00458   }
00459   Float_t var2=0;
00460   Float_t var3=0;
00461 
00462 
00463   Float_t phtrack=ntpTrack->ph.sigcor;
00464   Float_t phevent=ntpEvent->ph.sigcor;
00465 
00466   if (phevent>0) {var2=phtrack/phevent;}
00467 
00468   if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor)
00469                               /float(ntpTrack->plane.n);
00470 
00471 
00472   Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00473   Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00474   Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00475 
00476   if (prob1==0) {prob1=0.0001;}
00477   if (prob2==0) {prob2=0.0001;}
00478   if (prob3==0) {prob3=0.0001;}
00479  
00480   ProbCC=prob1*prob2*prob3;
00481 
00482   prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00483   prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00484   prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00485  
00486   if (prob1==0) {prob1=0.0001;}
00487   if (prob2==0) {prob2=0.0001;}
00488   if (prob3==0) {prob3=0.0001;}
00489  
00490   ProbNC=prob1*prob2*prob3;
00491 
00492   PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00493 
00494   return PidCC;
00495 }

void MadTestAnalysis::ReadPIDFile ( std::string  tag  ) 

Definition at line 1398 of file MadTestAnalysis.cxx.

References fLikeliFile, and fLikeliHist.

01398                                               {
01399 
01400   fLikeliFile = new TFile(tag.c_str(),"READ");
01401  
01402   if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found
01403  
01404     fLikeliFile->cd();
01405 
01406     fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
01407     fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
01408     fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
01409     fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
01410     fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
01411     fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
01412     
01413     for(int i=0;i<6;i++){
01414       float integ = fLikeliHist[i]->Integral(); 
01415       fLikeliHist[i]->Scale(1./integ);
01416     }
01417   }
01418   else fLikeliFile = NULL;
01419 }

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

Reimplemented from MadAnalysis.

Definition at line 498 of file MadTestAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadQuantities::RecoMuEnergy(), and MadQuantities::RecoShwEnergy().

00498                                                       {
00499 
00500   if(!LoadEvent(event)) return 0;
00501   int largestTrack = -1;
00502   LoadLargestTrackFromEvent(event,largestTrack);
00503   float emu = RecoMuEnergy(0,largestTrack);
00504   float ehad = RecoShwEnergy(event);
00505   float ereco = emu + ehad; 
00506   return ereco;
00507 
00508 }

Double_t MadTestAnalysis::Sigmoid ( Double_t  x  )  [protected]

Definition at line 281 of file MadTestAnalysis.cxx.

Referenced by GetAnnPid().

00282 {
00283   double sig;
00284       if(x>37.){
00285          sig = 1.;
00286       }
00287       else if(x<-37.){
00288          sig = 0.;
00289       }
00290       else {
00291          sig = 1./(1.+exp(-x));
00292       }
00293       return sig;
00294 }

Float_t MadTestAnalysis::SingleTimeFrame ( Int_t  snarlentry,
Int_t  nentries 
) [protected]

Definition at line 699 of file MadTestAnalysis.cxx.

References MadBase::GetEntry(), RecPhysicsHeader::GetTimeFrame(), and MadBase::ntpHeader.

Referenced by CreatePAN().

00699                                                                        {
00700 
00701  Int_t tf=0,tfbef=0,tfaft=0;
00702    
00703    
00704    if(this->GetEntry(snarlentry))                                tf      = ntpHeader->GetTimeFrame();
00705    if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1))   tfaft   = ntpHeader->GetTimeFrame();   
00706    if(snarlentry>0        && this->GetEntry(snarlentry-1))       tfbef   = ntpHeader->GetTimeFrame();
00707  
00708  Float_t singletimeframe=0;
00709  if(snarlentry<(nentries-1) && snarlentry>0 ){
00710   if(tf!=tfaft && tf!=tfbef) singletimeframe=1;
00711  }
00712  if(snarlentry==(nentries-1)) {
00713   if(tf!=tfbef) singletimeframe=1; 
00714  }  
00715  if(snarlentry==0) {
00716   if(tf!=tfaft) singletimeframe=1;
00717  }
00718 
00719  this->GetEntry(snarlentry);
00720  return singletimeframe;  
00721  
00722 }


Member Data Documentation

Definition at line 50 of file MadTestAnalysis.h.

Referenced by CreatePAN().

TFile* MadTestAnalysis::fLikeliFile [protected]

Definition at line 28 of file MadTestAnalysis.h.

Referenced by MadTestAnalysis(), ReadPIDFile(), and ~MadTestAnalysis().

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

Definition at line 29 of file MadTestAnalysis.h.

Referenced by MadTestAnalysis(), PID(), and ReadPIDFile().

Definition at line 49 of file MadTestAnalysis.h.

Referenced by CreatePAN().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1