MadPIDAnalysis Class Reference

#include <MadPIDAnalysis.h>

Inheritance diagram for MadPIDAnalysis:
MadAnalysis MadQuantities MadBase

List of all members.

Public Member Functions

 MadPIDAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadPIDAnalysis (JobC *, string, int)
 ~MadPIDAnalysis ()
void MakeMyFile (std::string, int, int)
void MakeAbIDFile (std::string)
void ReadPIDFile (std::string)
void ReadAbPIDFile (std::string)
void ConfigureRoID (std::string, std::string)
void CreatePAN (std::string, Int_t, Int_t, EnergyCorrections::WhichCorrection_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)

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)
Bool_t PassFidND (Int_t event=0)
Bool_t PassFidNew (Int_t event=0)
Bool_t PassFidNewN (Int_t event=0)
Float_t RecoAnalysisEnergy (Int_t event=0)
Bool_t PassTimingCuts (Double_t timediff)
AnnInputBlock AnnVar (Int_t event=0)
Double_t GetAnnPid (AnnInputBlock anninput, Int_t det, Int_t tar, Int_t fid, Int_t prior, Bool_t first_ann)
Double_t Sigmoid (Double_t x)
Float_t SingleTimeFrame (Int_t snarlentry, Int_t nentries)

Protected Attributes

TFile * fLikeliFile
TH1F * fLikeliHist [6]
MadAbID andy_id
Anp::Interface ro_interface
Registry ro_registry

Private Attributes

Double_t Ann_weight_vector [226]
Float_t potall
Float_t potcut

Detailed Description

Definition at line 13 of file MadPIDAnalysis.h.


Constructor & Destructor Documentation

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

Definition at line 50 of file MadPIDAnalysis.cxx.

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

00052 {
00053 
00054   if(!chainSR) {
00055     record = 0;
00056     strecord = 0;
00057     emrecord = 0;
00058     mcrecord = 0;
00059     threcord = 0;
00060     Clear();
00061     whichSource = -1;
00062     isMC = true;
00063     isTH = true;
00064     isEM = true;
00065     Nentries = -1;
00066     return;
00067   }
00068   
00069   InitChain(chainSR,chainMC,chainTH,chainEM);
00070   whichSource = 0;
00071   fLikeliFile = NULL;
00072   for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00073 
00074   potall=0.;
00075   potcut=0.;
00076 
00077 }

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

Definition at line 80 of file MadPIDAnalysis.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.

00081 {
00082   record = 0;
00083   strecord = 0;
00084   emrecord = 0;
00085   mcrecord = 0;
00086   threcord = 0;
00087   Clear();
00088   isMC = true;
00089   isTH = true;
00090   isEM = true;
00091   Nentries = entries;
00092   whichSource = 1;
00093   jcPath = path;
00094   fJC = j;
00095   fLikeliFile = NULL;
00096   for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00097 }

MadPIDAnalysis::~MadPIDAnalysis (  ) 

Definition at line 100 of file MadPIDAnalysis.cxx.

References fLikeliFile, potall, and potcut.

00101 {
00102 
00103 
00104   if(fLikeliFile) fLikeliFile->Close();
00105 
00106   cout << "pot (all)     =" << potall << endl;
00107   cout << "pot (beamcuts)=" << potcut << endl;
00108 
00109 }


Member Function Documentation

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

Definition at line 257 of file MadPIDAnalysis.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::aShwstpu, AnnInputBlock::aShwstpv, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, 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, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower_Jim(), MadBase::LoadStrip(), NtpSRPlane::n, NtpSRShower::ndigit, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, 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, NtpSRStrip::strip, NtpSRStrip::time0, and NtpSRStrip::time1.

Referenced by CreatePAN().

00258 {
00259 
00260   AnnInputBlock anninput;
00261 
00262   // Initialize
00263 
00264   anninput.aTotrk       =0.;
00265   anninput.aTotstp      =0.;
00266   anninput.aTotph       =0.;
00267   anninput.aTotlen      =0.;
00268   anninput.aPhperpl     =0.;
00269   anninput.aPhperstp    =0.;
00270   
00271   
00272   anninput.aTrkpass     =0.;
00273   anninput.aTrkph       =0.;
00274   anninput.aTrklen      =0.;
00275   anninput.aTrkphperpl  =0.;
00276   anninput.aTrkphper    =0.;
00277   anninput.aTrkplu      =0.;
00278   anninput.aTrkplv      =0.;
00279   anninput.aTrkstp      =0.;
00280   
00281   anninput.aShwph       =0.;
00282   anninput.aShwstp      =0.;
00283   anninput.aShwdig      =0.;
00284   anninput.aShwpl       =0.;
00285   anninput.aShwphper    =0.;
00286   anninput.aShwphperpl  =0.;
00287   anninput.aShwphperdig =0.;
00288   anninput.aShwphperstp =0.;
00289   anninput.aShwplu      =0.;
00290   anninput.aShwplv      =0.;
00291   anninput.aShwstpu     =0.;
00292   anninput.aShwstpv     =0.;
00293   anninput.aPh3         =0.;
00294   anninput.aPh6         =0.;
00295   anninput.aPhlast      =0.;
00296   anninput.aPhcommon    =0.;
00297   anninput.aTimemax     =0.;
00298   anninput.aTimemin     =0.;
00299     
00300   
00301   // Calculate variables needed for ANN (and a few more)
00302 
00303   Int_t detector = 0;
00304 
00305   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) detector=1;
00306   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar)  detector=2;
00307   
00308 
00309   LoadEvent(event) ;  
00310   int track_index   = -1;
00311   LoadLargestTrackFromEvent(event,track_index);
00312   int shower_index  = -1;
00313   //  LoadLargestShowerFromEvent(event,shower_index);
00314   LoadShower_Jim(event,shower_index,detector);
00315 
00316   anninput.aTotrk       =ntpEvent->ntrack;
00317   anninput.aTotstp      =ntpEvent->nstrip;
00318   anninput.aTotph       =ntpEvent->ph.sigcor;
00319   anninput.aTotlen      =ntpEvent->plane.end-ntpEvent->plane.beg+1; // prepei na to alla3w
00320   
00321   if (ntpEvent->plane.n>0) anninput.aPhperpl     =ntpEvent->ph.sigcor/ntpEvent->plane.n;
00322   if (ntpEvent->nstrip>0) anninput.aPhperstp    =ntpEvent->ph.sigcor/ntpEvent->nstrip;
00323   
00324 
00325   if(track_index!=-1) {
00326     anninput.aTrkpass     =ntpTrack->fit.pass;
00327     anninput.aTrkph       =ntpTrack->ph.sigcor;
00328     anninput.aTrklen      =ntpTrack->plane.ntrklike;
00329 
00330     if(ntpTrack->plane.ntrklike>0) anninput.aTrkphperpl  =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
00331     if(ntpEvent->ph.sigcor>0)      anninput.aTrkphper    =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
00332 
00333     anninput.aTrkplu      =ntpTrack->plane.nu;
00334     anninput.aTrkplv      =ntpTrack->plane.nv;
00335     anninput.aTrkstp      =ntpTrack->nstrip;
00336   }
00337 
00338   if(shower_index!=-1) {
00339     anninput.aShwph       =ntpShower->ph.sigcor;
00340     anninput.aShwstp      =ntpShower->nstrip;
00341     anninput.aShwdig      =ntpShower->ndigit;
00342     anninput.aShwpl       =ntpShower->plane.end-ntpShower->plane.beg+1;
00343 
00344     if (ntpEvent->ph.sigcor>0) anninput.aShwphper    =ntpShower->ph.sigcor/ntpEvent->ph.sigcor;
00345     anninput.aShwphperpl  =ntpShower->ph.sigcor/(ntpShower->plane.end-ntpShower->plane.beg+1);
00346     
00347     if (ntpShower->ndigit>0) anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit;
00348     if (ntpShower->nstrip>0) anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip;
00349 
00350     anninput.aShwplu      =ntpShower->plane.nu;
00351     anninput.aShwplv      =ntpShower->plane.nv; 
00352   }
00353 
00354   
00355   Int_t ssind,ssplind;
00356   Double_t ssphtot;
00357   Bool_t foundshw,foundtrk;
00358   Int_t planes=ntpEvent->plane.beg;
00359   
00360   Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
00361   ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
00362   
00363   // loop over strips to compute ph3 ph6 phlast phcommon
00364 
00365   Double_t timemin=9.*10e30; 
00366   Double_t timemax=-9999999; 
00367   
00368   Float_t shwstpu=0;
00369   Float_t shwstpv=0;
00370   
00371   if(shower_index!=-1) {
00372     for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00373       Int_t stp_indexs=((ntpShower->stp)[jj]);
00374       if(stp_indexs!=-1){
00375         LoadStrip(stp_indexs);  
00376         if((ntpStrip->plane)%2==1) shwstpu=shwstpu+1;
00377         if((ntpStrip->plane)%2==0) shwstpv=shwstpv+1;
00378       }
00379     }       
00380   }      
00381   
00382   
00383   
00384   for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
00385     Int_t stp_index=((ntpEvent->stp)[evss]);
00386     if(stp_index!=-1){
00387       LoadStrip(stp_index);     
00388       // TIMING INFO // 
00389       //  Find max min time of events by loading strips for ND 
00390       if(detector==1){             
00391         Double_t striptime=ntpStrip->time1;
00392         if(striptime<=timemin) timemin=striptime;
00393         if(striptime>=timemax) timemax=striptime;                          
00394       } 
00395       
00396       if(detector==2){          
00397         Double_t striptime1=ntpStrip->time1;
00398         Double_t striptime0=ntpStrip->time0;
00399         Double_t striptime=0;
00400         if(striptime1>0 && striptime0<0)  striptime=striptime1;
00401         if(striptime0>0 && striptime1<0)  striptime=striptime0;
00402         if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
00403         if(striptime<=timemin) timemin=striptime;
00404         if(striptime>=timemax) timemax=striptime;
00405         
00406       }
00407       
00408       ssind=ntpStrip->strip;
00409       ssplind=ntpStrip->plane;
00410       ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;        
00411       
00412       foundshw=false;
00413       foundtrk=false;
00414       
00415       Double_t phstrips=0;
00416       Double_t phstript=0;       
00417  
00418       // shower strips
00419       Int_t sshwind,sshwplind;
00420       Double_t sshwphtot;
00421       
00422       if(shower_index!=-1) {
00423         for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00424           Int_t stp_indexs=((ntpShower->stp)[jj]);
00425           if(stp_indexs!=-1){
00426             LoadStrip(stp_indexs);      
00427             
00428             if(!foundshw){
00429               sshwind=ntpStrip->strip;
00430               sshwplind=ntpStrip->plane;
00431               sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;                 
00432               if(sshwind==ssind && sshwplind==ssplind) {
00433                 foundshw=true;
00434                 phstrips=sshwphtot;
00435               }
00436             }
00437           }
00438         }       
00439       }
00440       
00441       // tracks strips
00442       Int_t strkind,strkplind;
00443       Double_t strkphtot;
00444 
00445       if(track_index!=-1) {
00446         for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
00447           Int_t stp_indext=((ntpTrack->stp)[jj]);
00448           
00449           if(stp_indext!=-1){
00450             LoadStrip(stp_indext);      
00451             if(!foundtrk){
00452               strkind=ntpStrip->strip;
00453               strkplind=ntpStrip->plane;
00454               strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;       
00455               if(strkind==ssind && strkplind==ssplind) {
00456                 foundtrk=true;
00457                 phstript=strkphtot;
00458               }
00459             }
00460           }       
00461         }
00462       }
00463       
00464       if(foundshw && foundtrk) {
00465         hitcommon=hitcommon+1;
00466         phcommon=phcommon+phstrips+phstript;
00467       }  
00468 
00469       if(!foundshw && ! foundtrk) {
00470         hitnowere=hitnowere+1; 
00471         phnowere=phnowere+ssphtot; 
00472       }          
00473 
00474       if(ssplind>=planes && ssplind<=(planes+3)){
00475         ph3=ph3+ssphtot;
00476       }
00477       else if(ssplind>(planes+3) && ssplind<=(planes+6)){
00478         ph6=ph6+ssphtot; 
00479       }
00480       else {
00481         phlast=phlast+ssphtot;
00482       }     
00483 
00484     } // end if strip != -1 (!!!)
00485   } // end of looping over event strips....    
00486   
00487   anninput.aTimemin = timemin;
00488   anninput.aTimemax = timemax; 
00489   
00490   //
00491 
00492   anninput.aPh3       =ph3;
00493   anninput.aPh6       =ph6;
00494   anninput.aPhlast    =phlast;
00495   anninput.aPhcommon  =phcommon; 
00496   anninput.aShwstpu   =shwstpu;
00497   anninput.aShwstpv   =shwstpv;
00498   
00499   
00500   return anninput;
00501 
00502 }

void MadPIDAnalysis::ConfigureRoID ( std::string  keyfile,
std::string  tag 
)

Definition at line 2976 of file MadPIDAnalysis.cxx.

References Anp::Interface::Config(), ro_interface, ro_registry, Registry::Set(), and Registry::UnLockValues().

02976                                                                   {
02977    
02978   // DP updated feb 11th 2008 - new access method
02979 
02980   ro_registry.UnLockValues();
02981 
02982 
02983   ro_registry.Set("InterfaceConfigPath", keyfile.c_str());
02984   ro_registry.Set("FillkNNFilePath", tag.c_str());
02985 
02986 
02987   // configure RoID interface with registry values
02988 
02989   ro_interface.Config(ro_registry);
02990 
02991 }

void MadPIDAnalysis::CreatePAN ( std::string  tag,
Int_t  scanfilter,
Int_t  mcneartype,
EnergyCorrections::WhichCorrection_t  corrver 
)

Definition at line 1090 of file MadPIDAnalysis.cxx.

