loon: NuDemoModule.cxx Source File

00001 #include "TClonesArray.h"
00002 
00003 #include "CandNtupleSR/NtpSRTrack.h"
00004 #include "CandNtupleSR/NtpSREvent.h"
00005 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE
00006 #include "JobControl/JobCResult.h"
00007 #include "MessageService/MsgService.h"
00008 #include "MinosObjectMap/MomNavigator.h"
00009 #include "StandardNtuple/NtpStRecord.h"
00010 #include "Conventions/ReleaseType.h"
00011 
00012 #include "NtupleUtils/LISieve.h"
00013 #include "NtupleUtils/NuEvent.h"
00014 #include "NtupleUtils/NuLibrary.h"
00015 
00016 #include "NtupleUtils/NuDemoModule.h"
00017 
00018 ClassImp(NuDemoModule)
00019 
00020 CVSID("$Id: NuDemoModule.cxx,v 1.8 2007/11/24 11:26:47 hartnell Exp $");
00021 
00022 JOBMODULE(NuDemoModule,"NuDemoModule","Demo Module");
00023 
00024 //____________________________________________________________________72
00025 
00026 NuDemoModule::NuDemoModule()
00027 {
00028 
00029 }
00030 
00031 //____________________________________________________________________72
00032 
00033 NuDemoModule::~NuDemoModule()
00034 {
00035 
00036 }
00037 
00038 //____________________________________________________________________72
00039 
00040 JobCResult NuDemoModule::Ana(const MomNavigator *mom)
00041 { 
00042   //get the record from mom
00043   NtpStRecord* pntp = dynamic_cast<NtpStRecord*>
00044     (mom->GetFragment("NtpStRecord"));
00045   const NtpStRecord& ntp = *pntp;
00046 
00047   const RecCandHeader& rec=ntp.GetHeader();
00048   MAXMSG("NuAnalysis",Msg::kInfo,5)
00049     <<"Found: run="<<rec.GetRun()
00050     <<", subrun="<<rec.GetSubRun()
00051     <<", detector="<<rec.GetVldContext().GetDetector()
00052     <<", simFlag="<<rec.GetVldContext().GetSimFlag()
00053     <<endl;
00054   
00055   //get an instance of the code library
00056   NuLibrary& lib=NuLibrary::Instance();
00057 
00058   static Int_t entry=-1;
00059   entry++;//entry is now zero on first iteration
00060   
00061   TClonesArray& evtTca=(*ntp.evt);
00062   const Int_t numEvts=evtTca.GetEntriesFast();
00063   
00065   //loop over the reconstructed events
00067   for (Int_t ievt=0;ievt<numEvts;ievt++){
00068     const NtpSREvent& evt=
00069       *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
00070     lib.cnt.evtCounter++;
00071     
00072     MAXMSG("NuAnalysis",Msg::kInfo,5)
00073       <<"Found event, "<<ievt+1<<"/"<<numEvts<<endl;
00074     
00075     const NtpSREventSummary& evthdr=ntp.evthdr;
00076     
00077     if (evthdr.litime!=-1) continue;
00078     lib.cnt.evtNotlitime++;
00079     
00080     //simple event object to hold important quantities
00081     //NuEvent& nu=output.GetNuEventToFill();
00082     NuEvent nu;
00083     //set the analysis (cuts and reco) version to use
00084     nu.entry=entry;//store the entry in the tree number
00085     nu.nevt=numEvts;
00086     nu.anaVersion=NuCuts::kCC0250Std;
00087     nu.releaseType=ReleaseType::kCedarDaikon;
00088     //this->CopyConfig(config,nu);
00089   
00090     //get the run, snarl, etc info
00091     lib.ext.ExtractGeneralInfo(ntp,nu);
00092     
00093     //get info from the evt
00094     lib.ext.ExtractEvtInfo(evt,nu);
00095 
00096     //ensure good number of tracks in the event
00097     if (!lib.cuts.IsGoodNumberOfTracks(nu)) continue;
00098     lib.cnt.evtWithTrkCounter++;
00099     
00100     //get the best track in the event
00101     Int_t bestTrack=lib.reco.GetBestTrack(nu);
00102     const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,
00103                                                    bestTrack-1);
00104     if (ptrk==0) {
00105       MSG("NuAnalysis",Msg::kWarning)
00106         <<endl<<endl<<"Ahhh, no track in this event!"<<endl<<endl;
00107     }//let the segv happen...
00108     const NtpSRTrack& trk=*ptrk;
00109     
00110     //extract trk/shw info (put this after ExtractEvtInfo to save cpu)
00111     lib.ext.ExtractTrkShwInfo(ntp,evt,nu);
00112     
00114     //RECONSTRUCT the neutrino energy
00115     lib.reco.GetEvtEnergy(nu);
00116     
00117     //check if the trk is in the fiducial volume
00118     if (!lib.cuts.IsInFidVolTrk(nu)) continue;
00119     lib.cnt.trkInFidVolCounter++;
00120     
00121     //get the timing info
00122     lib.ext.ExtractMinMaxEvtTimes(ntp,evt,nu);
00123     lib.ext.ExtractTimeToNearestSpill(nu);
00124     //get the truth info (reco quantities are done above)
00125     lib.reco.GetTruthInfo(ntp,evt,nu);
00126     //do the beam reweighting
00127     //zBeamReweight.ExtractZBeamReweight(ntp,evt,nu);
00128     //lib.reco.ApplyReweights(nu);
00129     //get the data quality
00130     lib.ext.ExtractDataQuality(ntp,nu);
00131     lib.ext.ExtractLITags(ntp,nu);
00132     lib.ext.ExtractBeamInfoDB(nu);
00133     
00134     //calculate PIDs
00135     //pid.GetDpID(ntp,evt,nu);
00136     //pid.GetAbID(ntp,evt,nu);
00137     //pid.GetRoID(ntp,evt,nu);//too slow
00138     //pid.GetJeID(ntp,evt,nu);//not yet implemented
00139     //pid.GetPoID(nu);
00140     
00141     //get RO's numubar pid (0 or 1)
00142     //pid.GetRoIDNuMuBar(ntp,evt,nu);//slow?
00143     
00144     //MAXMSG("NuAnalysis",Msg::kInfo,10)
00145     //<<"PID: dp="<<nu.dpID<<", ab="<<nu.abID<<", ro="<<nu.roID
00146     //<<", po="<<nu.poID<<", roNMB="<<nu.roIDNuMuBar<<endl;
00147     
00149     //fill the tree
00151     //output.FillTree();
00152       
00153     //cut on LI
00154     if (lib.cuts.IsLI(nu)) continue;
00155     lib.cnt.evtNotIsLI++;
00156     
00157     //cut on the data quality
00158     if (!lib.cuts.IsGoodDataQuality(nu)) continue;
00159     lib.cnt.goodDataQualityCounter++;
00160       
00161     //cut on the spill time
00162     //plots.FillEvtAndSpillTimingPlots(nu);
00163     if (!lib.cuts.IsGoodTimeToNearestSpill(nu)) continue;
00164     lib.cnt.goodTimeToNearestSpillCounter++;
00165     
00166     //cut on the beam
00167     if (!lib.cuts.IsGoodBeam(nu)) continue;
00168     lib.cnt.goodBeamInfoDBCounter++;
00169     
00170     //require a good trk fit
00171     if (!lib.cuts.IsGoodTrackFitPass(nu)) continue;
00172     lib.cnt.goodTrkPassCounter++;
00173     
00174     //require a forward going neutrino about beam direction
00175     if (!lib.cuts.IsGoodDirCos(nu)) continue;
00176     lib.cnt.goodDirectionCosineCounter++;
00177     
00178     //fill pid plots AFTER pre-selection
00179     //plots.FillTrueDpIDHistos(nu);
00180     //plots.FillTrueDpIDHistosPQNQ(nu);
00181     //plots.FillDPIdSigmaQPPassPreSelCutPlots(nu);
00182     //plots.FillSigmaQPPlots(nu);//sig, chi, prob plots
00183     //plots.FillTruePoIDHistos(nu);
00184     //plots.FillTruePoIDHistosPQNQ(nu);
00185     //plots.FillTrueAbIDHistos(nu);
00186     //plots.FillTrueAbIDHistosPQNQ(nu);      
00187     //plots.FillDPIdSigmaQPPassUVVtxCutPlots(nu);//fill these anyway
00188 
00189     //cut on the PID
00190     if (!lib.cuts.IsGoodPID(nu)) {
00191       //plots.FillDPIdSigmaQPFailDpIDCutPlots(nu);
00192       continue;
00193     }
00194     lib.cnt.goodPIDCounter++;   
00195     //plots.FillDPIdSigmaQPPassDpIDCutPlots(nu);
00196 
00197     //cut on the fractional track momentum and sign error      
00198     if (!lib.cuts.IsGoodSigmaQP_QP(nu)) {
00199       //plots.FillDPIdSigmaQPFailSigQPCutPlots(nu);
00200       continue;
00201     }
00202     lib.cnt.goodFitSigQPCounter++;
00203     //plots.FillDPIdSigmaQPPassSigQPCutPlots(nu);
00204     
00205     //plots.FillStpfitqpPlots(ntp,trk,nu);
00206     //plots.FillLowProbPlots(trk,nu);
00207     
00208     //cut on the track fit probability      
00209     if (!lib.cuts.IsGoodFitProb(nu)) {
00210       //plots.FillDPIdSigmaQPFailProbCutPlots(nu);
00211       continue;
00212     }
00213     lib.cnt.goodFitProbCounter++;
00214     
00215     //plots.FillStpfitqpPlots(ntp,trk,nu);
00216     
00217     if (nu.charge==-1) {
00218       //plots.PrintEventInfo(nmTxt,nu);
00219       //time.vTimeNuMuEvt.push_back(evTime);
00220       lib.cnt.nuNQCounter++;
00221     }
00222     else if (nu.charge==+1) {
00223       //plots.PrintEventInfo(nmbTxt,nu);
00224       lib.cnt.CountTrkStdhepId(ntp,evt,trk,nu);
00225       //time.vTimeNuMuBarEvt.push_back(evTime);
00226       lib.cnt.nuPQCounter++;
00227     }
00228     else cout<<"ahhh, bad charge"<<endl;
00229     
00230     //fill histograms etc here...
00231     
00232     //plots.FillRecoEnYHistosN(nu);
00233     //plots.FillDPIdSigmaQPPlotsN(nu);
00234     //plots.FillEnergyBinHistos(nu);
00235     
00236     //plots.FillShwHistos(ntp,evt,nu);
00237     //plots.FillXYZHistos(nu);
00238     //plots.FillContainmentHistos(nu);
00239     //plots.FillTruePIDHistos(nu);
00240     //plots.FillTrackResponseHistos(ntp,trk,nu);
00241     //plots.FillRangeCurvCompHistos(nu);
00242     //plots.FillKinematicsHistos(nu);
00243     
00244     
00245     if (nu.charge<0) {
00246       //it's a NuMu
00247     }
00248     else if (nu.charge>0) {
00249       //it's a NuMuBar
00250     }
00251     else cout<<"ahhh, zero qp"<<endl;
00252     
00253   }//end of for                                       
00254   
00258 
00259 
00260   return JobCResult::kPassed;
00261 }
00262 
00263 //____________________________________________________________________72
00264 
00265 const Registry& NuDemoModule::DefaultConfig() const
00266 {
00268 
00269   MSG("NuDemoModule", Msg::kDebug)
00270     << "Running NuDemoModule::DefaultConfig..." << endl;
00271     
00272   static Registry r;
00273 
00274   string name=this->JobCModule::GetName();
00275   name+=".config.default";
00276   r.SetName(name.c_str());
00277   r.UnLockValues();
00278 
00279   r.LockValues();
00280 
00281   return r;
00282 }
00283 
00284 //____________________________________________________________________72
00285 
00286 void NuDemoModule::Config(const Registry& /*r*/)
00287 {
00289 
00290   MSG("NuDemoModule", Msg::kDebug)
00291     << "Running NuDemoModule::Config..." << endl;
00292 }