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 ,
string  ,
int   
)

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, NtpSRShower::nstrip, NtpSRTrack::nstrip, NtpSREvent::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::plane, NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, and NtpSRShower::stp.

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  ,
std::string   
)

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  ,
Int_t  ,
Int_t  ,
EnergyCorrections::WhichCorrection_t   
)

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, NtpSREvent::end, NtpSRTrack::end, NtpSRPlane::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(), 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(), SpillTimeFinder::Instance(), CoilTools::Instance(), SpillServerMonFinder::Instance(), HvStatusFinder::Instance(), MadQuantities::INu(), MadQuantities::IResonance(), MadQuantities::IsFid_2008(), MadQuantities::IsFidAll(), MadQuantities::IsFidFD_AB(), MadQuantities::IsFidVtxEvt(), DataUtil::IsGoodFDData(), LISieve::IsLI(), CoilTools::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, NtpSRDate::month, month, NtpSRPlane::n, NtpSRShower::ncluster, NtpSRTrack::ndigit, NtpSRShower::ndigit, NtpSREvent::ndigit, NtpSRFitTrack::ndof, NtpMCFluxInfo::nenergyfar, NtpMCFluxInfo::nenergynear, MadBase::Nentries, NtpSREventSummary::nevent, NtpMCFluxInfo::nimpwt, NtpSREvent::nshower, NtpSREventSummary::nshower, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSREvent::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(), NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSREventSummary::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, PID(), NtpSRStrip::plane, NtpSRShower::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, NtpSRTrack::vtx, NtpSRShower::vtx, NtpSREvent::vtx, MadQuantities::W2(), NtpSRShowerPulseHeight::wtCCgev, NtpSRShowerPulseHeight::wtNCgev, NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::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   ) 

Definition at line 3029 of file MadPIDAnalysis.cxx.

References andy_id, MadBase::eventSummary, MadAbID::FillPDFs(), VldContext::GetDetector(), 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  ,
int  ,
int   
)

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(), 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, NtpSREvent::plane, NtpSRTrack::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, NtpSRTrack::vtx, NtpSREvent::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, NtpSRTrack::vtx, NtpSREvent::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, NtpSREvent::ph, NtpSRTrack::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   ) 

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   ) 

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

MadAbID MadPIDAnalysis::andy_id [protected]

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().

Anp::Interface MadPIDAnalysis::ro_interface [protected]

Definition at line 41 of file MadPIDAnalysis.h.

Referenced by ConfigureRoID(), and CreatePAN().

Registry MadPIDAnalysis::ro_registry [protected]

Definition at line 42 of file MadPIDAnalysis.h.

Referenced by ConfigureRoID().


The documentation for this class was generated from the following files:
Generated on Fri Oct 10 22:45:49 2014 for loon by  doxygen 1.4.7