References andy_id, AnnVar(), 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::aShwstpu, AnnInputBlock::aShwstpv, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, 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::begu, NtpSRPlane::begv, NtpSRMomentum::best, MadAbID::CalcPID(), NtpSRFitTrack::chi2, NtpSRDetStatus::coilcurrent1, NtpSRDetStatus::coilcurrent2, NtpSRDataQuality::coldchips, NtpTHEvent::completeall, NtpTHEvent::completeslc, NtpSRTrack::contained, NtpSRDataQuality::cratemask, NtpStRecord::dataquality, NtpSREventSummary::date, NtpSRDate::day, Munits::day, NtpSRVertex::dcosx, NtpSRVertex::dcosy, NtpSRVertex::dcosz, NtpStRecord::detstatus, NtpSRTrack::ds, NtpSRPlane::end, NtpSREvent::end, NtpSRTrack::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, MadBase::eventSummary, BeamMonSpill::fHornCur, Anp::Interface::FillSnarl(), NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, NtpMCFluxInfo::fluxevtno, NtpMCFluxInfo::fluxrun, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTrtgtd, BDSpillAccessor::Get(), GetAnnPid(), Dcs_Mag_Near::GetCurrent(), ScanList::GetDecision(), VldContext::GetDetector(), MadBase::GetEntry(), HvStatusFinder::GetHvStatus(), CoilTools::GetMagNear(), VldTimeStamp::GetNanoSec(), SpillTimeFinder::GetNearestSpill(), SpillServerMonFinder::GetNearestSpill(), MadAbID::GetRecoCosTheta(), MadAbID::GetRecoE(), MadAbID::GetRecoY(), NtpStRecord::GetRelease(), RecDataHeader::GetRun(), VldTimeStamp::GetSec(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), SpillServerMon::GetSpillTime(), SpillServerMon::GetSpillTimeError(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), SpillTimeFinder::GetTimeOfNearestSpill(), SpillTimeND::GetTimeStamp(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), MadAbID::GetTrackEMCharge(), MadAbID::GetTrackLikePlanes(), MadAbID::GetTrackPHfrac(), MadAbID::GetTrackPHmean(), MadAbID::GetTrackPlanes(), MadAbID::GetTrackQPsigmaQP(), RecPhysicsHeader::GetTrigSrc(), Anp::Interface::GetVar(), RecHeader::GetVldContext(), HvStatus::Good(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSREvent::index, MadQuantities::Initial_state(), SpillServerMonFinder::Instance(), CoilTools::Instance(), HvStatusFinder::Instance(), SpillTimeFinder::Instance(), MadQuantities::INu(), MadQuantities::IResonance(), MadQuantities::IsFid_2008(), MadQuantities::IsFidAll(), MadQuantities::IsFidFD_AB(), MadQuantities::IsFidVtxEvt(), DataUtil::IsGoodFDData(), LISieve::IsLI(), RunStatus::IsOK(), CoilTools::IsReverse(), Detector::kFar, Detector::kNear, NtpSRShowerPulseHeight::linCCgev, NtpSRShowerPulseHeight::linNCgev, NtpSREventSummary::litime, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), BDSpillAccessor::LoadSpill(), MadBase::LoadStrip(), MadBase::LoadTHEvent(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), NtpSRTrack::momentum, month, NtpSRDate::month, NtpSRPlane::n, NtpSRShower::ncluster, NtpSRTrack::ndigit, NtpSRShower::ndigit, NtpSREvent::ndigit, NtpSRFitTrack::ndof, NtpMCFluxInfo::nenergyfar, NtpMCFluxInfo::nenergynear, MadBase::Nentries, NtpSREventSummary::nevent, NtpMCFluxInfo::nimpwt, NtpSREventSummary::nshower, NtpSREvent::nshower, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTHEvent, MadBase::ntpTrack, MadBase::ntpTruth, NtpSREvent::ntrack, NtpSRPlane::nu, MadQuantities::Nucleus(), NtpSRPlane::nv, NtpMCFluxInfo::nwtfar, NtpMCFluxInfo::nwtnear, NtpSRFitTrack::pass, PassAnalysisCuts(), PassBasicCuts(), passfail(), PassFid(), PassFidND(), PassFidNew(), PassFidNewN(), PassTimingCuts(), NtpSREventSummary::ph, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, PID(), NtpSRShower::plane, NtpSRStrip::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, potall, potcut, NtpTHEvent::purity, MadQuantities::Q2(), NtpSRMomentum::qp, NtpSRMomentum::range, NtpSRPulseHeight::raw, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergyNew(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergyNew(), MadQuantities::RecoShwEnergySqrt(), ro_interface, run(), BMSpillAna::SelectSpill(), BMSpillAna::SetSnarlTime(), BMSpillAna::SetSpill(), NtpSRShower::shwph, NtpSRPulseHeight::sigcor, NtpSRDataQuality::spillstatus, BeamMonSpill::SpillTime(), NtpSRDataQuality::spilltimeerror, NtpSRDataQuality::spilltype, NtpSRShower::stp, MadBase::strecord, MadQuantities::Target4Vector(), NtpMCFluxInfo::tptype, NtpMCFluxInfo::tpx, NtpMCFluxInfo::tpy, NtpMCFluxInfo::tpz, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSRVertex::u, BMSpillAna::UseDatabaseCuts(), NtpSRVertex::v, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, MadQuantities::W2(), NtpSRShowerPulseHeight::wtCCgev, NtpSRShowerPulseHeight::wtNCgev, MadQuantities::X(), NtpSRVertex::x, MadQuantities::Y(), NtpSRVertex::y, NtpSRDate::year, Munits::year, NtpSRVertex::z, and SpillInfoBlock::Zero().

01090                                                                                                                            {
01091 
01092 // 1 = le-10 2 = pme 3 = phe there is no other way in an MC file to now what target positions it
01093 // was...!
01094   
01095   std::string savename = "PAN_" + tag + ".root";
01096   TFile save(savename.c_str(),"RECREATE"); 
01097   save.cd();
01098 
01099 
01100   SpillInfo pppinfo;
01101  
01102 
01103   // data quality variables from the standard ntuple
01104 
01105   NtpSRDataQuality *ntpDataQual;
01106   NtpSRDetStatus *ntpDetStatus;
01107 
01108   
01109  // Beam monitoring object to apply beam data quality cuts
01110   
01111   BMSpillAna bm;
01112   bm.UseDatabaseCuts();
01113 
01114  
01115   //ann quantities
01116 
01117   Float_t ann_aTotrk       =0.;
01118   Float_t ann_aTotstp      =0.;
01119   Float_t ann_aTotph       =0.;
01120   Float_t ann_aTotlen      =0.;
01121   Float_t ann_aPhperpl     =0.;
01122   Float_t ann_aPhperstp    =0.;
01123   
01124   
01125   Float_t ann_aTrkpass     =0.;
01126   Float_t ann_aTrkph       =0.;
01127   Float_t ann_aTrklen      =0.;
01128   Float_t ann_aTrkphperpl  =0.;
01129   Float_t ann_aTrkphper    =0.;
01130   Float_t ann_aTrkplu      =0.;
01131   Float_t ann_aTrkplv      =0.;
01132   Float_t ann_aTrkstp      =0.;
01133   
01134   Float_t ann_aShwph       =0.;
01135   Float_t ann_aShwstp      =0.;
01136   Float_t ann_aShwdig      =0.;
01137   Float_t ann_aShwpl       =0.;
01138   Float_t ann_aShwphper    =0.;
01139   Float_t ann_aShwphperpl  =0.;
01140   Float_t ann_aShwphperdig =0.;
01141   Float_t ann_aShwphperstp =0.;
01142   Float_t ann_aShwplu      =0.;
01143   Float_t ann_aShwplv      =0.;
01144   Float_t ann_aShwstpu     =0.; 
01145   Float_t ann_aShwstpv     =0.;  
01146   Float_t ann_aPh3         =0.;
01147   Float_t ann_aPh6         =0.;
01148   Float_t ann_aPhlast      =0.;
01149   Float_t ann_aPhcommon    =0.;
01150 
01151   //PAN Quantities 
01152   //Truth:
01153   Float_t true_enu;       // true neutrino energy (GeV)
01154   Float_t true_pxnu;      // true neutrino momentum-x (GeV)
01155   Float_t true_pynu;      // true neutrino momentum-y (GeV)
01156   Float_t true_pznu;      // true neutrino momentum-z (GeV)
01157   Float_t true_etgt;       // true target energy (GeV)
01158   Float_t true_pxtgt;      // true target momentum-x (GeV)
01159   Float_t true_pytgt;      // true target momentum-y (GeV)
01160   Float_t true_pztgt;      // true target momentum-z (GeV)
01161   Float_t true_emu;       // true muon energy
01162   Float_t true_eshw;      // true shower energy
01163   Float_t true_x;         // true x
01164   Float_t true_y;         // true y
01165   Float_t true_q2;        // true q2
01166   Float_t true_w2;        // true w2
01167   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
01168   Float_t true_dircosz;   // track z-direction cosine
01169   Float_t true_vtxx;      // true x vtx
01170   Float_t true_vtxy;      // true y vtx
01171   Float_t true_vtxz;      // true z vtx
01172 
01173   //For Neugen:
01174   Int_t flavour;          // true flavour: 1-e 2-mu 3-tau
01175   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
01176   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
01177   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
01178   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
01179   Int_t had_fs;           // hadronic final state, number between 200-300
01180   
01181   //For Beam Reweighting:
01182  // NuParent *nuparent = new NuParent();
01183 
01184   //Reco Quantities
01185   Int_t detector;         // Near = 1, Far = 2;
01186   Int_t run;              // run number
01187   Int_t snarl;            // snarl number
01188   Int_t nevent;           // number of events in snarl
01189   Int_t event;            // event index
01190   Int_t subrun;           // subrun number     
01191   Int_t trigsrc;          // trigger source    
01192   Int_t mcevent;          // mc event index
01193   Int_t ntrack;           // number of tracks in event
01194   Int_t nshower;          // number of showers in event
01195 
01196   Int_t is_fid;           // pass fid cut. true = 1, false = 0
01197   Int_t is_fid_dp;        // pass dp analysis fid cut. true = 1, false = 0
01198   Int_t is_fid_ns;        // pass  nd analysis fid cut true=1, false=0 // !! NEW
01199   Int_t is_fid_ab;        // andy's FD fiducial cuts - 2007
01200   Int_t is_fid_2008;      // 2008 fiducial cuts
01201   
01202 
01203   Int_t is_cev;           // fully contained. true = 1, false = 0 
01204   Int_t is_cont_trk;      // new containemnt support by IsContained() for the largest track
01205   Int_t is_mc;            // is a corresponding mc event found
01206   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
01207   Int_t pass_cccuts;      // FD: ntrack, fid and trkdircos requirement
01208   Int_t pass_cccuts_2007; // FD: as above, but with Andy's Fidvol
01209   Int_t pass_cccuts_2008; // FD: as above, but with Andy's Fidvol
01210 
01211   Float_t pid0;           // pid parameter (DP)
01212   Float_t pid1;           // pid parameter (AB)
01213   Float_t pid2;           // pid parameter (RO)
01214   Float_t annpid;         // ANN PID
01215   Float_t annpid_f1p1;    // ANN PID FAR FIDUCIAL DP PRIOR 1 // UN-OSICLLATED 
01216   Float_t annpid_f1p2;    // ANN PID FAR FIDUCIAL DP PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95
01217   Float_t annpid_f1p3;    // ANN PID FAR FIDUCIAL DP PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95
01218   Float_t annpid_f2p1;    // ANN PID FAR FIDUCIAL NS PRIOR 1  UN-OSICLLATED 
01219   Float_t annpid_f2p2;    // ANN PID FAR FIDUCIAL NS PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95
01220   Float_t annpid_f2p3;    // ANN PID FAR FIDUCIAL NS PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95
01221   
01222   Int_t pass;             // pass analysis cuts. true = 1, false = 0
01223   Int_t passcut;
01224 
01225   Float_t reco_enu;       // reco neutrino energy
01226   Float_t reco_emu;       // reco muon energy
01227   Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
01228   
01229   Int_t pass_fd_qualcuts;  // pass FD quality cuts (HV trips etc)
01230   Int_t litag;  
01231   Int_t lisieve;           // IS LI - LISieve
01232 
01233 
01234 
01235   Float_t  reco_eshwcc;    // linear cc version
01236   Float_t  reco_eshwnc;    // linear nc version 
01237   Float_t  reco_eshwccw;   // weighted cc version 
01238   Float_t  reco_eshwncw;   // weighted nc version 
01239 
01240   Float_t  reco_eshwcc_uncorr;    // linear cc version - uncorrected
01241   Float_t  reco_emu_uncorr;
01242 
01243  
01244   Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
01245   Float_t reco_qe_enu;    // reco qe neutrino energy
01246   Float_t reco_qe_q2;     // reco qe q2
01247   Float_t reco_y;         // reco y
01248   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
01249   Float_t reco_dircosz;   // reco track vtx z-dircos
01250   Float_t reco_vtxx;      // reco x vtx
01251   Float_t reco_vtxy;      // reco y vtx
01252   Float_t reco_vtxz;      // reco z vtx
01253 
01254   Float_t evtlength;      // reco event length
01255   Float_t evtph;          // reco event ph
01256   Float_t trklength;      // reco track length
01257   Float_t trkmomrange;    // reco track momentum from range
01258   Float_t trkqp;          // reco track q/p
01259   Float_t trkeqp;         // reco track q/p error
01260   Float_t trkmombest;     // reco best track momentum based on new containment support
01261   Float_t trkphfrac;      // reco pulse height fraction in track
01262   Float_t trkphplane;     // reco average track pulse height/plane
01263   Float_t trkds;          // track DS 
01264   Float_t rawph;          // Raw ph
01265 // Additional track info
01266   Float_t trkendz;  // track end Z position
01267   Float_t trkendx;  // track end X position
01268   Float_t trkendy;  // track end Y position
01269   Float_t trkendu;  // track end U position
01270   Float_t trkendv;  // track end V position
01271   Float_t trkvtxz;  // track begin Z position
01272   Float_t trkvtxx;  // track begin X position
01273   Float_t trkvtxy;  // track begin Y position
01274   Float_t trkvtxu;  // track begin U position
01275   Float_t trkvtxv;  // track begin V position
01276   Float_t trkplbu;  // track begin plane u 
01277   Float_t trkplbv;  // track begin plane v 
01278   Float_t trkpleu;  // track end   plane u 
01279   Float_t trkplev;  // track end   plane v 
01280   Float_t trkple;   // track end   plane  
01281   Float_t trkplb;   // track begin   plane 
01282   Float_t trkcosx;        // track x direction cosine
01283   Float_t trkcosy;        // track y direction cosine
01284   Float_t trkcosz;        // track z direction cosine
01285   Float_t evanglet;
01286   Int_t trkfitpass;     //
01287   Int_t totdigits,totstrips;
01288   Int_t trkdigits,trkstrips;
01289   Float_t trkph,trkchi2,trkndof;
01290   Int_t trkdiffuv;
01291 
01292   // shower variables
01293   Float_t shwph;
01294   Int_t   shwdigits, shwstrips;
01295   Float_t shwvtxz;  // shower begin Z position
01296   Float_t shwvtxx;  // shower begin X position
01297   Float_t shwvtxy;  // shower begin Y position
01298   Float_t shwvtxu;  // shower begin U position
01299   Float_t shwvtxv;  // shower begin V position
01300  
01301   // cluster variables
01302 
01303    Int_t shwncluster;
01304   
01305   // SCAN DECISION   // <- define an instance of ScanList class here
01306    Int_t scandecision;
01307    ScanList scan;
01308   
01309    // EVENT TIMING
01310    Double_t evttimemin;  //MIN STRIP TIME OF EVENT
01311    Double_t evttimemax;  //MAX STRIP TIME OF EVENT
01312    Double_t snarlt;
01313    Double_t litime; 
01314  
01315    // BEAM INFO
01316    Float_t  pot,horn,tar,vvpos,hhpos,magn;    
01317    // Beam INFO DATABASE
01318    Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
01319    Double_t timedb;
01320    // TIME INFO DATABASE
01321    Double_t neartdb,fartdb,neardifftdb;
01322 
01323    // mc flux info for xf,pt weighting
01324    Float_t mcflux_tpx;
01325    Float_t mcflux_tpy;
01326    Float_t mcflux_tpz;
01327    Float_t mcflux_tptype;
01328    Float_t mcflux_nenergyN;
01329    Float_t mcflux_nenergyF;
01330    Float_t mcflux_nwtN;
01331    Float_t mcflux_nwtF;
01332    Float_t mcflux_nimpwt;
01333 
01334    Int_t mcflux_fluxrun;
01335    Int_t mcflux_fluxevtno;
01336 
01337    // inu from MC
01338 
01339    Int_t inu;
01340 
01341 
01342    // ANDY's PID inputs
01343 
01344    Float_t abid_costheta;
01345    Float_t abid_eventenergy;
01346    Int_t abid_trackcharge;
01347    Int_t abid_trackenergy;
01348    Int_t abid_tracklikeplanes;
01349    Float_t abid_trackphfrac;
01350    Float_t abid_trackphmean;
01351    Float_t abid_trackqpsigmaqp;
01352    Float_t abid_y;
01353 
01354 
01355    // RUSTEM's PID inputs
01356 
01357 
01358    Float_t roid_numplanes;      // number of scintillator planes
01359    Float_t roid_meansigcor;     // mean sigcor ADC
01360    Float_t roid_lowphhiph;      // fraction of "lower ph"/"higher ph"
01361    Float_t roid_trkphfracclose; // fraction of track ph/ph of strips
01362 
01363 
01364    // not used
01365    //   Float_t roid_trkplanefrac;   // frac of track planes/total planes
01366                                      // around track
01367    //   Float_t roid_trkphfracend;   // fraction of track ph at end/track ph
01368 
01369 
01370 
01371    // new data quality variables
01372 
01373    Int_t pass_fd_qualcuts_db;   // data quality flag from DB quantities
01374 
01375    Int_t dataqual_cratemask;    // number of readout crates
01376    Int_t dataqual_spilltype;    // spill type 0-no, 1-true, 3-fake
01377    Int_t dataqual_spillerror;   // GPS error in ns (sum in quad ND+FD)
01378    Int_t dataqual_spillstatus;  // SpillServer status - 1=OK
01379    Int_t dataqual_coldchips;    // number of cold chips in FD (<20kHZ)
01380 
01381    Float_t dataqual_coilcurr_sm1; // SM1 coil current (A)  
01382    Float_t dataqual_coilcurr_sm2; // SM2 coil current (A)
01383 
01384 
01385    Int_t alec_hvstatus;    // HV status from offline DB   1=good
01386    Int_t alec_coilstatus;  // coil status from offline DV 1=good
01387    Int_t alec_spillstatus; // spillserver status from monitoring
01388                            // blocks - 1=good
01389 
01390 
01391 
01392  
01393    Int_t pass_fd_timing; // passes spill timing cuts 
01394 
01395    Int_t nd_reclaim;     // ND reclamation cuts
01396 
01397 
01398    Int_t coilok;         // coil current OK
01399    Int_t coilpol;        // coil polarity 1 - forward, -1 reversed
01400 
01401    Int_t numplanes;      // evt.plane.n
01402 
01403 
01404    Float_t evtcom;       // reco completeness and purity variables
01405    Float_t evtslc;
01406    Float_t evtpur;
01407 
01408 
01409 
01410 
01411 // IS SNARL CUT IN PIECES 
01412    Float_t singletimeframe;
01413 // Beam Info Block    
01414    SpillInfoBlock *spinfo = new SpillInfoBlock();  
01415    spinfo->Zero(); 
01416    singletimeframe=-999;
01417  //
01418   Int_t day,month,year;
01419   Int_t pass_beamcuts;     // pass default beam quality cuts   
01420 
01421 
01422 
01423   //Trees
01424 
01425   TTree *tree = new TTree("pan","pan");
01426 
01427   //Truth
01428   tree->Branch("true_enu",&true_enu,"true_enu/F",512000);
01429   tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000);
01430   tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000);
01431   tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000);
01432   tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000);
01433   tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000);
01434   tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000);
01435   tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000);
01436   tree->Branch("true_emu",&true_emu,"true_emu/F",512000);
01437   tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000);
01438   tree->Branch("true_x",&true_x,"true_x/F",512000);
01439   tree->Branch("true_y",&true_y,"true_y/F",512000);
01440   tree->Branch("true_q2",&true_q2,"true_q2/F",512000);
01441   tree->Branch("true_w2",&true_w2,"true_w2/F",512000);
01442   tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000);
01443   tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000);
01444   tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000);
01445   tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000);
01446   tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000);
01447   //Neugen
01448   tree->Branch("flavour",&flavour,"flavour/I",512000);
01449   tree->Branch("process",&process,"process/I",512000);
01450   tree->Branch("nucleus",&nucleus,"nucleus/I",512000);
01451   tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000);
01452   tree->Branch("initial_state",&initial_state,"initial_state/I",512000);
01453   tree->Branch("had_fs",&had_fs,"had_fs/I",512000);
01454   //NuParent
01455   //tree->Branch("nuparent","NuParent",&nuparent,8000,1);
01456   //Reco
01457   tree->Branch("detector",&detector,"detector/I",512000);
01458   tree->Branch("run",&run,"run/I",512000);
01459   tree->Branch("snarl",&snarl,"snarl/I",512000);
01460   tree->Branch("event",&event,"event/I",512000);
01461   tree->Branch("mcevent",&mcevent,"mcevent/I",512000);
01462   tree->Branch("ntrack",&ntrack,"ntrack/I",512000);
01463   tree->Branch("nshower",&nshower,"nshower/I",512000);
01464   tree->Branch("subrun",&subrun,"subrun/I",512000);   
01465   tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000); 
01466   tree->Branch("litime",&litime,"litime/D",512000);    
01467   tree->Branch("is_fid",&is_fid,"is_fid/I",512000);
01468   tree->Branch("is_fid_dp",&is_fid_dp,"is_fid_dp/I",512000);
01469   tree->Branch("is_fid_ns",&is_fid_ns,"is_fid_ns/I",512000); 
01470   
01471 
01472   tree->Branch("pass_fd_qualcuts",  &pass_fd_qualcuts, "pass_fd_qualcuts/I");
01473   tree->Branch("litag",             &litag,            "litag/I");
01474   tree->Branch("lisieve",             &lisieve,            "lisieve/I");
01475 
01476 
01477   tree->Branch("is_cev",&is_cev,"is_cev/I",512000);
01478   tree->Branch("is_cont_trk",&is_cont_trk,"is_cont_trk/I",512000);
01479   tree->Branch("is_mc",&is_mc,"is_mc/I",512000);
01480   tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000);
01481   tree->Branch("pass_cccuts",&pass_cccuts,"pass_cccuts/I",512000);
01482   tree->Branch("pass_cccuts_2007",&pass_cccuts_2007,"pass_cccuts_2007/I",512000);
01483   tree->Branch("pass_cccuts_2008",&pass_cccuts_2008,"pass_cccuts_2008/I",512000);
01484   tree->Branch("pid0",&pid0,"pid0/F",512000);
01485   tree->Branch("pid1",&pid1,"pid1/F",512000);
01486   tree->Branch("pid2",&pid2,"pid2/F",512000);
01487   tree->Branch("annpid",&annpid,"annpid/F",512000);    
01488   tree->Branch("annpid_f1p1",&annpid_f1p1,"annpid_f1p1/F",512000);  
01489   tree->Branch("annpid_f1p2",&annpid_f1p2,"annpid_f1p2/F",512000); 
01490   tree->Branch("annpid_f1p3",&annpid_f1p3,"annpid_f1p3/F",512000); 
01491   tree->Branch("annpid_f2p1",&annpid_f2p1,"annpid_f2p1/F",512000); 
01492   tree->Branch("annpid_f2p2",&annpid_f2p2,"annpid_f2p2/F",512000); 
01493   tree->Branch("annpid_f2p3",&annpid_f2p3,"annpid_f2p3/F",512000);   
01494   tree->Branch("pass",&pass,"pass/I",512000);
01495   tree->Branch("passcut",           &passcut,          "passcut/I");
01496 
01497   tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000);
01498   tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000);
01499   tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000);
01500   tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000);  
01501   tree->Branch("reco_eshwcc",  &reco_eshwcc,  "reco_eshwcc/F",512000);  
01502   tree->Branch("reco_eshwcc_uncorr",  &reco_eshwcc_uncorr,  "reco_eshwcc_uncorr/F",512000);  
01503   tree->Branch("reco_emu_uncorr",  &reco_emu_uncorr,  "reco_emu_uncorr/F",512000);  
01504   tree->Branch("reco_eshwnc",  &reco_eshwnc,  "reco_eshwnc/F",512000);  
01505   tree->Branch("reco_eshwccw", &reco_eshwccw, "reco_eshwccw/F",512000); 
01506   tree->Branch("reco_eshwncw", &reco_eshwncw, "reco_eshwncw/F",512000);
01507   tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000);
01508   tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000);
01509   tree->Branch("reco_y",&reco_y,"reco_y/F",512000);
01510   tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000);
01511   tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000);
01512   tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000);
01513   tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000);
01514   tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000);
01515 
01516   tree->Branch("evtlength",&evtlength,"evtlength/F",512000);
01517   tree->Branch("evtph",&evtph,"evtph/F",512000);
01518   tree->Branch("trklength",&trklength,"trklength/F",512000);
01519   tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000);
01520   tree->Branch("trkqp",&trkqp,"trkqp/F",512000);
01521   tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000);
01522   tree->Branch("trkmombest",&trkmombest,"trkmombest/F",512000);
01523   tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000);
01524   tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000);
01525   tree->Branch("rawph",&rawph,"rawph/F",512000);   
01526   tree->Branch("numplanes",&numplanes,"numplanes/I",512000);
01527    
01528   tree->Branch("trkendz",&trkendz,"trkendz/F",512000);
01529   tree->Branch("trkendu",&trkendu,"trkendu/F",512000);
01530   tree->Branch("trkendv",&trkendv,"trkendv/F",512000);
01531   tree->Branch("trkendx",&trkendx,"trkendx/F",512000);
01532   tree->Branch("trkendy",&trkendy,"trkendy/F",512000);
01533   tree->Branch("trkplbu",&trkplbu,"trkplbu/F",512000);
01534   tree->Branch("trkplbv",&trkplbv,"trkplbv/F",512000);
01535   tree->Branch("trkpleu",&trkpleu,"trkpleu/F",512000);
01536   tree->Branch("trkplev",&trkplev,"trkplev/F",512000);
01537   tree->Branch("trkple",&trkple,"trkple/F",512000); 
01538   tree->Branch("trkplb",&trkplb,"trkplb/F",512000); 
01539   tree->Branch("trkcosx",&trkcosx,"trkcosx/F",512000);
01540   tree->Branch("trkcosy",&trkcosy,"trkcosy/F",512000);
01541   tree->Branch("trkcosz",&trkcosz,"trkcosz/F",512000);
01542   
01543   tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000);
01544   tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000);
01545   tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000);
01546   tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000);
01547   tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000);
01548   tree->Branch("trkds",  &trkds,"trkds/F",512000);
01549   tree->Branch("evanglet",&evanglet,"evanglet/F",512000);
01550   
01551   tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
01552   tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
01553  
01554   tree->Branch("totdigits",&totdigits,"totdigits/I");
01555   tree->Branch("totstrips",&totstrips,"totstrips/I");
01556   tree->Branch("trkdigits",&trkdigits,"trkdigits/I");
01557   tree->Branch("trkstrips",&trkstrips,"trkstrips/I");
01558   tree->Branch("trkph",&trkph,"trkph/F");
01559   tree->Branch("trkdiffuv",&trkdiffuv,"trkdiffuv/I",512000);  
01560  
01561   tree->Branch("trkchi2",&trkchi2,"trkchi2/F");
01562   tree->Branch("trkndof",&trkndof,"trkndof/F");
01563   tree->Branch("trkfitpass",&trkfitpass,"trkfitpass/I",512000);  
01564 
01565   tree->Branch("shwdigits",&shwdigits,"shwdigits/I");
01566   tree->Branch("shwstrips",&shwstrips,"shwstrips/I");
01567   tree->Branch("shwph",&shwph,"shwph/F");
01568 
01569   tree->Branch("shwvtxz",&shwvtxz,"shwvtxz/F",512000);
01570   tree->Branch("shwvtxu",&shwvtxu,"shwvtxu/F",512000);
01571   tree->Branch("shwvtxv",&shwvtxv,"shwvtxv/F",512000);
01572   tree->Branch("shwvtxx",&shwvtxx,"shwvtxx/F",512000);
01573   tree->Branch("shwvtxy",&shwvtxy,"shwvtxy/F",512000);
01574   tree->Branch("shwncluster",&shwncluster,"shwncluster/I");   
01575    
01576   // SCAN DECISIONS 
01577   tree->Branch("scandecision",&scandecision,"scandecision/I");
01578 
01579   //BEAM INFO 
01580   tree->Branch("pot",    &pot,    "pot/F",512000);
01581   tree->Branch("horn",   &horn,   "horn/F",512000);
01582   tree->Branch("tar",    &tar,    "tar/F",512000);
01583   tree->Branch("vvpos",  &vvpos,  "vvpos/F",512000);
01584   tree->Branch("hhpos",  &hhpos,  "hhpos/F",512000);
01585   tree->Branch("magn",   &magn,   "magn/F",512000);
01586   tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000);
01587 
01588 
01589 
01590   // BEAM INFO DB
01591   tree->Branch("potdb",     &potdb,      "potdb/F",512000);
01592   tree->Branch("horndb",    &horndb,     "horndb/F",512000);
01593   tree->Branch("tardb",     &tardb,      "tardb/F",512000);
01594   tree->Branch("vvposdb",   &vvposdb,    "vvposdb/F",512000);
01595   tree->Branch("hhposdb",   &hhposdb,    "hhposdb/F",512000);
01596   tree->Branch("vvwidthdb", &vvwidthdb,  "vvwidthdb/F",512000);
01597   tree->Branch("hhwidthdb", &hhwidthdb,  "hhwidthdb/F",512000);
01598   tree->Branch("timedb",    &timedb,     "timedb/D");
01599 
01600   // TIME DB  
01601   tree->Branch("neartdb",      &neartdb,    "neartdb/D");
01602   tree->Branch("fartdb",       &fartdb,     "fartdb/D");
01603   tree->Branch("neardifftdb",  &neardifftdb,"neardifftdb/D"); 
01604   tree->Branch("snarlt",       &snarlt,     "snarlt/D");
01605  // mc flux info
01606 
01607   tree->Branch("mcflux_tpx", &mcflux_tpx, "mcflux_tpx/F");  
01608   tree->Branch("mcflux_tpy", &mcflux_tpy, "mcflux_tpy/F");  
01609   tree->Branch("mcflux_tpz", &mcflux_tpz, "mcflux_tpz/F");  
01610   tree->Branch("mcflux_tptype", &mcflux_tptype, "mcflux_tptype/F");  
01611   tree->Branch("mcflux_nenergyN", &mcflux_nenergyN, "mcflux_nenergyN/F");  
01612   tree->Branch("mcflux_nenergyF", &mcflux_nenergyF, "mcflux_nenergyF/F");  
01613   tree->Branch("mcflux_nwtN", &mcflux_nwtN, "mcflux_nwtN/F");  
01614   tree->Branch("mcflux_nwtF", &mcflux_nwtF, "mcflux_nwtF/F");  
01615   tree->Branch("mcflux_nimpwt", &mcflux_nimpwt, "mcflux_nimpwt/F");  
01616   tree->Branch("mcflux_fluxrun", &mcflux_fluxrun, "mcflux_fluxrun/I");  
01617   tree->Branch("mcflux_fluxevtno", &mcflux_fluxevtno, "mcflux_fluxevtno/I"); 
01618   
01619 // ann vars
01620 
01621    tree->Branch("ann_aTotrk",       &ann_aTotrk, "ann_aTotrk/F");
01622    tree->Branch("ann_aTotstp",      &ann_aTotstp, "ann_aTotstp/F");
01623    tree->Branch("ann_aTotph",       &ann_aTotph, "ann_aTotph/F");
01624    tree->Branch("ann_aTotlen",      &ann_aTotlen, "ann_aTotlen/F");
01625    tree->Branch("ann_aPhperpl",     &ann_aPhperpl, "ann_aPhperpl/F");
01626    tree->Branch("ann_aPhperstp",    &ann_aPhperstp, "ann_aPhperstp/F"); 
01627    tree->Branch("ann_aTrkpass",     &ann_aTrkpass, "ann_aTrkpass/F");  
01628    tree->Branch("ann_aTrkph",       &ann_aTrkph, "ann_aTrkph/F");    
01629    tree->Branch("ann_aTrklen",      &ann_aTrklen, "ann_aTrklen/F");   
01630    tree->Branch("ann_aTrkphperpl",  &ann_aTrkphperpl, "ann_aTrkphperpl/F");
01631    tree->Branch("ann_aTrkphper",    &ann_aTrkphper, "ann_aTrkphper/F"); 
01632    tree->Branch("ann_aTrkplu",      &ann_aTrkplu, "ann_aTrkplu/F");   
01633    tree->Branch("ann_aTrkplv",      &ann_aTrkplv, "ann_aTrkplv/F");   
01634    tree->Branch("ann_aTrkstp",      &ann_aTrkstp, "ann_aTrkstp/F");    
01635    tree->Branch("ann_aShwph",       &ann_aShwph, "ann_aShwph/F");      
01636    tree->Branch("ann_aShwstp",      &ann_aShwstp, "ann_aShwstp/F");   
01637    tree->Branch("ann_aShwdig",      &ann_aShwdig, "ann_aShwdig/F");     
01638    tree->Branch("ann_aShwpl",       &ann_aShwpl,      "ann_aShwpl/F");      
01639    tree->Branch("ann_aShwphper",    &ann_aShwphper,   "ann_aShwphper/F");  
01640    tree->Branch("ann_aShwphperpl",  &ann_aShwphperpl, "ann_aShwphperpl/F"); 
01641    tree->Branch("ann_aShwphperdig", &ann_aShwphperdig,"ann_aShwphperdig/F");
01642    tree->Branch("ann_aShwphperstp", &ann_aShwphperstp,"ann_aShwphperstp/F");
01643    tree->Branch("ann_aShwplu",      &ann_aShwplu,     "ann_aShwplu/F");     
01644    tree->Branch("ann_aShwplv",      &ann_aShwplv,     "ann_aShwplv/F");    
01645    tree->Branch("ann_aShwstpu",     &ann_aShwstpu,    "ann_aShwstpu/F");    
01646    tree->Branch("ann_aShwstpv",     &ann_aShwstpv,    "ann_aShwstpv/F");    
01647    tree->Branch("ann_aPh3",         &ann_aPh3,        "ann_aPh3/F");       
01648    tree->Branch("ann_aPh6",         &ann_aPh6,        "ann_aPh6/F");        
01649    tree->Branch("ann_aPhlast",      &ann_aPhlast,     "ann_aPhlast/F");     
01650    tree->Branch("ann_aPhcommon",    &ann_aPhcommon,   "ann_aPhcommon/F");   
01651 
01652   tree->Branch("day",     &day,    "day/I");
01653   tree->Branch("month",   &month,  "month/I");
01654   tree->Branch("year",    &year,   "year/I");
01655   tree->Branch("pass_beamcuts",  &pass_beamcuts, "pass_beamcuts/I");
01656 
01657 
01658   // ANDY's PID inputs
01659 
01660   tree->Branch("abid_costheta", &abid_costheta, "abid_costheta/F");
01661   tree->Branch("abid_eventenergy", &abid_eventenergy, "abid_eventenergy/F");
01662   tree->Branch("abid_trackcharge", &abid_trackcharge, "abid_trackcharge/I");
01663   tree->Branch("abid_trackenergy", &abid_trackenergy, "abid_trackenergy/i");
01664   tree->Branch("abid_tracklikeplanes", &abid_tracklikeplanes, "abid_tracklikeplanes/I");
01665   tree->Branch("abid_trackphfrac", &abid_trackphfrac, "abid_trackphfrac/F");
01666   tree->Branch("abid_trackphmean", &abid_trackphmean, "abid_trackphmean/F");
01667   tree->Branch("abid_trackqpsigmaqp", &abid_trackqpsigmaqp, "abid_trackqpsigmaqp/F");
01668   tree->Branch("abid_y", &abid_y, "abid_y/F");
01669 
01670 
01671   // Rustem's PID inputs
01672 
01673   tree->Branch("roid_numplanes", &roid_numplanes, "roid_numplanes/F");
01674   tree->Branch("roid_meansigcor", &roid_meansigcor, "roid_meansigcor/F");
01675   tree->Branch("roid_lowphhiph", &roid_lowphhiph, "roid_lowphhiph/F");
01676   tree->Branch("roid_trkphfracclose", &roid_trkphfracclose, "roid_trkphfracclose/F");
01677 
01678 
01679 // not used
01680   //  tree->Branch("roid_trkplanefrac", &roid_trkplanefrac, "roid_trkplanefrac/F");
01681   //  tree->Branch("roid_trkphfracend", &roid_trkphfracend, "roid_trkphfracend/F");
01682 
01683 
01684 
01685   tree->Branch("inu", &inu, "inu/I");
01686 
01687 
01688 
01689   // new data quality variables
01690 
01691   tree->Branch("pass_fd_qualcuts_db", &pass_fd_qualcuts_db, "pass_fd_qualcuts_db/I");
01692 
01693   tree->Branch("dataqual_cratemask", &dataqual_cratemask, "dataqual_cratemask/I");
01694   tree->Branch("dataqual_spilltype", &dataqual_spilltype, "dataqual_spilltype/I");
01695   tree->Branch("dataqual_spillerror", &dataqual_spillerror, "dataqual_spillerror/I");
01696   tree->Branch("dataqual_spillstatus", &dataqual_spillstatus, "dataqual_spillstatus/I");
01697   tree->Branch("dataqual_coldchips", &dataqual_coldchips, "dataqual_coldchips/I");
01698   tree->Branch("dataqual_coilcurr_sm1", &dataqual_coilcurr_sm1, "dataqual_coilcurr_sm1/F");
01699   tree->Branch("dataqual_coilcurr_sm2", &dataqual_coilcurr_sm2, "dataqual_coilcurr_sm2/F");
01700 
01701 
01702   tree->Branch("alec_hvstatus", &alec_hvstatus, "alec_hvstatus/I");
01703   tree->Branch("alec_coilstatus", &alec_coilstatus, "alec_coilstatus/I");
01704   tree->Branch("alec_spillstatus", &alec_spillstatus, "alec_spillstatus/I");
01705 
01706 
01707   // andy's FD fidvol
01708 
01709   tree->Branch("is_fid_ab",&is_fid_ab,"is_fid_ab/I",512000);
01710   tree->Branch("is_fid_2008",&is_fid_2008,"is_fid_2008/I",512000);
01711 
01712   // spill timing cuts
01713 
01714   tree->Branch("pass_fd_timing",&pass_fd_timing,"pass_fd_timing/I",512000);
01715 
01716 
01717   tree->Branch("nd_reclaim",&nd_reclaim,"nd_reclaim/I",512000);
01718 
01719   tree->Branch("coilok",&coilok,"coilok/I",512000);
01720   tree->Branch("coilpol",&coilpol,"coilpol/I",512000);
01721 
01722   tree->Branch("evtcom", &evtcom,"evtcom/F");
01723   tree->Branch("evtslc", &evtslc,"evtslc/F");
01724   tree->Branch("evtpur", &evtpur,"evtpur/F");
01725 
01726 
01727 
01728  //POT counting Tree
01729 
01730   TTree *potcount = new TTree("potcount","potcount");
01731   potcount->Branch("potall",&potall,"potall/F",512000);
01732   potcount->Branch("potcut",&potcut,"potcut/F",512000);
01733  
01734 
01735 
01736 
01737 
01738   Bool_t       first_ann = true;
01739 
01740   // *** NtpSt LOOP
01741 
01742   for(int i=0;i<Nentries;i++){
01743          
01744     if(i%400==0) std::cout << float(i)*100./float(Nentries) 
01745                             << "% done" << std::endl;
01746 
01747     if(!this->GetEntry(i)) continue;
01748 
01749   
01750     nevent = eventSummary->nevent;
01751     //    if(nevent==0) continue;    
01752 
01753 
01754     // header info
01755     
01756     run      = ntpHeader->GetRun();
01757     snarl    = ntpHeader->GetSnarl();
01758     subrun   = ntpHeader->GetSubRun();   
01759     trigsrc  = ntpHeader->GetTrigSrc(); 
01760     day      = eventSummary->date.day; 
01761     month    = eventSummary->date.month; 
01762     year     = eventSummary->date.year; 
01763     
01764 
01765     Double_t snarltime =ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);       
01766     // Int_t    month    = eventSummary->date.month; 
01767     litime            = eventSummary->litime; 
01768 
01769 
01770     // old data quality flag
01771     
01772 
01773     pass_fd_qualcuts   =passfail(snarltime);
01774 
01775 
01776 
01777     // NEW -> set release type for energy corrections code
01778 
01779 
01780     ReleaseType::Release_t reltype=strecord->GetRelease();
01781 
01782 
01783 
01784    // start of data quality code
01785 
01786 
01787     ntpDataQual = &(strecord->dataquality);
01788     ntpDetStatus = &(strecord->detstatus);
01789 
01790 
01791     pass_fd_qualcuts_db=1;
01792 
01793     dataqual_cratemask=0;
01794     dataqual_spilltype=0;
01795     dataqual_spillerror=0;
01796     dataqual_spillstatus=0;
01797     dataqual_coldchips=0;
01798     dataqual_coilcurr_sm1=0;
01799     dataqual_coilcurr_sm2=0;
01800 
01801     // data quality from db - following Alec.
01802 
01803     alec_hvstatus=1;
01804     alec_coilstatus=1;
01805     alec_spillstatus=1;
01806 
01807     coilok=0;
01808     coilpol=0;
01809 
01810 
01811 
01812     // NEW DATA QUALITY CHECKS - only for FD data
01813 
01814     if(ntpHeader->GetVldContext().GetSimFlag()==1 && ntpHeader->GetVldContext().GetDetector()==2) {
01815    
01816       // use singleton version of HvStatusFinder (AB 14/08/07)
01817       // HvStatusFinder hv;
01818       HvStatusFinder& hv = HvStatusFinder::Instance();
01819       HvStatus::HvStatus_t hv_ok=hv.GetHvStatus(ntpHeader->GetVldContext(),60,1);
01820 
01821       SpillServerMonFinder& smon=SpillServerMonFinder::Instance();
01822       const SpillServerMon& spill_near=smon.GetNearestSpill(ntpHeader->GetVldContext());
01823       VldTimeStamp dt=spill_near.GetSpillTime()-ntpHeader->GetVldContext().GetTimeStamp();
01824 
01825       Int_t dt_sec=abs(dt.GetSec());
01826       Int_t gps_error=spill_near.GetSpillTimeError();
01827 
01828       Bool_t coil_ok=CoilTools::IsOK(ntpHeader->GetVldContext());
01829 
01830 
01831       if(!HvStatus::Good(hv_ok)) alec_hvstatus=0;  // bad hv
01832       if(!coil_ok) alec_coilstatus=0;  // bad coil
01833       if (dt_sec<360 && gps_error>1000) alec_spillstatus=0; // bad GPS
01834 
01835   
01836       dataqual_cratemask=ntpDataQual->cratemask;
01837       dataqual_spilltype=ntpDataQual->spilltype;
01838       dataqual_spillerror=ntpDataQual->spilltimeerror;
01839       dataqual_spillstatus=ntpDataQual->spillstatus;
01840       dataqual_coldchips=ntpDataQual->coldchips;
01841       dataqual_coilcurr_sm1=ntpDetStatus->coilcurrent1;
01842       dataqual_coilcurr_sm2=ntpDetStatus->coilcurrent2;
01843 
01844       // data quality cut
01845 
01846       if (!DataUtil::IsGoodFDData(strecord)) pass_fd_qualcuts_db=0;
01847      
01848 
01849     }
01850 
01851     // end of new data quality code
01852 
01853 
01854 
01855   
01856 
01857     // Rustem's PID -  fill data for this snarl
01858     
01859     Bool_t ro_pass=false;
01860     ro_pass=ro_interface.FillSnarl(strecord);
01861     
01862 
01863 
01864     // OLD LITAG VARIABLE - no longer used
01865 
01866     litag=0; 
01867 
01868     Int_t numshower=eventSummary->nshower;
01869     Float_t allph=eventSummary->ph.raw;
01870     Int_t plbeg=eventSummary->plane.beg;
01871     Int_t plend=eventSummary->plane.end;
01872      
01873     if (numshower) {
01874      if (allph>1e6) litag+=10;
01875 
01876      if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) || 
01877          ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) || 
01878          ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) || 
01879          ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;     
01880      if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) || 
01881          ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) || 
01882          ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) || 
01883          ((plbeg==442 || plbeg==443) && (plend==484 || plend==485)))  litag++;
01884    }
01885 
01886 
01887 
01888    // NEW - check for FD LI using LISieve 
01889 
01890    lisieve=0;
01891    
01892    if (ntpHeader->GetVldContext().GetDetector()==2) {
01893    
01894      if (LISieve::IsLI(*strecord)) lisieve=1;
01895      
01896    }
01897    
01898 
01899    // coil current 
01900 
01901    magn=0;
01902 
01903 
01904 
01905 
01906     // SCAN INFO  
01907     
01908     scandecision=scan.GetDecision(run,snarl);
01909 
01910     // only pass through scanned snarls if scanfilter=1
01911     if (scanfilter==0 || scandecision>0) {  
01912 
01913       // Get Beam Monitoring Info if NOT MC: 
01914      
01915       if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01916 
01917 
01918         /*      
01919         pppinfo.GetSpillInfo(month,year,snarltime,*spinfo); 
01920                
01921         pot   =spinfo->GetPot();   
01922         horn  =spinfo->GetHorn();
01923         hhpos =spinfo->GetHpos();
01924         vvpos =spinfo->GetVpos();
01925         tar   =spinfo->GetTgtpos();
01926         magn  =spinfo->GetMagnet();   
01927         */
01928 
01929 
01930 
01931 
01932 
01933         potdb       = -999;
01934         horndb      = -999; 
01935         tardb       = -999; 
01936         vvposdb     = -999; 
01937         hhposdb     = -999; 
01938         vvwidthdb   = -999; 
01939         hhwidthdb   = -999; 
01940         timedb      = -999; 
01941         neartdb     = -999; 
01942         fartdb      = -999; 
01943         neardifftdb = -999; 
01944         
01945         // BEAM INFO DB
01946 
01947         Detector::Detector_t detectort;
01948         SimFlag::SimFlag_t simflag;
01949         VldTimeStamp  timestamp;
01950         VldContext    cx;
01951         cx        = ntpHeader->GetVldContext(); 
01952         detectort = ntpHeader->GetVldContext().GetDetector();
01953         simflag   = ntpHeader->GetVldContext().GetSimFlag();
01954         timestamp = ntpHeader->GetVldContext().GetTimeStamp();
01955         
01956         pass_beamcuts=0;
01957 
01958 
01959         // coil access tools - new, DP - Feb 1 2008
01960 
01961         Bool_t coil_ok=CoilTools::IsOK(cx);
01962         Bool_t coil_pol=CoilTools::IsReverse(cx);
01963 
01964         coilpol=1;
01965 
01966         if (coil_pol) coilpol=-1;
01967         if (coil_ok) coilok=1;
01968 
01969 
01970 
01971         const Dcs_Mag_Near* magnear =
01972           CoilTools::Instance().GetMagNear(cx);        // NearDet only
01973 
01974         if (magnear) {
01975           magn=magnear->GetCurrent();
01976         }
01977 
01978 
01979 
01980         /*
01981           DbiResultPtr<Dcs_Mag_Near> ptr(cx);
01982           if(ptr.GetNumRows()){
01983           const Dcs_Mag_Near* row = ptr.GetRow(0);
01984           magn = row->GetCurrent();
01985           }
01986         */
01987 
01988      
01989 
01990 
01991         BDSpillAccessor &spillAccess=BDSpillAccessor::Get();
01992         const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp);
01993         if(spillInfoDB) { 
01994           if(spillInfoDB->fTortgt !=0.)      potdb     = spillInfoDB->fTortgt;
01995           else if(spillInfoDB->fTrtgtd !=0.) potdb     = spillInfoDB->fTrtgtd;
01996           else if(spillInfoDB->fTor101 !=0.) potdb     = spillInfoDB->fTor101;
01997           else potdb = -999;
01998           
01999           horndb    = spillInfoDB->fHornCur;
02000           tardb     = (int)spillInfoDB->GetStatusBits().target_in;
02001           vvposdb   = spillInfoDB->fTargBpmY[0];
02002           hhposdb   = spillInfoDB->fTargBpmX[0];
02003           vvwidthdb = spillInfoDB->fProfWidY;
02004           hhwidthdb = spillInfoDB->fProfWidX;
02005           timedb    = spillInfoDB->SpillTime();
02006          
02007 
02008           // apply beam quality cuts
02009 
02010           bm.SetSpill(*spillInfoDB);
02011           bm.SetSnarlTime(timestamp);
02012 
02013           if (bm.SelectSpill()==1) pass_beamcuts=1;
02014 
02015 
02016 
02017           /*
02018           if (potdb>0.5 && fabs(snarltime-timedb)<1 && 
02019               (hhposdb<0 && hhposdb>-0.002) && 
02020               (vvposdb>0 && vvposdb<0.002) &&
02021               (hhwidthdb>0.0001 && hhwidthdb<0.0022) &&
02022               (vvwidthdb>0.0001 && vvwidthdb<0.005) &&
02023               (horndb<-155 && horndb>-200)) pass_beamcuts=1;
02024           */
02025 
02026           if (potdb>0.2 && fabs(snarltime-timedb)<1) potall=potall+potdb;
02027 
02028           //      if (pass_beamcuts && horndb<-155 && magn<-1000) potcut=potcut+potdb;
02029         
02030 
02031           // beam+magnet quality cuts for calculating POT
02032 
02033           if (pass_beamcuts && horndb<-155 && coilok==1 && coilpol==1) potcut=potcut+potdb;
02034 
02035         }
02036 
02037         
02038         
02039         // TIME DB
02040         SpillTimeFinder& spillfinder = SpillTimeFinder::Instance();
02041         fartdb                       = spillfinder.GetTimeOfNearestSpill(cx);
02042         const SpillTimeND& spnd      = spillfinder.GetNearestSpill(cx);
02043         neartdb                      = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9);
02044         neardifftdb                  = spillfinder.GetTimeToNearestSpill(cx);
02045         snarlt                       = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
02046         
02047       }
02048       
02049       //singletimeframe=SingleTimeFrame(i,Nentries);
02050       
02051      Int_t trkcount=0;
02052      Int_t shwcount=0;
02053      Float_t totph=0;
02054       
02055      Int_t evtmin=0;
02056      Int_t evtmax=nevent;
02057 
02058       if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
02059         //      cout << "MULTIPLE events!!\n";
02060 
02061         Int_t maxevent=0;
02062         Float_t maxph=0;
02063          for(int j=0;j<nevent;j++){
02064            if(!LoadEvent(j)) continue; 
02065            //  cout << "evtlength=" << ntpEvent->plane.n << " totph=" 
02066            //   << ntpEvent->ph.sigcor <<  "ntrack=" << ntpEvent->ntrack
02067            //   << " nshower=" << ntpEvent->nshower << endl;
02068            Float_t evtpheight=ntpEvent->ph.sigcor;
02069            if (evtpheight>maxph) {
02070              maxph=evtpheight;
02071              maxevent=j;
02072            }
02073          }
02074          //      cout << "BIGGEST event=" << maxevent << endl;
02075          evtmin=maxevent;
02076          evtmax=maxevent+1;
02077       }
02078 
02079 
02080 
02081       // *** EVENT LOOP - fill tree once for each reconstructed event that is the biggest in a snarl
02082 
02083       for (int j=evtmin;j<evtmax;j++){ 
02084 
02085         if(!LoadEvent(j)) continue; //no event found
02086 
02087         totph = 0;
02088 
02089         //set event index:
02090         event = j;
02091         ntrack = ntpEvent->ntrack;
02092         nshower = ntpEvent->nshower;
02093         totdigits=ntpEvent->ndigit;
02094         totstrips=ntpEvent->nstrip;
02095                 
02096         // Initialize everything
02097         //nuparent->Zero();
02098   
02099         mcflux_tpx=0; mcflux_tpy=0; mcflux_tpz=0; mcflux_tptype=0;   
02100         mcflux_nenergyN=0; mcflux_nenergyF=0; mcflux_nwtN=0; mcflux_nwtF=0;
02101         mcflux_nimpwt=0; mcflux_fluxrun=0; mcflux_fluxevtno=0; 
02102         true_enu = 0; true_emu = 0; true_eshw = 0; 
02103         true_pxnu = 0; true_pynu = 0; true_pznu = 0;
02104         true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
02105         true_dircosneu = 0; true_dircosz = 0;
02106         true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
02107         flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
02108         initial_state = 0; had_fs = 0;
02109         true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;     
02110         detector = -1; mcevent = -1;
02111         is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; 
02112         pid0 = -999; pid1 = -999; pid2 = -999; pass = 0; passcut=0;
02113         reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0; 
02114 
02115         reco_eshwcc_uncorr=0;
02116         reco_emu_uncorr=0;
02117 
02118         reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0; 
02119         reco_dircosz = 0;
02120         reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999;
02121         evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1;
02122         trkmomrange = -999; trkqp = -999; trkeqp = -999;
02123         trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.;   
02124         trkcosx=0.; trkcosy=0.; trkcosz=0.; evanglet=0.;
02125         trkendz=999; trkendu=999; trkendv=999; trkendx=999;
02126         trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999;
02127         trkds=-1; evttimemin=-1; evttimemax=-1;
02128         evtph=0.;
02129         annpid=-999.;
02130         annpid_f1p1=-999.;annpid_f1p2=-999.;annpid_f1p3=-999.;
02131         annpid_f2p1=-999.;annpid_f2p2=-999.;annpid_f2p3=-999.;
02132         reco_eshwcc=0.;reco_eshwnc=0.;reco_eshwccw=0.;reco_eshwncw=0.;
02133         is_fid_dp=0;is_fid_ns=0;
02134         is_cont_trk=0; trkmombest=-999.;
02135 
02136         trkfitpass=0;
02137         trkdiffuv=0;
02138         trkdigits=0; trkstrips=0; trkph=0; trkndof=0; trkchi2=0;
02139         shwdigits=0; shwstrips=0; shwph=0;      
02140         shwvtxu=0;
02141         shwvtxv=0;
02142         shwvtxx=0;
02143         shwvtxy=0;
02144         shwvtxz=0;
02145         shwncluster=0;
02146         numplanes=0;
02147         
02148         ann_aTotrk       =0.;
02149         ann_aTotstp      =0.;
02150         ann_aTotph       =0.;
02151         ann_aTotlen      =0.;
02152         ann_aPhperpl     =0.;
02153         ann_aPhperstp    =0.;
02154         
02155         
02156         ann_aTrkpass     =0.;
02157         ann_aTrkph       =0.;
02158         ann_aTrklen      =0.;
02159         ann_aTrkphperpl  =0.;
02160         ann_aTrkphper    =0.;
02161         ann_aTrkplu      =0.;
02162         ann_aTrkplv      =0.;
02163         ann_aTrkstp      =0.;
02164         
02165         ann_aShwph       =0.;
02166         ann_aShwstp      =0.;
02167         ann_aShwdig      =0.;
02168         ann_aShwpl       =0.;
02169         ann_aShwphper    =0.;
02170         ann_aShwphperpl  =0.;
02171         ann_aShwphperdig =0.;
02172         ann_aShwphperstp =0.;
02173         ann_aShwplu      =0.;
02174         ann_aShwplv      =0.;
02175         ann_aShwstpu     =0.;
02176         ann_aShwstpv     =0.;
02177         ann_aPh3         =0.;
02178         ann_aPh6         =0.;
02179         ann_aPhlast      =0.;
02180         ann_aPhcommon    =0.;
02181         
02182         
02183         abid_costheta         =0.;
02184         abid_eventenergy      =0.;
02185         abid_trackcharge      =0;
02186         abid_trackenergy      =0;
02187         abid_tracklikeplanes  =0;
02188         abid_trackphfrac      =0.;
02189         abid_trackphmean      =0.;
02190         abid_trackqpsigmaqp   =0.;
02191         abid_y                =0.;
02192 
02193 
02194 
02195         roid_numplanes        =0;
02196         roid_meansigcor       =0.;
02197         roid_lowphhiph        =0.;
02198         roid_trkphfracclose   =0.;
02199  
02200         //      roid_trkplanefrac     =0.;
02201         //      roid_trkphfracend     =0.;
02202 
02203   
02204 
02205         inu=0;
02206         pass_cccuts=0;
02207         pass_cccuts_2007=0;
02208         pass_cccuts_2008=0;
02209 
02210         is_fid_ab=0;
02211         is_fid_2008=0;
02212 
02213         nd_reclaim=0;
02214 
02215 
02216         evtcom=0.; evtpur=0.; evtslc=0.;
02217 
02218         // end of initialisation...
02219 
02220 
02221         rawph=ntpEvent->ph.raw;
02222 
02223         //get detector type:
02224         if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02225           detector = 2;
02226           totph=eventSummary->ph.sigcor;
02227           //fiducial criteria on vtx for far detector
02228           is_fid = 1;
02229           if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
02230         }
02231         else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02232           detector = 1;
02233         
02234           //get total charge in trk/shw
02235           for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
02236             LoadTrack(ii);
02237             totph+=ntpTrack->ph.sigcor;
02238           }
02239           trkcount+=ntrack;
02240           
02241           for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
02242             LoadShower(ii);
02243             totph+=ntpShower->ph.sigcor;
02244           }
02245           shwcount+=nshower;
02246           
02247           //fiducial criteria on vtx for near detector
02248           is_fid = 1;
02249 
02250           // ND vertex
02251           if(!PassFidND(ntpEvent->index)) is_fid = 0;
02252         }
02253 
02254 
02255         // legacy fidvol definitions...
02256 
02257         if (detector==2) is_fid_dp=PassFidNew(ntpEvent->index);
02258         else is_fid_dp=PassFidND(ntpEvent->index);
02259         
02260         if (detector==2) is_fid_ns=PassFidNewN(ntpEvent->index);
02261         else is_fid_ns=PassFidND(ntpEvent->index);
02262 
02263 
02264         if (detector==2) is_fid_ab=IsFidFD_AB(ntpEvent->index);
02265         else is_fid_ab=PassFidND(ntpEvent->index);
02266 
02267 
02268 
02269         // NEW 2008 FIDVOL DEFINITION:
02270 
02271         is_fid_2008=IsFid_2008(ntpEvent->index);
02272 
02273 
02274                 
02275         //check for a corresponding mc event 
02276     
02277         if(ntpHeader->GetVldContext().GetSimFlag()==4){ 
02278           if(LoadTruthForRecoTH(j,mcevent)) {
02279 
02280           if(LoadTHEvent(j)) { 
02281             evtcom = (ntpTHEvent->completeall);
02282             evtslc = (ntpTHEvent->completeslc);
02283             evtpur = (ntpTHEvent->purity);
02284           }
02285 
02286           Float_t *vtx = TrueVtx(mcevent);
02287           Float_t *nu3mom = TrueNuMom(mcevent);
02288           Float_t *targ4vec = Target4Vector(mcevent);
02289           //true info for tree:
02290           is_mc             = 1;
02291           nucleus           = Nucleus(mcevent);
02292           flavour           = Flavour(mcevent);
02293           initial_state     = Initial_state(mcevent);
02294           had_fs            = HadronicFinalState(mcevent);
02295           process           = IResonance(mcevent); 
02296           cc_nc             = IAction(mcevent);
02297           true_enu          = TrueNuEnergy(mcevent); 
02298           true_pxnu         = nu3mom[0];
02299           true_pynu         = nu3mom[1];
02300           true_pznu         = nu3mom[2];
02301           true_emu          = TrueMuEnergy(mcevent); 
02302           true_eshw         = TrueShwEnergy(mcevent); 
02303           true_x            = X(mcevent);
02304           true_y            = Y(mcevent);
02305           true_q2           = Q2(mcevent);
02306           true_w2           = W2(mcevent);
02307           true_dircosz      = TrueMuDCosZ(mcevent);
02308           true_dircosneu    = TrueMuDCosNeu(mcevent);
02309           true_vtxx         = vtx[0];
02310           true_vtxy         = vtx[1];
02311           true_vtxz         = vtx[2];
02312           true_etgt         = targ4vec[3];
02313           true_pxtgt        = targ4vec[0];
02314           true_pytgt        = targ4vec[1];
02315           true_pztgt        = targ4vec[2];
02316           inu               = INu(mcevent);      
02317 
02318           if (LoadTruth(mcevent)) {
02319             mcflux_tpx=ntpTruth->flux.tpx;
02320             mcflux_tpy=ntpTruth->flux.tpy;
02321             mcflux_tpz=ntpTruth->flux.tpz;
02322             mcflux_tptype=ntpTruth->flux.tptype;
02323             mcflux_nenergyN=ntpTruth->flux.nenergynear;
02324             mcflux_nenergyF=ntpTruth->flux.nenergyfar;
02325             mcflux_nwtN=ntpTruth->flux.nwtnear;
02326             mcflux_nwtF=ntpTruth->flux.nwtfar;
02327             mcflux_nimpwt=ntpTruth->flux.nimpwt;
02328             mcflux_fluxrun=ntpTruth->flux.fluxrun;
02329             mcflux_fluxevtno=ntpTruth->flux.fluxevtno;
02330           }       
02331 
02332           delete [] vtx;
02333           delete [] nu3mom;
02334           delete [] targ4vec;
02335         }
02336         }
02337         
02338         if(PassBasicCuts(event)) {pass_basic=1;} // make this reflect pass/fail flag
02339         
02340 
02341         //reco info for tree:
02342 
02343         reco_vtxx         = ntpEvent->vtx.x;
02344         reco_vtxy         = ntpEvent->vtx.y;
02345         reco_vtxz         = ntpEvent->vtx.z;
02346         evtlength         = ntpEvent->plane.end-ntpEvent->plane.beg+1;
02347         evtph             = ntpEvent->ph.sigcor; 
02348 
02349         numplanes         = ntpEvent->plane.n;
02350         
02351         // different definition for neardet - accounts for 
02352         // uninstrumented planes
02353 
02354         AnnInputBlock anninput=AnnVar(event);
02355 
02356         ann_aTotrk       =   anninput.aTotrk;    
02357         ann_aTotstp      =   anninput.aTotstp;       
02358         ann_aTotph       =   anninput.aTotph;        
02359         ann_aTotlen      =   anninput.aTotlen;       
02360         ann_aPhperpl     =   anninput.aPhperpl;      
02361         ann_aPhperstp    =   anninput.aPhperstp;     
02362         
02363         
02364         ann_aTrkpass     =   anninput.aTrkpass;      
02365         ann_aTrkph       =   anninput.aTrkph;        
02366         ann_aTrklen      =   anninput.aTrklen;       
02367         ann_aTrkphperpl  =   anninput.aTrkphperpl;   
02368         ann_aTrkphper    =   anninput.aTrkphper;     
02369         ann_aTrkplu      =   anninput.aTrkplu;       
02370         ann_aTrkplv      =   anninput.aTrkplv;       
02371         ann_aTrkstp      =   anninput.aTrkstp;       
02372         
02373         ann_aShwph       =   anninput.aShwph;        
02374         ann_aShwstp      =   anninput.aShwstp;       
02375         ann_aShwdig      =   anninput.aShwdig;       
02376         ann_aShwpl       =   anninput.aShwpl;        
02377         ann_aShwphper    =   anninput.aShwphper;     
02378         ann_aShwphperpl  =   anninput.aShwphperpl;   
02379         ann_aShwphperdig =   anninput.aShwphperdig ; 
02380         ann_aShwphperstp =   anninput.aShwphperstp;  
02381         ann_aShwplu      =   anninput.aShwplu;       
02382         ann_aShwplv      =   anninput.aShwplv;       
02383         ann_aShwstpu     =   anninput.aShwstpu;       
02384         ann_aShwstpv     =   anninput.aShwstpv; 
02385         ann_aPh3         =   anninput.aPh3;          
02386         ann_aPh6         =   anninput.aPh6;          
02387         ann_aPhlast      =   anninput.aPhlast;       
02388         ann_aPhcommon    =   anninput.aPhcommon;     
02389 
02390         //unused:       Double_t trgtime=eventSummary->trigtime;
02391         //      evttimemin = anninput.aTimemin-trgtime;
02392         //      evttimemax = anninput.aTimemax-trgtime;
02393          
02394 
02395         // cedar - use snarl time rather than trigger time
02396         evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
02397         evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
02398         
02399         // apply spill timing cut for FD data - DP NEW 27/11/05
02400 
02401         pass_fd_timing=1;
02402 
02403         if (detector==2 && ntpHeader->GetVldContext().GetSimFlag()==1) {      
02404           pass_fd_timing=PassTimingCuts(snarlt-fartdb+evttimemin);        
02405         }
02406 
02407 
02408 
02409         // ANN PID // 
02410         
02411      
02412         Int_t target=1;
02413         if(ntpHeader->GetVldContext().GetSimFlag()!=4){
02414           //if(tar<5000.  && tar>3500.)    target=1;
02415           //if(tar>5000.  && tar<40000.)   target=2;
02416           //if(tar>40000. && tar<100000.)  target=3; // Used to work just fine for the pre shutdown data
02417           // after the shutdown this variable was not properly filled from ACNET data. in fact it is not
02418           // recoverable since it is corrupted in ACNET data.
02419           if(mcneartype==1) target=1;   
02420           if(mcneartype==2) target=2; 
02421           if(mcneartype==3) target=3; 
02422         }
02423         if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){
02424           if(mcneartype==1) target=1;   
02425           if(mcneartype==2) target=2; 
02426           if(mcneartype==3) target=3; 
02427         }
02428         if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){
02429           target=0; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC 
02430         }                       
02431         // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING
02432         
02433         Int_t fid    =1;
02434         Int_t prior = 1;
02435         
02436          if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02437            if(!PassFidND(event) && !PassFid(event) && is_fid_2008==0) annpid=-999; // changed - feb 2 2008 !!!
02438           else {
02439            annpid=GetAnnPid(anninput,detector,target,fid,prior,first_ann);      
02440            first_ann=false; 
02441           }
02442          }
02443          if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02444            if(!is_fid_ns){
02445             annpid_f2p1=-999;
02446             annpid_f2p2=-999;
02447             annpid_f2p3=-999;
02448            }
02449            if(is_fid_ns){
02450             fid=2;
02451             prior=1;
02452             annpid_f2p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02453             prior=2;
02454             annpid_f2p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02455             prior=3;
02456             annpid_f2p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02457             first_ann=false;
02458            }
02459            //      if(!is_fid_dp){
02460            
02461            if(!is_fid_2008){ // changed - feb 2 2008 !!
02462              
02463              annpid_f1p1=-999;
02464              annpid_f1p2=-999;
02465              annpid_f1p3=-999;
02466            }
02467 
02468            //      if(is_fid_dp){
02469 
02470            if(is_fid_2008){ // changed - feb 2 2008 !!
02471            
02472             fid=1;
02473             prior=1;
02474             annpid_f1p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02475             prior=2;
02476             annpid_f1p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02477             prior=3;
02478             annpid_f1p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);         
02479              first_ann=false;
02480            }            
02481          }
02482          // MAK March 8, 2005
02483          // data or MC?
02484          // needed by energy corrections below
02485          bool isdata=true;
02486          if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false;
02487          
02488         //index will -1 if no track/shower present in event
02489         int track_index   = -1;
02490         LoadLargestTrackFromEvent(j,track_index);
02491         int shower_index  = -1;
02492         //LoadLargestShowerFromEvent(j,shower_index); // Change
02493         LoadShower_Jim(j,shower_index,detector);
02494         if(shower_index==-1) nshower=0;
02495         if(shower_index!=-1){ //check shower is present before using ntpShower
02496          
02497            shwdigits    = ntpShower->ndigit;
02498            shwstrips    = ntpShower->nstrip;
02499            shwph        = ntpShower->ph.sigcor;
02500            shwvtxu      = ntpShower->vtx.u;
02501            shwvtxv      = ntpShower->vtx.v;
02502            shwvtxx      = ntpShower->vtx.x;
02503            shwvtxy      = ntpShower->vtx.y;
02504            shwvtxz      = ntpShower->vtx.z;
02505            shwncluster  = ntpShower->ncluster;
02506 
02507            Int_t shwtype=0; 
02508 
02509            /*
02510            // OLD CODE
02511            // MAK March 8, 2005
02512 
02513            const Int_t shw_correction_mode=1; // Niki's re-tuning
02514            shwtype=-1;
02515            reco_eshw         = RecoShwEnergy(shower_index,shwtype,
02516                                              detector,isdata,
02517                                              shw_correction_mode);
02518            shwtype=0;
02519            reco_eshwcc       = RecoShwEnergy(shower_index,shwtype,
02520                                              detector,isdata,
02521                                              shw_correction_mode);
02522            shwtype=2;
02523            reco_eshwnc       = RecoShwEnergy(shower_index,shwtype,
02524                                              detector,isdata,
02525                                              shw_correction_mode); 
02526            shwtype=1;
02527            reco_eshwccw      = RecoShwEnergy(shower_index,shwtype,
02528                                              detector,isdata,
02529                                              shw_correction_mode); 
02530            shwtype=3;
02531            reco_eshwncw      = RecoShwEnergy(shower_index,shwtype,
02532                                              detector,isdata,
02533                                              shw_correction_mode);
02534 
02535            */
02536 
02537 
02538            // NEW CODE
02539            // DP May 19th 2007
02540 
02541            // note - release type now determined from NtpStRecord
02542          
02543            shwtype=-1;
02544            reco_eshw         = RecoShwEnergyNew(shower_index,shwtype,
02545                                                 ntpHeader->GetVldContext(),
02546                                                 reltype,corrver);
02547            shwtype=0;
02548            reco_eshwcc       = RecoShwEnergyNew(shower_index,shwtype,
02549                                                 ntpHeader->GetVldContext(),
02550                                                 reltype,corrver);
02551 
02552            // other correction factors not yet supported
02553 
02554            shwtype=2;
02555            reco_eshwnc       = ntpShower->shwph.linNCgev;
02556 
02557            shwtype=1;
02558            reco_eshwccw      = ntpShower->shwph.wtCCgev;
02559 
02560            shwtype=3;
02561            reco_eshwncw      = ntpShower->shwph.wtNCgev;
02562 
02563 
02564 
02565            // uncorrected shower energy
02566            reco_eshwcc_uncorr      = ntpShower->shwph.linCCgev;
02567 
02568 
02569 
02570 
02571            reco_eshw_sqrt    = RecoShwEnergySqrt(event);           
02572 
02573         }
02574          
02575          
02576          if(track_index!=-1){ //check track is present before using ntpTrack
02577            // MAK March 8, 2005    
02578 
02579            //     reco_emu = RecoMuEnergy(0,track_index,isdata,
02580            //                     ntpHeader->GetVldContext().GetDetector());
02581         
02582            // DP May 19th 2007 - new version for Cedar/Daikon
02583            reco_emu = RecoMuEnergyNew(0,track_index,ntpHeader->GetVldContext(),
02584                                      reltype,corrver);
02585 
02586 
02587 
02588            // muon energy without corrections
02589 
02590            const float mumass=0.10566;
02591  
02592            Float_t mr=ntpTrack->momentum.range;
02593          
02594            float mc=0.0; // signed momentum from curvature
02595            if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
02596 
02597            if(IsFidAll(track_index) || ntpTrack->momentum.qp==0) {
02598              reco_emu_uncorr=sqrt(mr*mr+ mumass*mumass);
02599            }
02600            else {
02601              if(ntpTrack->momentum.qp == 0.0) reco_emu_uncorr=10000.0; 
02602              else reco_emu_uncorr=sqrt(mc*mc+ mumass*mumass);
02603            }
02604 
02605   
02606 
02607           trklength         = ntpTrack->plane.end-ntpTrack->plane.beg+1; 
02608           trkmomrange       = ntpTrack->momentum.range;
02609           trkqp             = ntpTrack->momentum.qp;
02610           trkeqp            = ntpTrack->momentum.eqp;
02611           trkendz           = ntpTrack->end.z;
02612           trkendu           = ntpTrack->end.u;
02613           trkendv           = ntpTrack->end.v;
02614           trkendx           = ntpTrack->end.x;
02615           trkendy           = ntpTrack->end.y;
02616           
02617           trkvtxz           = ntpTrack->vtx.z;
02618           trkvtxu           = ntpTrack->vtx.u;
02619           trkvtxv           = ntpTrack->vtx.v;
02620           trkvtxx           = ntpTrack->vtx.x;
02621           trkvtxy           = ntpTrack->vtx.y;
02622           trkds             = ntpTrack->ds;             
02623           trkplbu           = ntpTrack->plane.begu;
02624           trkplbv           = ntpTrack->plane.begv;
02625           trkpleu           = ntpTrack->plane.endu;
02626           trkplev           = ntpTrack->plane.endv;
02627           trkple            = ntpTrack->plane.end;
02628           trkplb            = ntpTrack->plane.beg;
02629           trkfitpass        = ntpTrack->fit.pass;
02630           trkdiffuv         = abs(ntpTrack->plane.nu-ntpTrack->plane.nv);
02631           trkdigits         = ntpTrack->ndigit;
02632           trkstrips         = ntpTrack->nstrip;
02633           trkph             = ntpTrack->ph.sigcor;
02634           trkndof           = ntpTrack->fit.ndof;
02635           trkchi2           = ntpTrack->fit.chi2;
02636           trkcosx           = ntpTrack->vtx.dcosx;
02637           trkcosy           = ntpTrack->vtx.dcosy;
02638           trkcosz           = ntpTrack->vtx.dcosz;
02639 
02640           // account for +/-3.33 deg beam angle wrt z in near/far detectors
02641 
02642           Float_t pbeamy,pbeamz,pbeamx;
02643           
02644           pbeamy=0; pbeamz=0; pbeamx=0;
02645 
02646           if(detector==2) {
02647             pbeamx=0.;
02648             pbeamy=sin(-3.33*3.14/180.);
02649             pbeamz=cos(-3.33*3.14/180.);
02650           }
02651           
02652           if(detector==1) {
02653             pbeamx=0.;
02654             pbeamy=sin(3.33*3.14/180.);
02655             pbeamz=cos(3.33*3.14/180.);
02656           }
02657           
02658           evanglet = trkcosx*pbeamx+trkcosy*pbeamy+trkcosz*pbeamz;
02659 
02660             
02661           Float_t phtrack   =ntpTrack->ph.sigcor;
02662           Float_t phevent   =ntpEvent->ph.sigcor;
02663 
02664 
02665           // ND track reclamation cuts
02666 
02667           if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear && trkfitpass==0) {
02668 
02669             if (abs(trkplbu-trkplbv)<=5 && abs(trkpleu-trkplev)<=40 && (trkpleu<270&&trkplev<270)) nd_reclaim=1;
02670 
02671           }
02672 
02673 
02674           /* In case of Birch
02675           // apply 1.8% correction to pulse height variables
02676           if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 
02677             ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData  ){
02678             phtrack=phtrack*1.018;
02679             phevent=phevent*1.018;
02680           }
02681           */
02682 
02683           if (phevent>0) {trkphfrac=phtrack/phevent;}
02684           if(ntpTrack->plane.n>0) {
02685             trkphplane      = ntpTrack->ph.sigcor/ntpTrack->plane.n;
02686           }
02687 
02688           is_cont_trk = ntpTrack->contained;
02689           trkmombest =ntpTrack->momentum.best;
02690         }
02691         else {
02692           trklength         = 0;
02693           trkmomrange       = 0;
02694           trkqp             = 0;
02695           trkeqp            = 0;
02696           trkphfrac         = 0;
02697           trkphplane        = 0;
02698         }
02699         
02700         reco_enu          = reco_emu + reco_eshwcc;     
02701         reco_qe_enu       = RecoQENuEnergy(track_index);
02702         if (reco_enu>0) {
02703           reco_y          = reco_eshwcc/reco_enu;
02704         }
02705         reco_qe_q2        = RecoQEQ2(track_index);
02706         reco_dircosz      = RecoMuDCosZ(track_index);
02707         reco_dircosneu    = RecoMuDCosNeu(track_index);
02708         is_cev            = IsFidAll(track_index);
02709 
02710 
02711 
02712         // pre-selection:
02713 
02714         if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_dp==1) pass_cccuts=1;
02715         if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_ab==1) pass_cccuts_2007=1;
02716         if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_2008==1) pass_cccuts_2008=1;
02717 
02718 
02719 
02720 
02721       
02722         if(detector==2){
02723 
02724           Bool_t general=false;
02725           Bool_t tr1=false;
02726           Bool_t tr2=false;
02727           Bool_t tr3=false;
02728           Bool_t tr4=false;
02729           Bool_t tr5=false;
02730           Bool_t shw1=false;
02731           Bool_t shw2=false;
02732           Bool_t shw3=false;
02733           
02734           Float_t phtot  = ntpEvent->ph.sigcor;
02735           
02736           passcut=0;
02737           
02738           // substitute snarl time nsec for trgtime - due to change to definition of trgtime in Cedar
02739 
02740           if((ann_aPhperstp<2000.&& phtot>2000 && evtlength>2 && (evttimemax-evttimemin)*1e6>0.01 && fabs(litime-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9)>0.5 && nevent<=2)) general = true;
02741           // !!! changed nevent==1 to nevent < = 2
02742           if(general){
02743             if(track_index!=-1){
02744               // Double_t trkcosx= ntpTrack->vtx.dcosx;
02745               // Double_t trkcosy= ntpTrack->vtx.dcosy;
02746               // Double_t trkcosz= ntpTrack->vtx.dcosz;
02747               // Double_t trkplb = ntpTrack->plane.beg;
02748               if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.6 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>13)) tr1=true;
02749               if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.3 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>15)) tr2=true;
02750               if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.9 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>10)) tr3=true;
02751               
02752               
02753               if(trkplb>2 && trkplb!=250 && trkplb!=251 && fabs(trkcosy)<0.7 && fabs(trkvtxy-trkendy)<5.5&&fabs(trkcosz)>0.5) tr4=true; 
02754               // !!! change the fabs(trkvtxy-trkendy)<4 to fabs(trkvtxy-trkendy)<5.5
02755               if (trkcosx==0) {tr5=true;}
02756               else {
02757                 if (fabs(trkcosx)*(trkcosy/trkcosx)<0.8) {tr5=true;}
02758               }
02759               if(tr1&&tr2&&tr3&&tr4&&tr5) passcut=1;
02760               
02761               }
02762             
02763             if(track_index==-1){
02764               
02765               Double_t shwvx=ntpShower->vtx.x;
02766               Double_t shwvy=ntpShower->vtx.y;
02767               Double_t shwvz=ntpShower->vtx.z;
02768               
02769               Double_t shwcosx=ntpShower->vtx.dcosx;
02770               Double_t shwcosy=ntpShower->vtx.dcosy;
02771               Double_t shwcosz=ntpShower->vtx.dcosz;      
02772               Double_t shwplb =ntpShower->plane.beg;
02773               Double_t shwph  =ntpShower->ph.sigcor;
02774               
02775               Double_t eex=ntpEvent->end.x;
02776               Double_t eey=ntpEvent->end.y;
02777               Double_t evx=ntpEvent->vtx.x;
02778               Double_t evy=ntpEvent->vtx.y;
02779               
02780               Double_t shwstpphtemp;
02781               Double_t shwphmax=-1;
02782               Double_t shwphmin=10e9;
02783               Double_t shphpl[500];
02784             
02785               for(Int_t spl=0;spl<500;spl++) shphpl[spl]=0;
02786               
02787               Double_t shwstpu=0;
02788               Double_t shwstpv=0;
02789               Double_t shwphmaxpl;
02790               Double_t shphplmax;
02791               Double_t shphplmaxv;
02792               Double_t shwphminpl;
02793               
02794               for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
02795                 
02796                 Int_t stp_indexs=((ntpShower->stp)[jj]);
02797         
02798                 if (stp_indexs==-1) continue;  // protect against index -1!!
02799         
02800                 LoadStrip(stp_indexs);  
02801                 
02802                 if((ntpStrip->plane)%2==1) shwstpu=shwstpu+1;
02803                 if((ntpStrip->plane)%2==0) shwstpv=shwstpv+1;   
02804         
02805                 shwstpphtemp=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02806         
02807                 if(shwstpphtemp>shwphmax) {
02808                   shwphmax             =shwstpphtemp; 
02809                   shwphmaxpl           =ntpStrip->plane-ntpShower->plane.beg;
02810                 }
02811         
02812                 if(shwstpphtemp<shwphmin) {
02813                   shwphmin             =shwstpphtemp;
02814                   shwphminpl           =ntpStrip->plane-ntpShower->plane.beg;
02815                 }
02816                 
02817                 for(Int_t spl=ntpShower->plane.beg;spl<=ntpShower->plane.end;spl++){
02818                   if(spl==ntpStrip->plane){
02819                     shphpl[spl]        =shphpl[spl]+shwstpphtemp;
02820                   }
02821                 }        
02822               }    
02823               
02824               shphplmax=-1; 
02825               shphplmaxv=-1;
02826               
02827               for(Int_t spl=ntpShower->plane.beg;spl<=ntpShower->plane.end;spl++){
02828                 
02829                 if(shphpl[spl]>shphplmaxv){
02830                   shphplmax           =spl;
02831                   shphplmaxv            =shphpl[spl];
02832                 }
02833              
02834               }      
02835              
02836               shphplmax=shphplmax-ntpShower->plane.beg;      
02837               
02838               Double_t diff1=evx-eex;
02839               Double_t diff2=evy-eey;
02840               Double_t diff3=shwstpu-shwstpv;
02841               // !!
02842               if((shwvx*shwvx+shwvy*shwvy)<19. && shwvz>0.2 && shwvz<29.5 && fabs(diff1)<2. && fabs(diff2)<2. && fabs(shwvx)<3.5 &&fabs(shwvy)<3.5) shw1= true;
02843               
02844               if((shwcosx*shwcosx+shwcosy*shwcosy)<0.8 && shwcosz>0.6 && shwplb<480)       shw2=true;
02845               
02846               if(fabs(diff3)<30. && (shphplmaxv/shwph)<0.9 && (shwphmax/shwph)<0.8)        shw3=true;
02847               
02848               if(shw1&&shw2&&shw3) passcut=1;   
02849               
02850             }
02851           }
02852         }
02853         
02854 
02855 
02856 
02857 
02858 
02859 
02860 
02861         // calculate PID for all events with a track
02862         
02863         if(ntrack){
02864 
02865 
02866           // DP PID       
02867 
02868           pid0              = PID(event,0);
02869 
02870 
02871 
02872           // Andy's PID
02873 
02874           pid1=-1.;
02875 
02876           // don't apply pre-selection cuts before calculating PID
02877           //      if(andy_id.PassCuts(ntpEvent,strecord)){
02878           
02879           pid1 = andy_id.CalcPID(ntpEvent,strecord);
02880           
02881           abid_costheta         =andy_id.GetRecoCosTheta(ntpEvent,strecord);
02882           abid_eventenergy      =andy_id.GetRecoE(ntpEvent,strecord);
02883           abid_trackcharge      =andy_id.GetTrackEMCharge(ntpEvent,strecord);
02884           abid_trackenergy      =andy_id.GetTrackPlanes(ntpEvent,strecord);
02885           abid_tracklikeplanes  =andy_id.GetTrackLikePlanes(ntpEvent,strecord);
02886           abid_trackphfrac      =andy_id.GetTrackPHfrac(ntpEvent,strecord);
02887           abid_trackphmean      =andy_id.GetTrackPHmean(ntpEvent,strecord);
02888           abid_trackqpsigmaqp   =andy_id.GetTrackQPsigmaQP(ntpEvent,strecord);
02889           abid_y                =andy_id.GetRecoY(ntpEvent,strecord);
02890           
02891 
02892                   
02893             
02894           // Rustem's PID - new access method
02895           
02896           if (ro_pass) {
02897 
02898 
02899             pid2=ro_interface.GetVar("knn_pid",ntpEvent);
02900             
02901             roid_numplanes      = ro_interface.GetVar("knn_01",ntpEvent);
02902             roid_meansigcor     = ro_interface.GetVar("knn_10",ntpEvent);
02903             roid_lowphhiph      = ro_interface.GetVar("knn_20",ntpEvent);
02904             roid_trkphfracclose = ro_interface.GetVar("knn_40",ntpEvent);
02905 
02906                     
02907           }
02908             
02909            
02910                   
02911           if(PassAnalysisCuts(event)) pass = 1;
02912         }
02913         else{ 
02914           pid0            = -999.;
02915           pid1            = -999.;
02916           pid2            = -999.;
02917           pass              =  0;
02918         }
02919         cout << "\t\tFilling tree\n";
02920         tree->Fill();
02921 
02922       } // end of EVENT LOOP
02923     } //  end of SCANINFO loop 
02924   } // end of NTUPLE ENTRIES loop
02925 
02926   // delete nuparent;
02927   delete spinfo;
02928 
02929   potcount->Fill();
02930   
02931   save.cd();
02932   tree->Write();
02933   potcount->Write();
02934   save.Write();
02935   save.Close();
02936 
02937 }

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

Definition at line 523 of file MadPIDAnalysis.cxx.

References Ann_weight_vector, AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aShwplu, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, and Sigmoid().

Referenced by CreatePAN().

00524 {
00525       
00526      
00527    // DET IS THE DETECTOR 1 = NEAR 2 = FAR
00528    // TAR IS THE ENERGY   1 = LE 2 = PME 3 = PHE
00529    // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS  1 = MINE (MORE STRICT)
00530    // PRIOR IS THE FAR TRAINED METHOD    , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 
00531    //                                      3 = DM2 0.002 SINTHETA2 = 0.95
00532       
00533       
00534 //      int    inneuron  = 13;
00535       int    inneuron  = 7;
00536       int    hidneuron = 15;
00537       double var;
00538      
00539       double out[15]; // input neurons
00540       double rin[15]; // first hidden layer
00541      
00542 //      double weight1[13][15];// weights first layer
00543       double weight1[7][15];
00544       double constant1[15];   // constant first layer
00545       
00546       double weighto[15];// weights second  (output) layer
00547       double constanto[1];   // constant second (output) layer  
00548       double prob;
00549       
00550 // Initialize    
00551 
00552        for(Int_t i=0;i<hidneuron;i++) {
00553          out[i]=0.;
00554          rin[i]=0.; 
00555          constant1[i]=0.;        
00556        }
00557             
00558        for(Int_t i=0;i<inneuron;i++) {
00559         for(Int_t j=0;j<hidneuron;j++){
00560           weight1[i][j]=0.;
00561         }
00562        }
00563 
00564         constanto[0]=0.;
00565         for(Int_t i=0;i<hidneuron;i++){
00566           weighto[i]=0.;
00567         }
00568       Int_t all  = inneuron*hidneuron+hidneuron+hidneuron+1;
00569       Int_t all1 = inneuron*hidneuron+hidneuron;
00570       Int_t n1  =0;
00571       Int_t nw1 =0;
00572       Int_t no  =0;
00573       Int_t nwo =0;
00574       ifstream weightfile;
00575       
00576 //   INPUT VARIABLES
00577      if(first_ann){
00578    
00579       if(det==2){  
00580       if(prior==1){ 
00581         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");    
00582         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");        
00583        }
00584        if(prior==2){ 
00585         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");    
00586         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00587        }
00588        if(prior==3){ 
00589         if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");    
00590         if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00591        }    
00592        }
00593        
00594         if(det==1){
00595          if(tar==1)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle_cedar_daikon7v2.dat");    
00596          if(tar==2)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme_cedar_daikon7v2.dat");    
00597          if(tar==3)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe_cedar_daikon7v2.dat");    
00598        }
00599       for(Int_t k=0;k<all;k++){  
00600         weightfile >> var ;
00601         Ann_weight_vector[k]=var;
00602       } 
00603       weightfile.close();   
00604       }
00605       if(det==2){ 
00606      
00607       out[0]  =  (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.;
00608       out[1]  =  (anninput.aPhperpl/10000.); 
00609     //out[1]  =  (anninput.aTotph/50000.); 
00610       out[2]  =  (anninput.aTotlen/40);
00611       out[3]  =  (anninput.aPh3/(anninput.aTotph+1.));
00612       out[4]  =  (anninput.aPh6/(anninput.aTotph+1.));
00613       out[5]  =  (anninput.aPhlast/(anninput.aTotph+1.));
00614       out[6]  =  ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
00615       /*
00616       out[0]  = anninput.aTrkphperpl/2000.;
00617       out[1]  = anninput.aTotrk*anninput.aTrkpass;
00618       out[2]  = (anninput.aTotstp/100.);
00619       out[3]  = (anninput.aTotph/50000.);   
00620       out[4]  = (anninput.aTotlen/40);
00621       out[5]  = (anninput.aPhperpl/5000.);      
00622       out[6]  = (anninput.aPhperstp/1000.);
00623       out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.;
00624       out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
00625       */
00626       }
00627       
00628      
00629       if(det==1){ 
00630        
00631        out[0]  =  (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.;
00632        out[1]  =  (anninput.aPhperpl/10000.); 
00633     // out[1]  =  (anninput.aTotph/100000.); 
00634        out[2]  =  (anninput.aTotlen/40);
00635        out[3]  =  (anninput.aPh3/(anninput.aTotph+1.));
00636        out[4]  =  (anninput.aPh6/(anninput.aTotph+1.));
00637        out[5]  =  (anninput.aPhlast/(anninput.aTotph+1.));
00638        out[6]  =  ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00639       /*
00640        out[0]  = anninput.aTrkphperpl/5000.;
00641        out[1]  = anninput.aTotrk*anninput.aTrkpass;
00642        out[2]  = (anninput.aTotstp/100.);
00643        out[3]  = (anninput.aTotph/100000.);   
00644        out[4]  = (anninput.aTotlen/40);
00645        out[5]  = (anninput.aPhperpl/10000.);    
00646        out[6]  = (anninput.aPhperstp/2000.);  
00647        out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.;
00648        out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00649       */
00650       }
00651       /*
00652       out[7]  = (anninput.aPh3/(anninput.aTotph+1.));
00653       out[8]  = (anninput.aPh6/(anninput.aTotph+1.));
00654       out[9]  = (anninput.aPhlast/(anninput.aTotph+1.));
00655       out[12] = (anninput.aShwphperdig/800.);
00656       */
00657 // INPUT FILE
00658  
00659      
00660      for(Int_t k=0;k<all;k++){  
00661        if(k<all1){
00662          if(k%(inneuron+1)==0){
00663           nw1=0;
00664           constant1[n1]=Ann_weight_vector[k];     
00665           n1=n1+1;         
00666          }
00667          else {
00668            weight1[nw1][n1-1]=Ann_weight_vector[k];        
00669            nw1=nw1+1;      
00670          }
00671         }
00672         else if(k>=all1){
00673           if((k-all1)%(hidneuron+1)==0){
00674           nwo=0;
00675           constanto[no]=Ann_weight_vector[k];     
00676           no=no+1;        
00677          }
00678          else {
00679            weighto[nwo]=Ann_weight_vector[k];      
00680            nwo=nwo+1;      
00681          }
00682         }          
00683        }
00684        
00685 
00686 // end of reading the weights
00687    
00688        // first layer         
00689        for(Int_t i=0;i<hidneuron;i++){
00690         rin[i]=constant1[i];
00691         for(Int_t j=0;j<inneuron;j++){
00692          rin[i]=rin[i]+weight1[j][i]*out[j];
00693         }
00694        }
00695        // output                        
00696         for(Int_t i=0;i<hidneuron;i++){  
00697          out[i]=Sigmoid(rin[i]);         
00698         }               
00699         
00700         prob=constanto[0];      
00701         for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
00702         
00703                         
00704         if(anninput.aTotlen>50) prob=0.;
00705                 
00706  return prob;
00707 }

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

Definition at line 2997 of file MadPIDAnalysis.cxx.

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

02999                                               {
03000   // Mike Kordosky: 11 August 2005
03001   // Defines a fiducial volume used to pass events which fill the PDFs
03002   // Transcribed from MakeMyFile() above
03003   
03004   bool is_fid=false;
03005   
03006   if(det==Detector::kFar){
03007     
03008 
03009     is_fid = true;
03010     if(z<0.5 || z>29.4 ||   //ends
03011        (z<16.5&&z>14.5) ||  //between SMs
03012        sqrt((x*x)           //radial cut
03013             +(y*y))>3.5) is_fid = false;
03014   }
03015   else if(det==Detector::kNear){
03016     
03017     //fiducial criteria on vtx for near detector
03018     is_fid = true;
03019     if(z<1 || z>5 ||
03020        sqrt(((x-1.4828)*(x-1.4828)) +
03021             ((y-0.2384)*(y-0.2384)))>1) is_fid = false;
03022   }
03023 
03024   return is_fid;
03025   
03026 }

void MadPIDAnalysis::MakeAbIDFile ( std::string  tag  ) 

Definition at line 3029 of file MadPIDAnalysis.cxx.

References andy_id, MadBase::eventSummary, MadAbID::FillPDFs(), VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), MadQuantities::IAction(), MadQuantities::INu(), Detector::kFar, MadBase::LoadEvent(), MadBase::LoadTruthForRecoTH(), MadBase::Nentries, NtpSREventSummary::nevent, MadAbID::NormalizePDFs(), MadBase::ntpEvent, MadBase::ntpHeader, MadAbID::PassCuts(), NtpSREvent::ph, NtpSRPulseHeight::sigcor, MadBase::strecord, and MadAbID::WritePDFs().

03029                                               {
03030 
03031 
03032   // name of output PDF file
03033 
03034   std::string savename = "ABHistos_" + tag + ".root";
03035 
03036 
03037   for(int i=0;i<Nentries;i++){  // loop over SNARLS
03038 
03039     
03040     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
03041                             << "% done" << std::endl;
03042     if(!this->GetEntry(i)) continue;
03043 
03044     Int_t nevent = eventSummary->nevent;
03045     if(nevent==0) continue;
03046     
03047     
03048     Int_t evtmin=0;
03049     Int_t evtmax=nevent;
03050     
03051 
03052     // for FD - if there are multiple events per snarl, choose the one with the biggest pulse height
03053 
03054     if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
03055       
03056       Int_t maxevent=0;
03057       Float_t maxph=0;
03058       for(int j=0;j<nevent;j++){
03059         if(!LoadEvent(j)) continue; 
03060         Float_t evtpheight=ntpEvent->ph.sigcor;
03061         if (evtpheight>maxph) {
03062           maxph=evtpheight;
03063           maxevent=j;
03064         }
03065       }
03066       evtmin=maxevent;
03067       evtmax=maxevent+1;
03068     }
03069     
03070     // EVENT LOOP 
03071 
03072     for(int j=evtmin;j<evtmax;j++){       
03073 
03074 
03075       if(!LoadEvent(j)) continue; //no event found
03076       
03077       
03078       Int_t ccnc = -1;
03079   
03080       //check for a corresponding mc event
03081  
03082       Int_t mcevent=-1;
03083 
03084       if(LoadTruthForRecoTH(j,mcevent)) {
03085         //truth info for event:
03086         
03087         Int_t iact             = IAction(mcevent);
03088         Int_t inu              = INu(mcevent);
03089         // unused:      Int_t ires             = IResonance(mcevent);  
03090         
03091         if(abs(inu)==14 && iact==1) ccnc=1;
03092         else ccnc=0;
03093    
03094       }
03095         
03096       // fill PDFs for this event
03097       if (andy_id.PassCuts(ntpEvent,strecord)) {
03098         andy_id.FillPDFs(ntpEvent,strecord,ccnc);
03099       }
03100       
03101     } // end of EVENT loop
03102 
03103   } // end of SNARL loop
03104 
03105 
03106   // finished looping through snarls, now normalise PDFs and write them out
03107 
03108   andy_id.NormalizePDFs();
03109   andy_id.WritePDFs(savename.c_str());
03110 
03111 
03112 }

void MadPIDAnalysis::MakeMyFile ( std::string  tag,
int  tarpos,
int  pdfbinning 
)

Definition at line 789 of file MadPIDAnalysis.cxx.

References NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, MadBase::eventSummary, NtpSRTrack::fit, MadQuantities::Flavour(), VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), MadQuantities::IAction(), NtpSREvent::index, MadQuantities::IsFid_2008(), 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, and MadQuantities::TrueNuEnergy().

00789                                                                        {
00790 
00791 
00792 
00793   // pdf binning: 0 - old binning, max 50 bins, 1 - new finer binning
00794 
00795   std::string savename = "DPHistos_" + tag + ".root";
00796   TFile save(savename.c_str(),"RECREATE"); 
00797   save.cd();
00798   
00799  
00800 
00801   //#declare lots of histos:
00802 
00803   Int_t evlbins=0;
00804 
00805   if (pdfbinning==0) evlbins=60;
00806   else evlbins=300;
00807 
00808 
00809   TH1F *evlength_nc = new TH1F("Event length - nc",
00810                                "Event length - nc",
00811                                evlbins,0,600);
00812   evlength_nc->SetXTitle("Event length (planes)");
00813   TH1F *evlength_cc = new TH1F("Event length - cc",
00814                                "Event length - cc",
00815                                evlbins,0,600);
00816   evlength_cc->SetXTitle("Event length (planes)");
00817 
00818 
00819 
00820   TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
00821                                "Track ph frac - nc",
00822                                50,0,1);
00823   phfrac_nc->SetXTitle("pulse height fraction");
00824 
00825 
00826   TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
00827                                "Track ph frac - cc",
00828                                50,0,1);
00829   phfrac_cc->SetXTitle("pulse height fraction");
00830 
00831 
00832 
00833   Int_t phplanebins=0;
00834   Float_t phplanemax=0;
00835  
00836   if (pdfbinning==0) {
00837     phplanebins=50;
00838     phplanemax=5000;
00839   }
00840   else {
00841     phplanebins=100;
00842     phplanemax=3000;
00843   }
00844 
00845   TH1F *phplane_nc = new TH1F("ph per plane (nc)",
00846                               "ph per plane (nc)",
00847                                phplanebins,0,phplanemax);
00848   phplane_nc->SetXTitle("pulse height per plane");
00849   TH1F *phplane_cc = new TH1F("ph per plane (cc)",
00850                               "ph per plane (cc)",
00851                                phplanebins,0,phplanemax);
00852   phplane_cc->SetXTitle("pulse height per plane");
00853 
00854 
00855   evlength_nc->SetDirectory(&save);
00856   evlength_cc->SetDirectory(&save);
00857 
00858   phfrac_nc->SetDirectory(&save);
00859   phfrac_cc->SetDirectory(&save);
00860 
00861   phplane_nc->SetDirectory(&save);
00862   phplane_cc->SetDirectory(&save);
00863 
00864 
00865 
00866 
00867 
00868   //  gDirectory->ls();
00869 
00870   for(int i=0;i<Nentries;i++){
00871     
00872     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00873                             << "% done" << std::endl;
00874     if(!this->GetEntry(i)) continue;
00875 
00876     Int_t nevent = eventSummary->nevent;
00877     if(nevent==0) continue;
00878     
00879     Int_t trkcount=0;
00880     Int_t shwcount=0;
00881 
00882     Int_t evtmin=0;
00883     Int_t evtmax=nevent;
00884 
00885     if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
00886       
00887       Int_t maxevent=0;
00888       Float_t maxph=0;
00889       for(int j=0;j<nevent;j++){
00890         if(!LoadEvent(j)) continue; 
00891         Float_t evtpheight=ntpEvent->ph.sigcor;
00892         if (evtpheight>maxph) {
00893           maxph=evtpheight;
00894           maxevent=j;
00895         }
00896       }
00897       evtmin=maxevent;
00898       evtmax=maxevent+1;
00899     }
00900 
00901     // EVENT LOOP - fill tree once for each reconstructed event
00902 
00903     for(int j=evtmin;j<evtmax;j++){       
00904 
00905 
00906       if(!LoadEvent(j)) continue; //no event found
00907       
00908       Int_t ntrack = 0; 
00909       ntrack=ntpEvent->ntrack;
00910       Int_t nshower = 0;
00911       nshower=ntpEvent->nshower;
00912 
00913       Float_t totph=0;
00914 
00915 
00916       for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00917         LoadTrack(ii);
00918 
00919         totph+=ntpTrack->ph.sigcor;
00920         LoadTHTrack(ii);
00921 
00922       }
00923 
00924       trkcount+=ntrack;
00925 
00926 
00927        for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00928         LoadShower(ii);
00929 
00930         totph+=ntpShower->ph.sigcor;
00931         }
00932 
00933        shwcount+=nshower;
00934  
00935 
00936 
00937       Int_t cc_nc = 0;
00938       Int_t flavour = 0;
00939       Int_t detector = -1;
00940       Int_t is_fid = 0;
00941       Double_t neu_e  = -1; 
00942       Double_t weight = -1; 
00943  
00944       // new definition of is_fid
00945 
00946       is_fid = IsFid_2008(ntpEvent->index);
00947 
00948      
00949       //get detector type:
00950 
00951       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00952         detector = 2;
00953 
00954         if (evlength_cc->GetNbinsX()==evlbins) {
00955           if (pdfbinning==0) {
00956             evlength_cc->SetBins(50,0,500);
00957             evlength_nc->SetBins(50,0,500);
00958           }
00959           else {
00960             evlength_cc->SetBins(250,0,500);
00961             evlength_nc->SetBins(250,0,500);
00962           }
00963         }
00964 
00965       }
00966       else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00967         detector = 1;   
00968 
00969         if (pdfbinning==1 && evlength_cc->GetBinLowEdge(250)>200) {
00970           evlength_cc->SetBins(150,0,300);
00971           evlength_nc->SetBins(150,0,300);
00972         }
00973         if (pdfbinning==0 && evlength_cc->GetBinLowEdge(50)>200) {
00974           evlength_cc->SetBins(60,0,300);
00975           evlength_nc->SetBins(60,0,300);
00976         }
00977 
00978 
00979       }
00980       
00981       //check for a corresponding mc event
00982  
00983       Int_t mcevent=-1;
00984 
00985       if(LoadTruthForRecoTH(j,mcevent)) {
00986         //true info for tree:
00987 
00988         cc_nc             = IAction(mcevent);
00989         flavour           = Flavour(mcevent);
00990         neu_e             = TrueNuEnergy(mcevent);  
00991                           
00992         // tarpos   1 = LE 2 = ME 3 = HE     
00993         if(tarpos==2)  weight=1; 
00994         if(tarpos==1)  weight=1;        
00995         if(tarpos==3)  weight=1; 
00996       }
00997 
00998       int track_index   = -1;
00999       LoadLargestTrackFromEvent(j,track_index);
01000 
01001       if (track_index!=-1) {
01002 
01003 
01004         Int_t trkpass=ntpTrack->fit.pass;
01005         Int_t nd_reclaim=0;
01006 
01007         Int_t trkplbu    = ntpTrack->plane.begu;
01008         Int_t trkplbv    = ntpTrack->plane.begv;
01009         Int_t trkpleu    = ntpTrack->plane.endu;
01010         Int_t trkplev    = ntpTrack->plane.endv;
01011 
01012 
01013         if (abs(trkplbu-trkplbv)<=5 && abs(trkpleu-trkplev)<=40 && (trkpleu<270&&trkplev<270)) nd_reclaim=1;
01014 
01015 
01016          // track reclamation - ND only
01017         if (detector==1 && trkpass==0 && nd_reclaim==1) trkpass=1;
01018  
01019 
01020         if (is_fid && trkpass) {
01021         
01022 
01023         Float_t evlength=0;
01024         Float_t phfrac=0;
01025         Float_t phplane=0;
01026 
01027         evlength=ntpEvent->plane.n;
01028         if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01029           evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
01030         }
01031 
01032         Float_t phtrack=ntpTrack->ph.sigcor;
01033         Float_t phevent=ntpEvent->ph.sigcor;
01034         if (phevent>0) {phfrac=phtrack/phevent;}
01035 
01036         Float_t trkplane=ntpTrack->plane.n;
01037         if (trkplane>0) {phplane=phtrack/trkplane;}
01038        
01039         if(flavour==2 && cc_nc==1) {
01040           evlength_cc->Fill(evlength,weight);
01041           phfrac_cc->Fill(phfrac,weight);
01042           phplane_cc->Fill(phplane,weight);
01043         }
01044 
01045         else if(cc_nc==0) {
01046           evlength_nc->Fill(evlength,weight);
01047           phfrac_nc->Fill(phfrac,weight);
01048           phplane_nc->Fill(phplane,weight);      
01049         }
01050 
01051         
01052         }
01053       }
01054     }
01055   }
01056   //#close file  
01057   //  std::cout<<"writing hist "<<endl;
01058   //  save.ls();
01059   save.Write();
01060   save.Close();
01061 }

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

Reimplemented from MadAnalysis.

Definition at line 120 of file MadPIDAnalysis.cxx.

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

Referenced by CreatePAN().

00121 {
00122   if(!LoadEvent(event)) return false;
00123   //  if(PassBasicCuts(event)) return true;
00124   Float_t pidcut=-0.4;
00125   if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
00126   
00127   if(!(PID(event,0)>pidcut)) return false; 
00128   return true;
00129 }

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

Reimplemented from MadAnalysis.

Definition at line 188 of file MadPIDAnalysis.cxx.

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

Referenced by CreatePAN().

00189 { 
00190   if(!LoadEvent(event)) return false;
00191    if(ntpEvent->ntrack==0) return false;
00192    //   Int_t track = -1;
00193    // LoadLargestTrackFromEvent(event,track);  
00194    // if(!PassTrackCuts(track)) return false;  
00195    //   if(!IsFidVtxEvt(ntpEvent->index)) return false;  
00196 //  if(!IsFidVtx(track)) return false;
00197   return true;
00198 }

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

Definition at line 199 of file MadPIDAnalysis.cxx.

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

Referenced by CreatePAN().

00200 { 
00201   if(!LoadEvent(event)) return false;
00202   if(!IsFidVtxEvt(ntpEvent->index)) return false;  
00203   return true;
00204 }

Bool_t MadPIDAnalysis::PassFidND ( Int_t  event = 0  )  [protected]

Definition at line 229 of file MadPIDAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by CreatePAN().

00230 { 
00231   if(!LoadEvent(event)) return false;
00232  
00233   // no track - return false...
00234 
00235   if(ntpEvent->ntrack==0) return false;
00236 
00237   Int_t track = -1;
00238   LoadLargestTrackFromEvent(event,track);
00239 
00240   // use track vertex
00241   
00242   if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>5.0 ||   //ends
00243                  
00244                     (pow(ntpTrack->vtx.x-1.4828,2)+           //radial cut
00245                     pow(ntpTrack->vtx.y-0.2384,2))>(1.0))) return false;  
00246       
00247     //      (pow(ntpTrack->vtx.x-1.4885,2)+           //radial cut
00248          //          pow(ntpTrack->vtx.y-0.1397,2))>(1.0))) return false;  
00249 
00250   return true; 
00251 
00252 }

Bool_t MadPIDAnalysis::PassFidNew ( Int_t  event = 0  )  [protected]

Definition at line 132 of file MadPIDAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSREvent::vtx, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by CreatePAN().

00133 { 
00134 
00135   if(!LoadEvent(event)) return false;
00136 
00137   Int_t track = -1;
00138   LoadLargestTrackFromEvent(event,track);
00139 
00140   // for events with no track, use event vertex
00141 
00142   if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 ||   //ends
00143            (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) ||  //between SMs
00144            (pow(ntpEvent->vtx.x,2)+           //radial cut
00145             pow(ntpEvent->vtx.y,2))>14)) return false;
00146 
00147   // otherwise, use track vertex
00148 
00149   if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 ||   //ends
00150            (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) ||  //between SMs
00151            (pow(ntpTrack->vtx.x,2)+           //radial cut
00152             pow(ntpTrack->vtx.y,2))>14)) return false;  
00153 
00154   return true; 
00155 
00156 }

Bool_t MadPIDAnalysis::PassFidNewN ( Int_t  event = 0  )  [protected]

Definition at line 160 of file MadPIDAnalysis.cxx.

References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSREvent::vtx, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by CreatePAN().

00161 { 
00162 
00163 
00164   if(!LoadEvent(event)) return false;
00165 
00166   Int_t track = -1;
00167   LoadLargestTrackFromEvent(event,track);
00168 
00169   // for events with no track, use event vertex
00170 
00171   if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 ||   //ends
00172            (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) ||  //between SMs
00173            (pow(ntpEvent->vtx.x,2)+           //radial cut
00174             pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false;
00175 
00176   // otherwise, use track vertex
00177 
00178   if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 ||   //ends
00179            (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) ||  //between SMs
00180            (pow(ntpTrack->vtx.x,2)+           //radial cut
00181             pow(ntpTrack->vtx.y,2))>(3.5*3.5))) return false;  
00182 
00183   return true;  
00184 
00185 }

Bool_t MadPIDAnalysis::PassTimingCuts ( Double_t  timediff  )  [protected]

Definition at line 214 of file MadPIDAnalysis.cxx.

References MuELoss::e.

Referenced by CreatePAN().

00215 { 
00216  
00217  
00218   if (timediff<-20e-6) return false;
00219   if (timediff>30e-6) return false;
00220 
00221 
00222   return true;
00223  
00224 }

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

Reimplemented from MadAnalysis.

Definition at line 112 of file MadPIDAnalysis.cxx.

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

00113 {
00114   if(!LoadTruth(mcevent)) return false;
00115   if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
00116   return false;
00117 }

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

Reimplemented from MadAnalysis.

Definition at line 711 of file MadPIDAnalysis.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().

00712 {
00713 
00714 
00715   if(!LoadEvent(event)) return -10;
00716   Int_t track = -1;
00717   if(!LoadLargestTrackFromEvent(event,track)) return -10;
00718   if(method!=0) return -10;
00719   
00720   Float_t ProbNC = 1.;
00721   Float_t ProbCC = 1.;
00722   Float_t PidCC = 0.;
00723   
00724 
00725   Float_t var1=eventSummary->plane.n;
00726   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00727     var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00728   }
00729   Float_t var2=0;
00730   Float_t var3=0;
00731 
00732 
00733   Float_t phtrack=ntpTrack->ph.sigcor;
00734   Float_t phevent=ntpEvent->ph.sigcor;
00735  
00736   /* In case of Birch
00737   // apply 1.8% correction to pulse height variables
00738   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 
00739      ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData  ){
00740     phtrack=phtrack*1.018;
00741     phevent=phevent*1.018;
00742   }
00743   */
00744 
00745   if (phevent>0) {var2=phtrack/phevent;}
00746 
00747   if (ntpTrack->plane.n!=0) var3 = float(phtrack)
00748                               /float(ntpTrack->plane.n);
00749 
00750 
00751   Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00752   Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00753   Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00754 
00755   if (prob1==0) {prob1=0.0001;}
00756   if (prob2==0) {prob2=0.0001;}
00757   if (prob3==0) {prob3=0.0001;}
00758  
00759   ProbCC=prob1*prob2*prob3;
00760 
00761   prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00762   prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00763   prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00764  
00765   if (prob1==0) {prob1=0.0001;}
00766   if (prob2==0) {prob2=0.0001;}
00767   if (prob3==0) {prob3=0.0001;}
00768  
00769   ProbNC=prob1*prob2*prob3;
00770 
00771   PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00772 
00773   return PidCC;
00774 }

void MadPIDAnalysis::ReadAbPIDFile ( std::string  tag  ) 

Definition at line 2966 of file MadPIDAnalysis.cxx.

References andy_id, and MadAbID::ReadPDFs().

02966                                                {
02967 
02968 
02969   cout << "Reading AbID pdfs from " << tag.c_str() << endl;
02970   andy_id.ReadPDFs(tag.c_str());
02971 
02972 }

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

Definition at line 2941 of file MadPIDAnalysis.cxx.

References fLikeliFile, and fLikeliHist.

02941                                              {
02942 
02943   fLikeliFile = new TFile(tag.c_str(),"READ");
02944  
02945   if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found
02946  
02947     fLikeliFile->cd();
02948 
02949     fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
02950     fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
02951     fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
02952     fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
02953     fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
02954     fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
02955     
02956     for(int i=0;i<6;i++){
02957       float integ = fLikeliHist[i]->Integral(); 
02958       fLikeliHist[i]->Scale(1./integ);
02959     }
02960   }
02961   else fLikeliFile = NULL;
02962 }

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

Reimplemented from MadAnalysis.

Definition at line 777 of file MadPIDAnalysis.cxx.

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

00777                                                      {
00778 
00779   if(!LoadEvent(event)) return 0;
00780   int largestTrack = -1;
00781   LoadLargestTrackFromEvent(event,largestTrack);
00782   float emu = RecoMuEnergy(0,largestTrack);
00783   float ehad = RecoShwEnergy(event);
00784   float ereco = emu + ehad; 
00785   return ereco;
00786 
00787 }

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

Definition at line 506 of file MadPIDAnalysis.cxx.

Referenced by GetAnnPid().

00507 {
00508   double sig;
00509       if(x>37.){
00510          sig = 1.;
00511       }
00512       else if(x<-37.){
00513          sig = 0.;
00514       }
00515       else {
00516          sig = 1./(1.+exp(-x));
00517       }
00518       return sig;
00519 }

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

Definition at line 1064 of file MadPIDAnalysis.cxx.

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

01064                                                                       {
01065 
01066   Int_t tf=0,tfbef=0,tfaft=0;
01067    
01068   
01069   if(this->GetEntry(snarlentry))                                tf      = ntpHeader->GetTimeFrame();
01070   if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1))   tfaft   = ntpHeader->GetTimeFrame();   
01071   if(snarlentry>0        && this->GetEntry(snarlentry-1))       tfbef   = ntpHeader->GetTimeFrame();
01072   
01073   Float_t singletimeframe=0;
01074   if(snarlentry<(nentries-1) && snarlentry>0 ){
01075     if(tf!=tfaft && tf!=tfbef) singletimeframe=1;
01076   }
01077   if(snarlentry==(nentries-1)) {
01078     if(tf!=tfbef) singletimeframe=1; 
01079   }  
01080   if(snarlentry==0) {
01081     if(tf!=tfaft) singletimeframe=1;
01082   }
01083   
01084   this->GetEntry(snarlentry);
01085   return singletimeframe;  
01086   
01087 }


Member Data Documentation

Definition at line 39 of file MadPIDAnalysis.h.

Referenced by CreatePAN(), MakeAbIDFile(), and ReadAbPIDFile().

Double_t MadPIDAnalysis::Ann_weight_vector[226] [private]

Definition at line 46 of file MadPIDAnalysis.h.

Referenced by GetAnnPid().

TFile* MadPIDAnalysis::fLikeliFile [protected]

Definition at line 36 of file MadPIDAnalysis.h.

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

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

Definition at line 37 of file MadPIDAnalysis.h.

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

Float_t MadPIDAnalysis::potall [private]

Definition at line 47 of file MadPIDAnalysis.h.

Referenced by CreatePAN(), MadPIDAnalysis(), and ~MadPIDAnalysis().

Float_t MadPIDAnalysis::potcut [private]

Definition at line 48 of file MadPIDAnalysis.h.

Referenced by CreatePAN(), MadPIDAnalysis(), and ~MadPIDAnalysis().

Definition at line 41 of file MadPIDAnalysis.h.

Referenced by ConfigureRoID(), and CreatePAN().

Definition at line 42 of file MadPIDAnalysis.h.

Referenced by ConfigureRoID().


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

Generated on 13 Sep 2017 for loon by  doxygen 1.6.1