#include <NuReco.h>
Public Member Functions | |
| NuReco () | |
| ~NuReco () | |
| void | ApplyReweights (NuEvent &nu) const |
| void | CalcEvtVariables (NuEvent &nu) const |
| void | CalcExtraTruthVariables (NuEvent &nu) const |
| void | CalcKinematicVariables (NuEvent &nu) const |
| void | CalcShwVariables (NuEvent &nu, bool includekNN) const |
| void | CalcTrkReclamation (NuEvent &nu) const |
| void | CalcTrkVariables (NuEvent &nu) const |
| void | CalcResolution (NuEvent &nu) const |
| void | CapTrackMomentum (Double_t &trkMom) const |
| void | ChooseTruthIndexToUse (NuEvent &nu) const |
| TVector3 | xyz2uvz (Float_t x, Float_t y, Float_t z, VldContext vc) const |
| TVector3 | uvz2xyz (Float_t u, Float_t v, Float_t z, VldContext vc) const |
| Int_t | FDRCBoundary (const NuEvent &nu) const |
| Float_t | GetSmallestDeepDistToEdge (const NuEvent &nu) const |
| void | TestGetSmallestDeepDistToEdge () const |
| Int_t | GetBestShower (const NuEvent &nu) const |
| Int_t | GetBestTrack (const NuEvent &nu) const |
| void | GetBestTruthIndex (const NtpStRecord &ntp, const NtpSREvent &evt, NuEvent &nu) const |
| void | GetPrimaryTruthIndex (const NtpStRecord &ntp, const NtpSREvent &evt, NuEvent &nu) const |
| void | GetSecondaryTruthIndex (const NtpStRecord &ntp, const NtpSREvent &evt, NuEvent &nu) const |
| Int_t | GetContainmentFlag (const NuEvent &nu) const |
| Int_t | GetContainmentFlagCC0093Std (const NuEvent &nu) const |
| Int_t | GetContainmentFlagCC0250Std (const NuEvent &nu) const |
| Int_t | GetContainmentFlagPitt (const NuEvent &nu) const |
| Int_t | GetContainmentFlagStdReco (const NuEvent &nu) const |
| void | GetDirCosNu (const NtpSRTrack &trk, NuEvent &nu) const |
| Float_t | GetDirCosNu (const NuEvent &nu) const |
| Double_t | GetCosBetweenPr_Theta (const NtpSRTrack &trk) const |
| Retrieves the cos of the angle between the radial momentum and track vertex. | |
| Double_t | GetEvtMedianTime (const NtpStRecord &ntp, const NtpSREvent &evt) const |
| Double_t | GetShwMedianTime (const NtpStRecord &ntp, const NtpSRShower &shw) const |
| Int_t | GetShwMaxPlane (const NtpStRecord &ntp, const NtpSRShower &shw) const |
| Which plane does shw have the most ph deposited in? | |
| void | GetEvtDeltaTs (const NtpStRecord &ntp, std::map< Int_t, Double_t > &deltaTs) const |
| Double_t | GetEvtEnergy (NuEvent &nu, bool includekNN) const |
| void | GetEvtsPerSlc (const NtpStRecord &ntp, std::map< Int_t, Int_t > &evtsPerSlc) const |
| Int_t | GetHadronicFinalState (const NtpStRecord &ntp, const NuEvent &nu) const |
| void | GetPreInukeFinalStateVars (const NtpStRecord &ntp, NuEvent &nu) const |
| Int_t | GetInitialState (const NtpStRecord &ntp, const NuEvent &nu) const |
| const NtpSRTrack * | GetLongestTrack (const NtpStRecord &ntp, const NtpSREvent &evt) const |
| Int_t | GetLongestTrack (const NuEvent &nu) const |
| const NtpSRTrack * | GetLongestTrackTLikePlanes (const NtpStRecord &ntp, const NtpSREvent &evt) const |
| Int_t | GetNucleus (const NtpStRecord &ntp, const NuEvent &nu) const |
| Int_t | GetPrimaryShowerCCStd (const NuEvent &nu) const |
| Int_t | GetPrimaryShowerStdReco (const NuEvent &nu) const |
| const NtpSRTrack * | GetTrackWithIndexX (const NtpStRecord &ntp, const NtpSREvent &evt, Int_t index) const |
| Int_t | GetPrimaryTrack (const NuEvent &nu) const |
| Float_t | GetRadius (Float_t x, Float_t y, const NuEvent &nu) const |
| Double_t | GetShowerEnergyCC (const NuEvent &nu) const |
| Double_t | GetShowerEnergyNC (const NuEvent &nu) const |
| Double_t | GetShowerEnergyCor (Double_t shwEn, const CandShowerHandle::ShowerType_t &st, const NuEvent &nu) const |
| bool | InitializeShowerEnergykNN (NuKDTree< double, 3 > *kdTree, Detector::Detector_t det, std::vector< double > &C, double &cutoff_lo, double &cutoff_hi, int &neighbours) const |
| Helper function for GetShowerEnergykNN. | |
| Double_t | GetShowerEnergykNN (const NuEvent &num, Float_t *res=0) const |
| res will be filled with an energy resolution estimate | |
| Double_t | GetShowerEnergyNearTrack (const NtpStRecord &ntp, const NtpSRTrack &trk, Float_t *nearShowerEnergyDW=0) const |
| Calculates the shower energy 'near' the track vertex. | |
| const NtpSRShower * | GetShowerWithIndexX (const NtpStRecord &ntp, const NtpSREvent &evt, Int_t index) const |
| Int_t | GetTrackCharge (const NtpSRTrack &trk, bool isReverse) const |
| Int_t | GetTrackCharge (Double_t qp, bool isReverse) const |
| Double_t | GetTrackEnergy (const NuEvent &nu) const |
| Double_t | GetTrackEnergyFromRange (const NuEvent &nu) const |
| Double_t | GetTrackEnergyFromRangeCor (Double_t trkMomentumRange, const NuEvent &nu) const |
| Double_t | GetTrackEnergyFromCurv (const NuEvent &nu) const |
| Double_t | GetTrackEnergyFromCurvCor (Double_t qp, const NuEvent &nu) const |
| Double_t | GetTrackTransverseMomentum (const NuEvent &nu) const |
| Double_t | GetTrkQPFraction (const NtpSRTrack &trk) const |
| void | GetTruthInfo (const NtpStRecord &ntp, const NtpSREvent &evt, NuEvent &nu) const |
| void | GetTruthInfo (const NtpStRecord &ntp, const NtpMCTruth &mc, NuEvent &nu) const |
| VldContext | GetVldContext (const NuEvent &nu) const |
| Bool_t | IsTrkWithDimuon (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | IsTrkWithMuonFromKaonDecay (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | IsTrkWithMuonFromPionDecay (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | IsTrkWithMuonFromNotNuNotPiNotKaonNotCharm (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | PrintCharm (const NtpStRecord &ntp, const NtpSREvent &evt) const |
| void | PrintStdhep (const NtpStRecord &ntp, const NtpMCTruth &mc) const |
| void | PrintTrueEnergy (const NtpStRecord &ntp, const NtpSREvent &evt, Double_t recoEn=0) const |
| void | SetBestShw (NuEvent &nu) const |
| void | SetBestTrk (NuEvent &nu) const |
| void | SetBestTrkMajorityCurvature (NuEvent &nu) const |
| void | SetBestTrkSAFit (NuEvent &nu) const |
| Bool_t | UseCurvature (const NuEvent &nu) const |
| Bool_t | UseRange (const NuEvent &nu) const |
| Int_t | GetPrimaryMuonStdHepIndex (const NtpStRecord &ntp, const NtpMCTruth &mc) const |
| Find the stdhep array entry that corresponds to the muon. | |
| Double_t | GetMuonFirstHitEnergy (const NtpStRecord &ntp, const NtpMCTruth &mc) const |
| Retrieves the energy of muon on the first detector hit. | |
| Double_t | GetMuonLastHitEnergy (const NtpStRecord &ntp, const NtpMCTruth &mc) const |
| Retrieves the energy of muon on the last detector hit. | |
| Bool_t | GetMuonContainmentMC (const NtpStRecord &ntp, const NtpMCTruth &mc) const |
| Experimental: Attempts to calculate the 'containment truth'. | |
| Int_t | GetTrackFirstStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Int_t | GetTrackFirstUStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Int_t | GetTrackFirstVStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | GetTrackFirstStripIsU (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Int_t | GetTrackLastStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Bool_t | GetTrackLastStripIsU (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Int_t | GetTrackLastUStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Int_t | GetTrackLastVStrip (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Float_t | GetPhi (const NtpSRVertex &vtx) const |
| Int_t | GetRegion (const NtpStRecord &rec, const NtpSRVertex &vtx) const |
| Retrieves the (rock muon) detector region from SNTP. | |
| Int_t | GetRegion (Detector::Detector_t det, SimFlag::SimFlag_t simFlag, Float_t vtxX, Float_t vtxY, Float_t vtxZ, Int_t vtxPlane) const |
| Int_t | GetEdgeRegion (const NtpStRecord &rec, const NtpSRVertex &vertex) const |
| Retrieves the FD edge region from the SNTP. | |
| Int_t | GetEdgeRegion (Detector::Detector_t det, const Float_t phi) const |
| Calculates the edge region based on direct variables. | |
| void | CalculateEdgeRegion (NuEvent &nu, Int_t track) const |
| Short_t | GetParallelStripInset (Detector::Detector_t det, Short_t edgeRegion, Short_t stripBegU, Short_t stripBegV) const |
| How many strips in, on a diagonal edge, is a hit? | |
| void | CalculateParallelStripInset (NuEvent &nu, Int_t track=0) const |
| Given a NuEvent, calculates from existing variables the parallel strip inset. | |
| Char_t | GetPerpendicularStripFlag (Detector::Detector_t det, Short_t edgeRegion, Bool_t IsHitu) const |
| Is the hit both on a diagonal edge and on the overhanging (perpendicular) strip? | |
| void | CalculatePerpendicularStripFlag (NuEvent &nu, Int_t track) const |
| Given a NuEvent, calculates the Perpendicular strip flags. | |
| Short_t | GetHorizontalVerticalStripNumber (Detector::Detector_t det, Short_t edgeRegion, Short_t stripTrku, Short_t stripTrkv) const |
| void | CalculateHorizontalVerticalStripNumber (NuEvent &nu, Int_t track) const |
| Given a NuEvent, calculates the HOVE strip numbers. | |
| void | RecalculateDerivativeTrackGeometryVariables (NuEvent &nu) const |
| Recalculates all of the base geometry variables it can, given only a NuEvent. | |
| Int_t | GetTrueTrackId (const NtpStRecord &ntp, const NtpSRTrack &trk) const |
| Get the true track identity. | |
|
|
Definition at line 57 of file NuReco.cxx. References MSG. 00058 {
00059 MSG("NuReco",Msg::kDebug)
00060 <<"Running NuReco Constructor..."<<endl;
00061
00062 MSG("NuReco",Msg::kDebug)
00063 <<"Finished NuReco Constructor"<<endl;
00064 }
|
|
|
Definition at line 67 of file NuReco.cxx. References MSG. 00068 {
00069 MSG("NuReco",Msg::kDebug)
00070 <<"Running NuReco Destructor..."<<endl;
00071
00072
00073 MSG("NuReco",Msg::kDebug)
00074 <<"Finished NuReco Destructor"<<endl;
00075 }
|
|
|
Definition at line 389 of file NuReco.cxx. References NuEvent::anaVersion, NuEvent::apply1SigmaWeight, NuEvent::applyBeamWeight, NuEvent::applyDetectorWeight, NuEvent::applyEnergyShifts, NuEvent::applyGeneratorWeight, NuEvent::beamWeight, CalcResolution(), NuEvent::detectorWeight, NuEvent::detectorWeightNM, NuEvent::detectorWeightNMB, NuEvent::energy, NuEvent::energyRw, NuEvent::fluxErr, NuEvent::fluxErrTotalErrorAfterTune, NuEvent::generatorWeight, NuCuts::IsNMBPQ(), MAXMSG, NuEvent::rw, NuEvent::rwActual, NuEvent::shwEn, NuEvent::shwEnNoRw, NuEvent::shwEnRw, NuEvent::shwEnWeight, NuEvent::simFlag, NuEvent::trkEn, NuEvent::trkEnNoRw, NuEvent::trkEnRw, and NuEvent::trkEnWeight. Referenced by NuDSTAna::BRevAna(), NuDSTAna::Contamination(), NuDSTAna::DPSystematic(), NuAnalysis::ExtractPIDsAndWeights(), NuDSTAna::FDTestAna(), NuDSTAna::FluxComponents(), NuDSTAna::JeffsTestAna(), NuAnalysis::LIRejectionTest(), NuAnalysis::LoopOverTruthInfo(), NuDSTAna::MakeFCTree(), NuDSTAna::MakeMicroDST(), NuDSTAna::MakeMicroDst2010(), NuDSTAna::MakeMicroDstFakeData(), NuDSTAna::MakeMicroDstForCSSSystematics(), NuDSTAna::MakeMicroDstHe(), NuDSTAna::MakeMicroDstJJEPresel(), NuDSTAna::MakeMicroDstWithStdCCRecoAndCuts(), NuDSTAna::MakeSelMicroDST(), NuDSTAna::MMRereco(), NuDSTAna::MMTransition(), NuDSTAna::NDOsc(), NuDSTAna::NDQPRB(), NuDSTAna::NDTestAna(), NuDSTAna::NewFieldAna(), NuDSTAna::QPStudy(), NuDSTAna::RHCTest(), NuDSTAna::StdCCAna(), and NuDSTAna::StdNMBAna(). 00390 {
00391 if (nu.simFlag==SimFlag::kData) {
00392 MAXMSG("NuReco",Msg::kInfo,1)
00393 <<"Not doing ApplyReweights for simFlag="<<nu.simFlag<<endl;
00394
00395 nu.trkEnRw=nu.trkEnNoRw;
00396 nu.shwEnRw=nu.shwEnNoRw;
00397 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00398 return;
00399 }
00400
00401 //shift the energies
00402 nu.trkEnRw=nu.trkEnNoRw*nu.trkEnWeight;
00403 nu.shwEnRw=nu.shwEnNoRw*nu.shwEnWeight;
00404 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00405
00406 //check whether to apply energy shifts
00407 if (nu.applyEnergyShifts) {
00408 MAXMSG("NuReco",Msg::kInfo,1)
00409 <<"Applying SKZP energy shifts"<<endl;
00410 nu.trkEn=nu.trkEnRw;
00411 nu.shwEn=nu.shwEnRw;
00412 nu.energy=nu.energyRw;
00413 }
00414 else {
00415 MAXMSG("NuReco",Msg::kInfo,1)
00416 <<"NOT applying SKZP energy shifts"<<endl;
00417 }
00418
00419 //set the detector weight to be either the NuMu or NuMuBar weight
00420 //because of wc.SetSampleSelection("NuMuBar");
00421 if (nu.anaVersion==NuCuts::kJJH1 ||
00422 nu.anaVersion==NuCuts::kJJE1 ||
00423 nu.anaVersion==NuCuts::kFullDST || NuCuts::IsNMBPQ(nu) ) {
00424 MAXMSG("NuReco",Msg::kInfo,1)
00425 <<"Setting the detector weight for NuMuBar selection"<<endl;
00426 nu.detectorWeight=nu.detectorWeightNMB;
00427 }
00428 else {
00429 MAXMSG("NuReco",Msg::kInfo,1)
00430 <<"Setting the detector weight for NuMu selection"<<endl;
00431 nu.detectorWeight=nu.detectorWeightNM;
00432 }
00433
00434 //reset to 1 at start
00435 //shouldn't have to do this except when re-weighting summary trees
00436 nu.rw=1;
00437
00438 //check whether to apply beam weight
00439 if (nu.applyBeamWeight) {
00440 MAXMSG("NuReco",Msg::kInfo,1)
00441 <<"Applying beam weight"<<endl;
00442 nu.rw*=nu.beamWeight;
00443 MAXMSG("NuReco",Msg::kInfo,1)
00444 <<"Beam weight is " << nu.rw <<endl;
00445 }
00446 else {
00447 MAXMSG("NuReco",Msg::kInfo,1)
00448 <<"NOT applying beam weight"<<endl;
00449 }
00450
00451 //check whether to apply detector weight
00452 if (nu.applyDetectorWeight) {
00453 MAXMSG("NuReco",Msg::kInfo,1)
00454 <<"Applying detector weight"<<endl;
00455 nu.rw*=nu.detectorWeight;
00456 }
00457 else {
00458 MAXMSG("NuReco",Msg::kInfo,1)
00459 <<"NOT applying detector weight"<<endl;
00460 }
00461
00462 //check whether to apply generator weight
00463 if (nu.applyGeneratorWeight) {
00464 MAXMSG("NuReco",Msg::kInfo,1)
00465 <<"Applying generator weight"<<endl;
00466 nu.rw*=nu.generatorWeight;
00467 }
00468 else {
00469 MAXMSG("NuReco",Msg::kInfo,1)
00470 <<"NOT applying generator weight"<<endl;
00471 }
00472
00473 //check whether to apply 1-sigma weight
00474 if (nu.apply1SigmaWeight) {
00475 MAXMSG("NuReco",Msg::kInfo,1)
00476 <<"Applying 1-sigma weight"<<endl;
00477 nu.rw*=(nu.fluxErrTotalErrorAfterTune-1.0)+1.0;
00478 }
00479 else {
00480 MAXMSG("NuReco",Msg::kInfo,1)
00481 <<"NOT applying 1-sigma weight"<<endl;
00482 }
00483 //this stores what was actually used
00484 nu.rwActual=nu.rw;
00485
00487 //now choose the flux error to use
00488 nu.fluxErr=nu.fluxErrTotalErrorAfterTune;
00489
00490 //Recalculate resolution with reweighted/shifted variables
00491 this->CalcResolution(nu);
00492
00493 }
|
|
|
Definition at line 202 of file NuReco.cxx. References NuEvent::distToEdgeEvtVtx, NuEvent::evtVtxUVDiffPl, GetRadius(), GetSmallestDeepDistToEdge(), NuEvent::planeEvtBegu, NuEvent::planeEvtBegv, NuEvent::rEvtEnd, NuEvent::rEvtVtx, NuEvent::xEvtEnd, NuEvent::xEvtVtx, NuEvent::yEvtEnd, and NuEvent::yEvtVtx. Referenced by GetEvtEnergy(). 00203 {
00204 //calculate vtx/end positions
00205 nu.rEvtVtx=this->GetRadius(nu.xEvtVtx,nu.yEvtVtx,nu);
00206 nu.rEvtEnd=this->GetRadius(nu.xEvtEnd,nu.yEvtEnd,nu);
00207 nu.evtVtxUVDiffPl=nu.planeEvtBegu-nu.planeEvtBegv;
00208 nu.distToEdgeEvtVtx=this->GetSmallestDeepDistToEdge(nu);
00209 }
|
|
|
Definition at line 366 of file NuReco.cxx. References PlexPlaneId::GetPlane(), UgliGeomHandle::GetPlaneIdFromZ(), GetRadius(), GetVldContext(), NuLibrary::Instance(), NuEvent::planeTrkVtxMC, NuLibrary::reco, NuEvent::rTrkVtxMC, NuEvent::vtxuMC, NuEvent::vtxvMC, NuEvent::vtxxMC, NuEvent::vtxyMC, NuEvent::vtxzMC, and xyz2uvz(). Referenced by GetTruthInfo(). 00367 {
00368 //get ugh
00369 const VldContext& vc=this->GetVldContext(nu);
00370 UgliGeomHandle ugh(vc);
00371
00372 //get an instance of the code library
00373 const NuLibrary& lib=NuLibrary::Instance();
00374
00375 //calculate the positions in UVZ space
00376 TVector3 uvz=lib.reco.xyz2uvz(nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,vc);
00377
00378 //write out the variables
00379 nu.vtxuMC=uvz.X();
00380 nu.vtxvMC=uvz.Y();
00381
00382 //calc the other variables
00383 nu.planeTrkVtxMC=ugh.GetPlaneIdFromZ(nu.vtxzMC).GetPlane();
00384 nu.rTrkVtxMC=lib.reco.GetRadius(nu.vtxxMC,nu.vtxyMC,nu);
00385 }
|
|
|
Definition at line 178 of file NuReco.cxx. References NuEvent::dirCosNu, NuEvent::energy, GetDirCosNu(), NuEvent::q2, NuEvent::shwEn, NuEvent::trkEn, NuEvent::w2, NuEvent::x, and NuEvent::y. Referenced by GetEvtEnergy(). 00179 {
00180 nu.dirCosNu=this->GetDirCosNu(nu);
00181
00183 //if(reco_emu>0) reco_Q2 =
00184 // 2*reco_enu*reco_emu*(1.0 - reco_dircosneu);
00185 //reco_W2 = M*M - reco_Q2 + 2*M*reco_eshw;
00186 //if(reco_enu>0) reco_y = reco_eshw/reco_enu;
00187 //if(reco_eshw>0 && reco_Q2>0) reco_x = reco_Q2/(2*M*reco_eshw);
00188
00189 //calc kinematics
00190 const Double_t M=(0.93827 + 0.93957)/2.0; // <nucleon mass>
00191
00192 //calc y, q2, w2 and x
00193 if (nu.trkEn+nu.shwEn) nu.y=nu.shwEn/(nu.trkEn+nu.shwEn);
00194 nu.q2=2*nu.energy*nu.trkEn*(1.0-nu.dirCosNu);
00195 nu.w2=M*M-nu.q2+2*M*nu.shwEn;
00196 if (nu.shwEn) nu.x=nu.q2/(2*M*nu.shwEn);
00197 //else nu.x=0;
00198 }
|
|
|
||||||||||||
|
|
Definition at line 271 of file NuReco.cxx. References NuEvent::detector, MAXMSG, NuEvent::trkEn, NuEvent::trkEnCurv, NuEvent::trkEnNoRw, NuEvent::trkEnRange, NuEvent::trkExists, NuEvent::trkfitpass, and NuEvent::usedCurv. Referenced by GetEvtEnergy(). 00272 {
00273 if (nu.trkExists && nu.trkEnCurv==0) {
00274 //sanity check for the FD
00275 if (nu.detector==Detector::kFar) {
00276 MAXMSG("NuReco",Msg::kDebug,100)
00277 <<"Looks like track reclamation for the FD..."
00278 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00279 <<", trkEnRange="<<nu.trkEnRange
00280 <<", trkfitpass="<<nu.trkfitpass
00281 <<endl;
00282 }
00283
00284 MAXMSG("NuReco",Msg::kDebug,10)
00285 <<endl<<"Setting curv=range"
00286 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00287 <<", trkEnRange="<<nu.trkEnRange<<", trkfitpass="<<nu.trkfitpass
00288 <<endl;
00289 nu.trkEnCurv=nu.trkEnRange;
00290 nu.trkEn=nu.trkEnRange;
00291 nu.trkEnNoRw=nu.trkEn;
00292 }
00293 else if (!nu.trkExists && nu.trkEnCurv==0) {
00294 MAXMSG("NuReco",Msg::kDebug,10)
00295 <<"nu.trkExists="<<nu.trkExists<<", nu.trkEnCurv="<<nu.trkEnCurv
00296 <<endl;
00297 }
00298 }
|
|
|
||||||||||||
|
||||||||||||
|
||||||||||||
|
||||||||||||
|
|
Definition at line 573 of file NuReco.cxx. Referenced by GetTrackEnergyFromCurvCor(), and GetTrackEnergyFromRangeCor(). 00574 {
00575 //cap the trkMom to avoid crazy values e.g. 1e19 causing fpes)
00576 if (trkMom>(1000*Munits::GeV)) trkMom=1000*Munits::GeV;
00577 else if (trkMom<(-1000*Munits::GeV)) trkMom=-1000*Munits::GeV;
00578 }
|
|
|
Definition at line 3127 of file NuReco.cxx. References NuEvent::mc, NuEvent::mcShw, and NuEvent::mcTrk. Referenced by GetBestTruthIndex(), and GetPrimaryTruthIndex(). 03128 {
03129 //use the track as the primary determination
03130 //else use the event
03131 //this can be changed for particular anaVersion as required
03132 //if (mcTrk>=0) nu.mc=nu.mcTrk;
03133 //else nu.mc=nu.mcEvt;
03134
03135 //UPDATE 2007/12/10: use the mcEvt as the best guess for the mc index
03136 //this is what the generator reweighting example uses
03137 //nu.mc=nu.mcEvt;
03138
03139 //UPDATE 2008/02/20: use trk or shw like the PANs do
03140 if (nu.mcTrk>=0) nu.mc=nu.mcTrk;
03141 else nu.mc=nu.mcShw;
03142 }
|
|
|
litag takes the value of 0,1,10,11, usually cut all >0 Definition at line 497 of file NuReco.cxx. References NuEvent::detector, NuEvent::nshw, NuEvent::planeEvtHdrBeg, NuEvent::planeEvtHdrEnd, and NuEvent::rawPhEvt. Referenced by NuExtraction::ExtractLITags(), and NuAnalysis::LIRejectionTest(). 00498 {
00500 Int_t litag=0;
00501
00502 if(nu.detector!=Detector::kFar) return litag;
00503
00504 //Int_t numshower=eventSummary->nshower;
00505 //Float_t allph=eventSummary->ph.raw;
00506 //Int_t plbeg=eventSummary->plane.beg;
00507 //Int_t plend=eventSummary->plane.end;
00508 Int_t numshower=nu.nshw;
00509 Float_t allph=nu.rawPhEvt;
00510 //Int_t plbeg=nu.planeEvtBeg;
00511 //Int_t plend=nu.planeEvtEnd;
00512 Int_t plbeg=nu.planeEvtHdrBeg;//evthdr.plane.beg (diff. to above)
00513 Int_t plend=nu.planeEvtHdrEnd;//evthdr.plane.end
00514
00515 if (numshower) {
00516 if (allph>1e6) litag+=10;
00517
00518 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
00519 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
00520 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
00521 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
00522 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
00523 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
00524 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
00525 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
00526 }
00527 return litag;
00528 }
|
|
|
work out the best shower to use in the case that there is more than one Definition at line 1900 of file NuReco.cxx. References NuEvent::anaVersion, GetPrimaryShowerCCStd(), GetPrimaryShowerStdReco(), NuCuts::IsNMBNQ(), NuCuts::IsNMBPQ(), and MAXMSG. Referenced by NuCounter::CountTrkStdhepId(), NuAnalysis::DemoInfidSRInterface(), NuPlots::FillShwHistos(), GetBestTruthIndex(), and SetBestShw(). 01901 {
01904
01905 if (nu.anaVersion==NuCuts::kJJH1 ||
01906 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMBPQ(nu) ) {
01907 MAXMSG("NuReco",Msg::kInfo,1)
01908 <<"Setting best shower as primary shower (from reconstruction)"
01909 <<endl;
01910 return this->GetPrimaryShowerStdReco(nu);
01911 }
01912 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01913 MAXMSG("NuReco",Msg::kInfo,1)
01914 <<"Setting best shower according to CC std"<<endl;
01915 return this->GetPrimaryShowerCCStd(nu);
01916 }
01917 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01918 nu.anaVersion==NuCuts::kCC0325Std ||
01919 nu.anaVersion==NuCuts::kCC0720Std ||
01920 nu.anaVersion==NuCuts::kCC0720Test ||
01921 NuCuts::IsNMBNQ(nu) ) {
01922 MAXMSG("NuReco",Msg::kInfo,1)
01923 <<"Setting best shower according to CC std"<<endl;
01924 return this->GetPrimaryShowerCCStd(nu);
01925 }
01926 else {
01927 MAXMSG("NuReco",Msg::kInfo,1)
01928 <<"Setting best shower as primary shower (from reconstruction)"
01929 <<endl;
01930 return this->GetPrimaryShowerStdReco(nu);
01931 }
01932 }
|
|
|
can be either the longest track or the one defined as the primary track by the reconstruction (evt.trk[0]) Definition at line 1828 of file NuReco.cxx. References NuEvent::anaVersion, GetLongestTrack(), GetPrimaryTrack(), NuCuts::IsNMB(), and MAXMSG. Referenced by NuDemoModule::Ana(), NuAnalysis::ChargeSeparationOneSnarl(), NuAnalysis::ChargeSignCut(), NuAnalysis::DemoInfidSRInterface(), NuAnalysis::Efficiencies(), NuAnalysis::EnergySpect(), NuPlots::FillShwHistos(), GetBestTruthIndex(), NuAnalysis::N_1(), SetBestTrk(), SetBestTrkMajorityCurvature(), and SetBestTrkSAFit(). 01829 {
01830 //work out the best track to use in the case that there is
01831 //more than one
01834
01835 if (nu.anaVersion==NuCuts::kJJH1 ||
01836 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMB(nu) ) {
01837 MAXMSG("NuReco",Msg::kInfo,1)
01838 <<"Setting best track as longest track"<<endl;
01839 return this->GetLongestTrack(nu);
01840 }
01841 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01842 MAXMSG("NuReco",Msg::kInfo,1)
01843 <<"Setting best track as longest track"<<endl;
01844 return this->GetLongestTrack(nu);
01845 }
01846 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01847 nu.anaVersion==NuCuts::kCC0325Std ||
01848 nu.anaVersion==NuCuts::kCC0720Std ||
01849 nu.anaVersion==NuCuts::kCC0720Test ) {
01850 MAXMSG("NuReco",Msg::kInfo,1)
01851 <<"Setting best track as longest track"<<endl;
01852 return this->GetLongestTrack(nu);
01853 }
01854 else {
01855 MAXMSG("NuReco",Msg::kInfo,1)
01856 <<"Setting best track as primary track (from reconstruction)"
01857 <<endl;
01858 return this->GetPrimaryTrack(nu);
01859 }
01860 }
|
|
||||||||||||||||
|
Definition at line 2887 of file NuReco.cxx. References ChooseTruthIndexToUse(), GetBestShower(), GetBestTrack(), GetShowerWithIndexX(), GetTrackWithIndexX(), NtpSRTrack::index, NtpSRShower::index, NtpSREvent::index, MAXMSG, NuEvent::mcEvt, NuEvent::mcShw, NuEvent::mcTrk, NtpTHTrack::neumc, NtpTHShower::neumc, NtpTHEvent::neumc, NtpStRecord::thevt, NtpStRecord::thshw, and NtpStRecord::thtrk. 02889 {
02890 TClonesArray& thevtTca=(*ntp.thevt);//TruthHelper Event TCA
02891 const Int_t numthevts=thevtTca.GetEntriesFast();
02892 if (numthevts<=0) {
02893 MAXMSG("NuReco",Msg::kDebug,1)
02894 <<"No THEvents, so can't GetBestTruthIndex..."<<endl;
02895 return;
02896 }
02897 MAXMSG("NuReco",Msg::kDebug,1)
02898 <<"Found THEvent, GetBestTruthIndex..."<<endl;
02899
02900 Int_t mcTrk=-1;//track mc index
02901 Int_t mcShw=-1;//shower mc index
02902 Int_t mcEvt=-1;//event mc index
02903
02904 const NtpTHEvent& thevt=
02905 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02906 mcEvt=thevt.neumc;
02907
02908 //get an instance of the code library
02909 //const NuLibrary& lib=NuLibrary::Instance();
02910
02911 //get the shower mc
02912 TClonesArray& thshwTca=(*ntp.thshw);//TruthHelper Shower TCA
02913 const Int_t numthshws=thshwTca.GetEntriesFast();
02914 if (numthshws>0) {
02915 Int_t bestShower=this->GetBestShower(nu);
02916 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02917 bestShower-1);
02918 if (pshw) {
02919 const NtpSRShower& shw=*pshw;
02920 const NtpTHShower& thshw=
02921 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02922 mcShw=thshw.neumc;
02923 }
02924 }
02925
02926 //get the track mc
02927 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Track TCA
02928 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02929 if (numthtrks>0) {
02930 Int_t bestTrack=this->GetBestTrack(nu);
02931 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02932 bestTrack-1);
02933 if (ptrk) {
02934 const NtpSRTrack& trk=*ptrk;
02935 const NtpTHTrack& thtrk=
02936 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02937 mcTrk=thtrk.neumc;
02938 }
02939 }
02940
02941 //set the variables in the object
02942 nu.mcTrk=mcTrk;
02943 nu.mcShw=mcShw;
02944 nu.mcEvt=mcEvt;
02945 this->ChooseTruthIndexToUse(nu);
02946 }
|
|
|
Definition at line 3444 of file NuReco.cxx. References NuEvent::anaVersion, NuEvent::containmentFlagCC0093Std, NuEvent::containmentFlagCC0250Std, NuEvent::containmentFlagPitt, NuEvent::detector, GetContainmentFlagStdReco(), ReleaseType::IsBirch(), NuCuts::IsNMBNQ(), MAXMSG, MSG, NuEvent::recoVersion, and NuEvent::xTrkEnd. Referenced by CalcTrkVariables(). 03445 {
03446 Int_t flag=-1;
03447
03448 //now set the containment flag
03449 if (nu.anaVersion==NuCuts::kCC0093Std) {
03450 MAXMSG("NuReco",Msg::kInfo,1)
03451 <<"Using CC0093Std containment criteria"<<endl;
03452 flag=nu.containmentFlagCC0093Std;
03453 }
03454 else if (nu.anaVersion==NuCuts::kCC0250Std) {
03455 MAXMSG("NuReco",Msg::kInfo,1)
03456 <<"Using CC0250Std containment criteria"<<endl;
03457 flag=nu.containmentFlagCC0250Std;
03458 }
03459 else if (nu.anaVersion==NuCuts::kCC0325Std ||
03460 nu.anaVersion==NuCuts::kCC0720Std ||
03461 nu.anaVersion==NuCuts::kCC0720Test||
03462 NuCuts::IsNMBNQ(nu)) {
03463 MAXMSG("NuReco",Msg::kInfo,1)
03464 <<"Using CC0325Std hybrid containment criteria"<<endl;
03465 flag=nu.containmentFlagCC0250Std;
03466 //in ND use the hybrid approach to reduce the bias that occurs for
03467 //events exiting the side of the detector and being reco'd by range
03468 if (nu.detector==Detector::kNear) {
03469 if (nu.xTrkEnd>1.3*Munits::m) {
03470 flag=this->GetContainmentFlagStdReco(nu);
03471 }
03472 }
03473 // MAXMSG("NuReco",Msg::kInfo,1)
03474 // <<"Using SNTP SR Containment Flag"<<endl;
03475 // flag = this->GetContainmentFlagStdReco(nu);
03476 //
03477 // // For events that stop in the coil hole, prefer to use the CC flag
03478 // if (nu.detector == Detector::kFar && nu.rTrkEnd < 0.4) {
03479 // flag=nu.containmentFlagCC0250Std;
03480 // }
03481 }
03482 else {
03483 if (ReleaseType::IsBirch(nu.recoVersion)) {
03484 MAXMSG("NuReco",Msg::kInfo,1)
03485 <<"Using Pitt containment criteria"<<endl;
03486 //note that trk.contained is not available before Cedar
03487 flag=nu.containmentFlagPitt;
03488 }
03489 else {
03490 MAXMSG("NuReco",Msg::kInfo,1)
03491 <<"Using std reco containment criteria"<<endl;
03492 flag=this->GetContainmentFlagStdReco(nu);
03493 //MAXMSG("NuReco",Msg::kInfo,1)
03494 //<<"Using Pitt containment criteria"<<endl;
03495 //flag=nu.containmentFlagPitt;
03496 }
03497 }
03498
03499 //sanity check
03500 if (!(flag==1 || flag==2 || flag==3 || flag==4)) {
03501 MSG("NuReco",Msg::kWarning)<<"flag="<<flag<<endl;
03502 }
03503
03504 //return the containment flag
03505 return flag;
03506 }
|
|
|
Definition at line 3608 of file NuReco.cxx. References NuEvent::detector, NuEvent::drTrkFidall, NuEvent::dzTrkFidall, NuEvent::xTrkEnd, NuEvent::yTrkEnd, and NuEvent::zTrkEnd. Referenced by CalcTrkVariables(). 03609 {
03610 const Float_t x=nu.xTrkEnd;
03611 const Float_t y=nu.yTrkEnd;
03612 const Float_t z=nu.zTrkEnd;
03613 const Float_t r=sqrt(x*x+y*y);
03614
03615 const Float_t xMin=-1.65;
03616 const Float_t xMax=+2.7;
03617 const Float_t yMin=-1.65;
03618 const Float_t yMax=+1.65;
03619 const Float_t zMax=+16;
03620 const Float_t rMin=+0.4;
03621 const Float_t rMinFD=+0.5;
03622
03623 Bool_t contained=false;
03624 Int_t flag=-1;
03625
03626 //taken from the code below
03627 /*
03628 Bool_t MadQuantities::IsFidAll(Int_t itrk){
03629 if (!(ntpTrack->end.z<16 && sqrt(pow(ntpTrack->end.x,2)+
03630 pow(ntpTrack->end.y,2))>0.4 &&
03631 ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 &&
03632 ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
03633
03634 ntpTrack->end.y>(-ntpTrack->end.x)-1.65 &&
03635 ntpTrack->end.y<ntpTrack->end.x+1.65 &&
03636 ntpTrack->end.y<(-ntpTrack->end.x)+3.55 &&
03637 ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
03638
03639
03640 }
03641 else if(ntpTrack->fidall.dr<0.5 ||
03642 ntpTrack->fidall.dz<0.5) return false;
03643 */
03644
03645 if (nu.detector==Detector::kNear) {
03646 //calculate if contained
03647 //don't understand the last 4 lines... ND shape?
03648 contained=z<zMax && r>rMin &&
03649 x>xMin && x<xMax &&
03650 y>yMin && y<yMax &&
03651 y>-x+xMin &&
03652 y<+x-xMin &&
03653 y<-x+3.55 &&
03654 y>+x-3.55;
03655
03656 if (contained) {
03657 //117=~6.95 m
03658 //118=~7.01 m
03659 //119=~7.07 m
03660 //120=~7.13 m (partial)
03661 //121=~7.18 m (full)
03662 //122=~7.24 m
03663 if (z<7) flag=1;//planes 1->117 inclusive
03664 else if (z>=7) flag=3;//planes 118 onwards
03665 }
03666 else {
03667 if (z<7) flag=2;
03668 else if (z>=7) flag=4;
03669 }
03670 }
03671 else if (nu.detector==Detector::kFar) {
03672 contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
03673 if (contained) flag=1;
03674 else flag=2;
03675 }
03676 else cout<<"Ahh, detector="<<nu.detector<<endl;
03677
03678 return flag;
03679 }
|
|
|
Definition at line 3543 of file NuReco.cxx. References NuEvent::detector, NuEvent::planeTrkEnd, NuEvent::xTrkEnd, NuEvent::yTrkEnd, and NuEvent::zTrkEnd. Referenced by CalcTrkVariables(). 03544 {
03545 const Float_t x=nu.xTrkEnd;
03546 const Float_t y=nu.yTrkEnd;
03547 const Float_t z=nu.zTrkEnd;
03548 const Float_t r=sqrt(x*x+y*y);//FD only
03549
03550 //minos-doc-3156 has these values:
03551 //ND cuts:
03552 const Float_t xMin=-1.65*(Munits::m);
03553 const Float_t xMax=+2.7*(Munits::m);
03554 const Float_t yMin=-1.65*(Munits::m);
03555 const Float_t yMax=+1.65*(Munits::m);
03556 //const Float_t uMin=-1.65*(Munits::m);
03557 const Float_t uMax=+3.55*(Munits::m);
03558 //const Float_t vMin=-3.55*(Munits::m);
03559 //const Float_t vMax=+1.65*(Munits::m);
03560 const Float_t zMax=+15*(Munits::m);
03561
03562 //FD cuts
03563 const Float_t planeMin=4;
03564 const Float_t planeMax=475;
03565 const Float_t rMax=TMath::Sqrt(14)*(Munits::m);
03566
03567 Bool_t contained=false;
03568 Int_t flag=-1;
03569
03570 if (nu.detector==Detector::kNear) {
03571 contained=z<zMax &&
03572 x>xMin && x<xMax &&
03573 y>yMin && y<yMax &&
03574 y>-x+xMin &&
03575 y<+x-xMin &&
03576 y<-x+uMax &&
03577 y>+x-uMax;
03578
03579 if (contained) {
03580 //117=~6.95 m
03581 //118=~7.01 m
03582 //119=~7.07 m
03583 //120=~7.13 m (partial)
03584 //121=~7.18 m (full)
03585 //122=~7.24 m
03586 if (z<7) flag=1;//planes 1->117 inclusive
03587 else if (z>=7) flag=3;//planes 118 onwards
03588 }
03589 else {
03590 if (z<7) flag=2;
03591 else if (z>=7) flag=4;
03592 }
03593 }
03594 else if (nu.detector==Detector::kFar) {
03595 contained=nu.planeTrkEnd>=planeMin &&
03596 nu.planeTrkEnd<=planeMax &&
03597 r<rMax;
03598 if (contained) flag=1;
03599 else flag=2;
03600 }
03601 else cout<<"Ahh, detector="<<nu.detector<<endl;
03602
03603 return flag;
03604 }
|
|
|
Calc Pittsburgh (Donna/Debdatta) containment flag containmentFlagPitt for ND: 1 = track is fully contained in the upstream region 2 = track is partially contained in the upstream region 3 = track is fully contained in the downstream region 4 = track is partially contained in the downstream region containmentFlagPitt for FD: 1 = track is fully contained 2 = track is partially contained Definition at line 3683 of file NuReco.cxx. References NuEvent::detector, NuEvent::drTrkFidall, NuEvent::dzTrkFidall, MSG, NuEvent::uTrkEnd, NuEvent::vTrkEnd, NuEvent::xTrkEnd, NuEvent::yTrkEnd, and NuEvent::zTrkEnd. Referenced by CalcTrkVariables(). 03684 {
03686
03692
03696
03697 Int_t containmentFlagPitt=0;
03698
03699 const Float_t x=nu.xTrkEnd;
03700 const Float_t y=nu.yTrkEnd;
03701 const Float_t z=nu.zTrkEnd;
03702 const Float_t u=nu.uTrkEnd;
03703 const Float_t v=nu.vTrkEnd;
03704 const Float_t rad=sqrt(x*x + y*y);
03705
03706 //the total number of instrumented planes will be 153 since
03707 //0 is a bookend and 32*4=128 are uninstrumented in the spectrometer
03708
03709 //Different regions in the ND:
03710 //Breakdown of number of planes:
03711 // Veto=21 planes numbered 0-20 (1st is steel bookend)
03712 // Target=40 planes numbered 21-60
03713 // Shower=60 planes numbered 61-120
03714 // Spectrometer=161 planes 121-281
03715 // First and last are instrumented
03716 // 33 have scintillator, 128 are steel
03717 // 4 steel planes for each one with scintillator,
03718 // Altogether: 5*32=160, 160+1 at end=161
03719
03720 //Number of strips:
03721 //Forward section:
03722 //96 planes with 64 strips
03723 //24 planes with 96 strips
03724 //Spectrometer section:
03725 //33 planes with 96 strips
03726 //96+24+33=153 instrumented in total
03727
03728 //strips TPos:
03729 //plane 6 (full) goes -2.64 -> 1.32 m
03730 //plane 11 (full) goes -1.32 -> 2.64 m
03731 //plane 4,8,10 (partial) goes -2.40 -> 0.24 m (V-view)
03732 //plane 5,7,9 (partial) goes -0.24 -> 2.40 m (U-view)
03733 //2.64 - 2.40 = 24 cm = 5.85 strips
03734 //the FI planes "stick out" by ~6 strips
03735
03736 //Tobi's code snippet:
03737 // in the near detector, a further check is needed:
03738 // partial U planes have strips 0-63
03739 // partial V planes have strips 4-67
03740 //if (det==static_cast<Detector::Detector_t>(s.Detector)) {
03741 // if (((pl-1)%5) && (pl%2) && st>63) continue;
03742 // if (((pl-1)%5) && (pl%2)==0 && st<4 ) continue;
03743 //}
03744
03745 //scintillator plane positions:
03746 //first full plane has z=~??? (plane 1)
03747 //last partial plane has z=~7.13 m (plane 120)
03748 //first full spect plane has z=~7.18 m (plane 121)
03749 //last spect plane z=~16.62 m (plane 281)
03750 //117=~6.95 m
03751 //118=~7.01 m
03752 //119=~7.07 m
03753 //120=~7.13 m (partial)
03754 //121=~7.18 m (full)
03755 //122=~7.24 m
03756
03757 if (nu.detector==Detector::kNear) {
03758 Float_t calZ=7*(Munits::m);
03759 Float_t specZ=15.6*(Munits::m);//same as Mad
03760 Float_t minR=0.8*(Munits::m);
03761
03762 Bool_t down_stop=false;
03763 Bool_t down_exit=false;
03764 Bool_t up_stop=false;
03765 Bool_t up_exit=false;
03766
03767 // downstream
03768 if (z>calZ) {
03769 if (u>-0.85 && u<2.1 &&
03770 v>-2.1 && v<0.85 &&
03771 x>-1.5 && x<2.4 &&
03772 y>-1.4 && y<1.4 &&
03773 rad>minR &&
03774 z>calZ && z<=specZ) down_stop=true;
03775 else down_exit=true;
03776 }
03777 //ustream
03778 else if (z<=calZ && z>=0) {
03779 if (u>0.3 && u<1.8 &&
03780 v>-1.8 && v<-0.3 &&
03781 x<2.4 &&
03782 rad>minR &&
03783 z<calZ && z>=0) up_stop=true;
03784 else up_exit=true;
03785 }
03786 else cout<<"Ahhh, z="<<z<<endl;
03787
03788 // set the track's pitt flag
03789 if (down_stop) containmentFlagPitt = 3;
03790 if (up_stop ) containmentFlagPitt = 1;
03791 if (down_exit) containmentFlagPitt = 4;
03792 if (up_exit ) containmentFlagPitt = 2;
03793 }
03794 else if (nu.detector==Detector::kFar) {
03795 //just use CC one for now...
03796 const Float_t rMinFD=+0.5;
03797 Bool_t contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
03798 if (contained) containmentFlagPitt=1;
03799 else containmentFlagPitt=2;
03800 }
03801 else cout<<"Ahhh, detector="<<nu.detector<<endl;
03802
03803 //sanity check
03804 if (!(containmentFlagPitt==1 || containmentFlagPitt==2 ||
03805 containmentFlagPitt==3 || containmentFlagPitt==4)) {
03806 MSG("NuReco",Msg::kWarning)
03807 <<"containmentFlagPitt="<<containmentFlagPitt<<endl;
03808 }
03809
03810 return containmentFlagPitt;
03811 }
|
|
|
Definition at line 3510 of file NuReco.cxx. References NuEvent::containedTrk, NuEvent::detector, and NuEvent::zTrkEnd. Referenced by GetContainmentFlag(). 03511 {
03512 Int_t flag=-1;
03513
03514 //117=~6.95 m
03515 //118=~7.01 m
03516 //119=~7.07 m
03517 //120=~7.13 m (partial)
03518 //121=~7.18 m (full)
03519 if (nu.detector==Detector::kNear) {
03520 if (nu.containedTrk) {
03521 //if (nu.zTrkEnd<7.15) flag=1;
03522 if (nu.zTrkEnd<7.0) flag=1;//make consistent with others
03523 else flag=3;
03524 }
03525 else {
03526 //if (nu.zTrkEnd<7.15) flag=2;
03527 if (nu.zTrkEnd<7.0) flag=2;//make consistent with others
03528 else flag=4;
03529 }
03530 }
03531 else if (nu.detector==Detector::kFar) {
03532 //either contained or not: no spectrometer in FD
03533 if (nu.containedTrk) flag=1;
03534 else flag=2;
03535 }
03536 else cout<<"Ahh, detector="<<nu.detector<<endl;
03537
03538 return flag;
03539 }
|
|
|
Retrieves the cos of the angle between the radial momentum and track vertex.
Definition at line 4666 of file NuReco.cxx. References NtpSRMomentum::best, NtpSRVertex::dcosx, NtpSRVertex::dcosy, NtpSRTrack::momentum, NtpSRTrack::vtx, NtpSRVertex::x, and NtpSRVertex::y. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04667 {
04668 Double_t x = trk.vtx.x,
04669 y = trk.vtx.y,
04670 // Need to reconsider this 'best' momentum - is this the right thing to do?
04671 px = trk.momentum.best * trk.vtx.dcosx,
04672 py = trk.momentum.best * trk.vtx.dcosy;
04673
04674 // Calculate |a||b|
04675 Double_t magxmagp = sqrt((x*x + y*y)*(px*px + py*py));
04676
04677 // Return the cosine value, or zero if the magnitudes are zero
04678 if(magxmagp != 0) return (x*px + y*py)/magxmagp;
04679 else return 0; // yes, this does happen!
04680 }
|
|
|
Code from MadMKAnalysis Definition at line 4394 of file NuReco.cxx. References NuEvent::detector, NuEvent::trkvtxdcosy, and NuEvent::trkvtxdcosz. 04395 {
04397
04398 Float_t dirCosNu=-999;
04399 if (nu.detector==Detector::kFar) {
04400 const Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
04401 const Float_t bl_y = sqrt(1. - bl_z*bl_z);
04402 dirCosNu=nu.trkvtxdcosz*bl_z+nu.trkvtxdcosy*bl_y;
04403 }
04404 else if (nu.detector==Detector::kNear) {
04405 const Float_t nu_cos = -5.799E-2;
04406 const Float_t nu_sin = sqrt(1 -nu_cos*nu_cos);
04407 dirCosNu=nu.trkvtxdcosz*nu_sin + nu.trkvtxdcosy*nu_cos;
04408 }
04409 else cout<<"Ahhh, detector="<<nu.detector<<endl;
04410
04411 //return dirCosNu
04412 return dirCosNu;
04413 }
|
|
||||||||||||
|
Definition at line 4417 of file NuReco.cxx. References NtpSRVertex::dcosy, NtpSRVertex::dcosz, NuEvent::detector, NuEvent::dirCosNu, and NtpSRTrack::vtx. Referenced by CalcKinematicVariables(). 04418 {
04419 // Returns the cosine of the angle between the beam and the track
04420
04421 // Original Code from MadMKAnalysis
04422
04423 // Reworked 19th Aug 2009by Nick Devenish
04424 // Matt Strait and I discovered that the definition of beam
04425 // angle to the detectors is inconsistent between PANS, DSTs and
04426 // others. Matt did some research and got the 'Authoritative'
04427 // answers from Wes Smart at Fermilab. These are now in here.
04428 // Far: 57.184957 millirad
04429 // Near: -58.297760 millirad
04430 // If we get a good, more citable source for these, they should go in
04431
04432 if (nu.detector==Detector::kFar) {
04433 const Float_t bl_y = 0.057184957; // Sin approximation
04434 const Float_t bl_z = sqrt(1 - bl_y*bl_y);
04435 nu.dirCosNu=trk.vtx.dcosz*bl_z+trk.vtx.dcosy*bl_y;
04436 }
04437 else if (nu.detector==Detector::kNear) {
04438 /*
04439 // simple correction based on the vertical position
04440 float vtxy=0;
04441 if(vtx) vtxy=vtx[1]; // in meters
04442 // cosine of the typical neutrino angle in the yz plane
04443 // calculated as py / sqrt( py^2 + pz^2)
04444 float nu_cos = -5.799E-2;
04445 if(vtxy>-2.0 && vtxy<2.0){ //prevents further nuttyness if vtxy is silly
04446 nu_cos += vtxy*1.23304E-3
04447 + vtxy*vtxy*1.08212E-5
04448 + vtxy*vtxy*vtxy*(-4.634E-5);
04449 }
04450 */
04451 const Float_t bl_y = -0.058297760; // Sin approximation
04452 const Float_t bl_z = sqrt(1 - bl_y*bl_y);
04453 nu.dirCosNu=trk.vtx.dcosz*bl_z+trk.vtx.dcosy*bl_y;
04454 }
04455 else cout<<"Ahhh, detector="<<nu.detector<<endl;
04456 }
|
|
||||||||||||
|
Calculates the edge region based on direct variables.
Definition at line 5122 of file NuReco.cxx. References det. 05123 {
05124 /* Description of edge regions:
05125 2
05126 3 ___ 1
05127 /+y \
05128 4 | L+x| 0
05129 \___/
05130 5 7
05131 6
05132 */
05133
05134 // Cannot currently define this for the near detector
05135 if (det == Detector::kNear) return -1;
05136
05137 Float_t positiveangle = phi;
05138 // Force angle to range [0, 2pi] and convert to radians
05139 if(positiveangle < 0) positiveangle += 2*M_PI;
05140 const Float_t torad = M_PI/180;
05141
05142 Float_t theoneang = atan(1.625/3.975); // points at corners
05143 if(positiveangle < theoneang) return 0;
05144 if(positiveangle < 90*torad - theoneang) return 1;
05145 if(positiveangle < 90*torad + theoneang) return 2;
05146 if(positiveangle < 180*torad - theoneang) return 3;
05147 if(positiveangle < 180*torad + theoneang) return 4;
05148 if(positiveangle < 270*torad - theoneang) return 5;
05149 if(positiveangle < 270*torad + theoneang) return 6;
05150 if(positiveangle < 360*torad - theoneang) return 7;
05151
05152 // back around to the bottom half of the first side!
05153 return 0;
05154 }
|
|
||||||||||||
|
Retrieves the FD edge region from the SNTP.
Definition at line 5109 of file NuReco.cxx. References det, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), GetPhi(), and RecHeader::GetVldContext(). Referenced by CalculateEdgeRegion(), NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 05110 {
05111 // Get the detector
05112 Detector::Detector_t det = rec.GetHeader().GetVldContext().GetDetector();
05113 // SimFlag::SimFlag_t simFlag = rec.GetHeader().GetVldContext().GetSimFlag();
05114
05115 // Technically shouldn't recalculate phi here, but no simple and consistent
05116 // way to get the already calculated value without passing in NuEvent, or
05117 return GetEdgeRegion(det, GetPhi(vtx));
05118 }
|
|
||||||||||||
|
Definition at line 1986 of file NuReco.cxx. References NtpStRecord::evt, GetEvtMedianTime(), and MAXMSG. Referenced by NuAnalysis::EnergySpect(), NuAnalysis::EnergySpectMC(), and NuAnalysis::N_1(). 01988 {
01989 Msg::LogLevel_t logLevel=Msg::kDebug;
01990
01991 const TClonesArray& evtTca=(*ntp.evt);
01992 const Int_t numEvts=evtTca.GetEntriesFast();
01993
01994 map<Int_t,Double_t> evtMedianTimes;
01995 multiset<Double_t> allMedianTimes;
01996 MsgFormat ffmt("%9.11f");
01997
01998 //get the median times for all the events
01999 for (Int_t i=0;i<numEvts;i++) {
02000 const NtpSREvent& evt=
02001 *dynamic_cast<NtpSREvent*>(evtTca[i]);
02002 Double_t medTime=this->GetEvtMedianTime(ntp,evt);
02003 evtMedianTimes[i]=medTime;
02004 allMedianTimes.insert(medTime);
02005 }
02006
02007 //now work out the smallest delta time for each event
02008 typedef map<Int_t,Double_t>::iterator evtTimesIt;
02009 for (evtTimesIt it=evtMedianTimes.begin();
02010 it!=evtMedianTimes.end();++it){
02011
02012 //get an iterator to this events time
02013 multiset<Double_t>::iterator tIt=allMedianTimes.find(it->second);
02014 if (tIt!=allMedianTimes.end()){
02015 //look at event time before and after
02016 multiset<Double_t>::iterator tBeforeIt=tIt;
02017 --tBeforeIt;
02018 multiset<Double_t>::iterator tAfterIt=tIt;
02019 ++tAfterIt;
02020
02021 Double_t deltaAfter=-1;
02022 if (tAfterIt!=allMedianTimes.end()){
02023 MAXMSG("NuReco",logLevel,500)
02024 <<"t="<<ffmt(*tIt)
02025 <<", tAfterIt="<<ffmt(*tAfterIt)<<endl;
02026 deltaAfter=fabs((*tAfterIt)-(*tIt));
02027 }
02028
02029 Double_t deltaBefore=-1;
02030 if (tBeforeIt!=allMedianTimes.end()){
02031 MAXMSG("NuReco",logLevel,500)
02032 <<"t="<<ffmt(*tIt)
02033 <<", tBeforeIt="<<ffmt(*tBeforeIt)<<endl;
02034 deltaBefore=fabs((*tBeforeIt)-(*tIt));
02035 }
02036
02037 Double_t smallestDelta=deltaAfter;
02038 if (deltaAfter!=-1 && deltaBefore!=-1){
02039 if (deltaBefore<deltaAfter) smallestDelta=deltaBefore;
02040 }
02041 if (deltaBefore==-1){
02042 smallestDelta=deltaAfter;
02043 }
02044 if (deltaAfter==-1){
02045 smallestDelta=deltaBefore;
02046 }
02047
02048 MAXMSG("NuReco",logLevel,500)
02049 <<"smallestDelta="<<ffmt(smallestDelta)
02050 <<", deltaBefore="<<ffmt(deltaBefore)
02051 <<", deltaAfter="<<ffmt(deltaAfter)
02052 <<endl;
02053
02054 //return everything - no
02055 //deltaTs[it->first]=smallestDelta;
02056 if (smallestDelta>=0) deltaTs[it->first]=smallestDelta;
02057 else cout<<"Ahhh bad delta T, dT="<<smallestDelta<<endl;
02058 }
02059 else cout<<"Ahh, not found time"<<endl;
02060 }
02061 }
|
|
||||||||||||
|
Definition at line 79 of file NuReco.cxx. References NuEvent::anaVersion, CalcEvtVariables(), CalcKinematicVariables(), CalcResolution(), CalcShwVariables(), CalcTrkReclamation(), CalcTrkVariables(), NuEvent::dirCosNu, NuEvent::energy, NuEvent::energyCC, NuEvent::energyNC, NuEvent::energyNoRw, NuEvent::energyRM, NuEvent::evt, MAXMSG, NuEvent::nshw, NuEvent::ntrk, NuEvent::primshw, NuEvent::primtrk, NuEvent::q2, SetBestShw(), SetBestTrk(), NuEvent::shwEn, NuEvent::shwEnCC, NuEvent::shwEnCor1, NuEvent::shwEnCor2, NuEvent::shwEnCor3, NuEvent::shwEnCor4, NuEvent::shwEnCor5, NuEvent::shwEnkNN, NuEvent::shwEnNC, NuEvent::shwExists1, NuEvent::shwExists2, NuEvent::shwExists3, NuEvent::trkEn, NuEvent::trkEnCorCurv1, NuEvent::trkEnCorCurv2, NuEvent::trkEnCorCurv3, NuEvent::trkEnCorRange1, NuEvent::trkEnCorRange2, NuEvent::trkEnCorRange3, NuEvent::trkExists1, NuEvent::trkExists2, NuEvent::trkExists3, NuEvent::w2, NuEvent::x, and NuEvent::y. Referenced by NuDemoModule::Ana(), NuDSTAna::BRevAna(), NuAnalysis::ChargeSeparationOneSnarl(), NuAnalysis::ChargeSignCut(), NuDSTAna::Contamination(), NuDSTAna::DPSystematic(), NuAnalysis::Efficiencies(), NuAnalysis::EnergySpect(), NuAnalysis::EnergySpectMC(), NuDSTAna::FDTestAna(), NuUtilities::FixDogwoodQP(), NuDSTAna::FluxComponents(), NuDSTAna::JeffsTestAna(), NuAnalysis::LIRejectionTest(), NuDSTAna::MakeFCTree(), NuAnalysis::MakeFullDST(), NuDSTAna::MakeMicroDST(), NuDSTAna::MakeMicroDst2010(), NuDSTAna::MakeMicroDstFakeData(), NuDSTAna::MakeMicroDstForCSSSystematics(), NuDSTAna::MakeMicroDstHe(), NuDSTAna::MakeMicroDstJJEPresel(), NuDSTAna::MakeMicroDstWithStdCCRecoAndCuts(), NuDSTAna::MakeSelMicroDST(), NuDSTAna::MMRereco(), NuDSTAna::MMTransition(), NuAnalysis::N_1(), NuDSTAna::NDOsc(), NuDSTAna::NDQPRB(), NuDSTAna::NDTestAna(), NuDSTAna::NewFieldAna(), NuDSTAna::QPStudy(), NuDSTAna::RHCTest(), NuDSTAna::StdCCAna(), and NuDSTAna::StdNMBAna(). 00080 {
00081 Msg::LogLevel_t logLevel=Msg::kDebug;
00082
00083 MAXMSG("NuReco",logLevel,200)
00084 <<"GetEvtEnergy: evt="<<nu.evt<<", trks="<<nu.ntrk
00085 <<", shws="<<nu.nshw<<endl;
00086
00087 //copy trk/shw-1,-2,-3 across to trk/shw
00088 this->SetBestTrk(nu);
00089
00090 this->SetBestShw(nu);
00091
00092 this->CalcEvtVariables(nu);
00093 this->CalcTrkVariables(nu);
00094 this->CalcTrkReclamation(nu);
00095
00096 this->CalcShwVariables(nu, includekNN);
00097
00098 //the final calculations
00099 nu.energyCC=nu.trkEn+nu.shwEnCC;
00100 nu.energyNC=nu.shwEnNC;
00101 nu.energyRM=nu.trkEn;
00102 nu.energy=nu.energyCC;//set as CC by default, analysis user should change this if necessary
00103
00104 // 7.2e20 analysis uses kNN shower energy by default
00105 if(nu.anaVersion == NuCuts::kCC0720Std){
00106 nu.shwEn = nu.shwEnkNN;
00107 nu.energy = nu.trkEn + nu.shwEn;
00108 }
00109
00110 //keep a record of the non-reweighted energy
00111 nu.energyNoRw=nu.energy;
00112 this->CalcKinematicVariables(nu);
00113
00114 //Set Cambridge resolution
00115 this->CalcResolution(nu);
00116
00117 MAXMSG("NuReco",Msg::kInfo,10)
00118 <<endl
00119 <<"GetEvtEnergy:"
00120 <<" nshw="<<nu.nshw
00121 <<", ntrk="<<nu.ntrk
00122 <<", primshw="<<nu.primshw
00123 <<", primtrk="<<nu.primtrk
00124 <<endl
00125 <<"shwExists1,2,3="<<nu.shwExists1<<","
00126 <<nu.shwExists2<<","<<nu.shwExists3
00127 <<", trkExists1,2,3="<<nu.trkExists1<<","
00128 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00129
00130 MAXMSG("NuReco",Msg::kInfo,10)
00131 <<"E="<<nu.energy<<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00132 <<endl
00133 <<"dircosnu="<<nu.dirCosNu
00134 <<", y="<<nu.y<<", q2="<<nu.q2
00135 <<", w2="<<nu.w2<<", x="<<nu.x<<endl;
00136
00137 //sanity check on numbers of trks/shws
00138 if (nu.ntrk>3 || nu.nshw>7) {
00139 MAXMSG("NuReco",Msg::kWarning,3)
00140 <<endl<<"Lots of trks/shws"
00141 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00142 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00143 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00144 <<endl
00145 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00146 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00147 <<", shw5="<<nu.shwEnCor5
00148 //<<", shw6="<<nu.shwEnCor6<<", shw7="<<nu.shwEnCor7
00149 <<endl
00150 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00151 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00152 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00153 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00154 if (nu.ntrk>3) {
00155 MAXMSG("NuReco",Msg::kError,300)
00156 <<endl<<"Lots of trks/shws"
00157 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00158 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00159 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00160 <<endl
00161 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00162 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00163 <<", shw5="<<nu.shwEnCor5
00164 //<<", shw6="<<nu.shwEnCor6<<", shw7="<<nu.shwEnCor7
00165 <<endl
00166 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00167 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00168 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00169 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00170 }
00171 }
00172
00173 return nu.energy;
00174 }
|
|
||||||||||||
|
Definition at line 2845 of file NuReco.cxx. References NtpSRStrip::index, MAXMSG, NtpSREvent::nstrip, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpStRecord::stp, NtpSRStrip::strip, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSRStrip::tpos. Referenced by NuAnalysis::EnergySpectMC(), and GetEvtDeltaTs(). 02847 {
02848 //Msg::LogLevel_t logLevel=Msg::kDebug;
02849
02850 multiset<Double_t> times;
02851 const TClonesArray& stpTca=(*ntp.stp);
02852 //const Int_t numStps=stpTca.GetEntriesFast();
02853
02854 MAXMSG("NuReco",Msg::kDebug,200)
02855 <<"evt.nstrip="<<evt.nstrip<<endl;
02856 for (Int_t i=0;i<evt.nstrip;i++) {
02857 const NtpSRStrip& stp=
02858 *(dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]));
02859
02860 Double_t time=stp.time0;
02861 if (time<0 || time>1) time=stp.time1;
02862
02863 if (time>0 && time<=1) {
02864 times.insert(time);
02865 }
02866 //else just don't put the time in the map - clearly rubbish
02867
02868 MAXMSG("NuReco",Msg::kDebug,500)
02869 <<" Strip index="<<stp.index<<", stp="<<stp.strip
02870 <<", t="<<time
02871 <<", pl="<<stp.plane
02872 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02873 }
02874
02875 //get the median time from the map
02876 multiset<Double_t>::const_iterator it=times.begin();
02877 advance(it,times.size()/2);
02878 Double_t medianTime=*it;
02879 MAXMSG("NuReco",Msg::kDebug,100)
02880 <<"Median time="<<medianTime<<endl;
02881
02882 return medianTime;
02883 }
|
|
||||||||||||
|
Definition at line 1967 of file NuReco.cxx. References NtpStRecord::evt, and NtpSREvent::slc. Referenced by NuAnalysis::EnergySpect(), NuAnalysis::EnergySpectMC(), and NuAnalysis::N_1(). 01969 {
01970 TClonesArray& evtTca=(*ntp.evt);
01971 const Int_t numEvts=evtTca.GetEntriesFast();
01972
01973 //loop over the events
01974 for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
01975 const NtpSREvent* pevt=
01976 dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
01977 const NtpSREvent& evt=(*pevt);
01978
01979 //count up the number of events per slice
01980 evtsPerSlc[evt.slc]++;
01981 }
01982 }
|
|
||||||||||||
|
This is a rewrite of the code in Mad/MadQuantities Definition at line 4349 of file NuReco.cxx. References abs(), NtpMCStdHep::IdHEP, NuEvent::iresonance, NtpMCStdHep::IstHEP, NtpMCStdHep::mc, NuEvent::mc, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 04351 {
04353
04354 Int_t hfs=0;
04355 Int_t proc=nu.iresonance;
04356 TClonesArray* pointStdhepArray = NULL;
04357 pointStdhepArray=ntp.stdhep;
04358 TClonesArray& stdhepArray = *pointStdhepArray;
04359 Int_t nStdHep = stdhepArray.GetEntries();
04360
04361 Int_t itr=nu.mc;
04362
04363 if(proc==1002){
04364 for(int i=0;i<nStdHep;i++){
04365 //LoadStdHep(i);
04366 const NtpMCStdHep* ntpStdHep=
04367 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
04368
04369 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3 &&
04370 !(abs(ntpStdHep->IdHEP)==15)){ //not a tau lepton
04371 hfs = ntpStdHep->IdHEP;
04372 break;
04373 }
04374 }
04375 hfs = hfs%1000;
04376 }
04377 else {
04378 for(int i=0;i<nStdHep;i++){
04379 //LoadStdHep(i);
04380 const NtpMCStdHep* ntpStdHep=
04381 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
04382 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3){
04383 hfs = ntpStdHep->IdHEP;
04384 break;
04385 }
04386 }
04387 hfs = hfs%1000;
04388 }
04389 return hfs;
04390 }
|
|
||||||||||||||||||||
|
Calculates the 'One Horizontal/Vertical strip number'. This number (referred to as HOVE) gives a single number for how far form the edge of the detector a given horizontal or vertical strip. Definition at line 5267 of file NuReco.cxx. References det. Referenced by CalculateHorizontalVerticalStripNumber(). 05268 {
05269 Short_t hoveNum = -999;
05270
05271 // This variable is not defined for the near detector
05272 if (det == Detector::kNear) return hoveNum;
05273
05274 if (edgeRegion == 0) hoveNum = 192 - stripTrku + stripTrkv;
05275 else if(edgeRegion == 2) hoveNum = 383 - stripTrku - stripTrkv;
05276 else if(edgeRegion == 4) hoveNum = 192 + stripTrku - stripTrkv;
05277 else if(edgeRegion == 6) hoveNum = 1 + stripTrku + stripTrkv;
05278
05279 return hoveNum;
05280 }
|
|
||||||||||||
|
This is a rewrite of the code in Mad/MadQuantities Definition at line 3937 of file NuReco.cxx. References abs(), NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, NtpMCStdHep::mc, NuEvent::mc, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 03939 {
03941
03942 //get the index of the true MC interaction
03943 //use the track one for now - have to be carefull if no track
03944 //in event
03945 Int_t itr=nu.mc;
03946
03947 Int_t initial_state=0;
03948 TClonesArray* pointStdhepArray = NULL;
03949 pointStdhepArray=ntp.stdhep;
03950 TClonesArray& stdhepArray = *pointStdhepArray;
03951 Int_t nStdHep = stdhepArray.GetEntries();
03952
03953 int protneut = -1; // 0 = neutron, 1 = proton, 2 = N, 3 = A
03954 int nubarnu = 0; // +1 = neutrino, -1 = antineutrino
03955
03956 for(int i=0;i<nStdHep;i++){
03957 //LoadStdHep(i);
03958 const NtpMCStdHep* ntpStdHep=
03959 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03960 if(ntpStdHep->mc==itr){
03961
03962 if(ntpStdHep->IstHEP==0){ //initial state particle
03963 if(abs(ntpStdHep->IdHEP)==12 ||
03964 abs(ntpStdHep->IdHEP)==14 ||
03965 abs(ntpStdHep->IdHEP)==16){ //(anti)neutrino
03966 nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP); //get sign
03967 }
03968 }
03969 if(ntpStdHep->IstHEP==11){ //target nucleon
03970 if(ntpStdHep->IdHEP==2212) protneut = 1; //proton
03971 else if(ntpStdHep->IdHEP==2112) protneut = 0; //neutron
03972 else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2; //nucleus
03973 else protneut = 3; //atom - probably never get here since IdHEP A>N?
03974 }
03975 }
03976 }
03977
03978 if(protneut==1 && nubarnu==1) initial_state=1; //p + v
03979 if(protneut==0 && nubarnu==1) initial_state=2; //n + v
03980 if(protneut==1 && nubarnu==-1) initial_state=3; //p + vbar
03981 if(protneut==0 && nubarnu==-1) initial_state=4; //n + vbar
03982 if(protneut==2 && nubarnu==1) initial_state=5; //N + v
03983 if(protneut==3 && nubarnu==1) initial_state=6; //A + v
03984 if(protneut==2 && nubarnu==-1) initial_state=7; //N + vbar
03985 if(protneut==3 && nubarnu==-1) initial_state=8; //A + vbar
03986
03987 return initial_state;
03988 }
|
|
|
Definition at line 923 of file NuReco.cxx. References NuEvent::ntrk, NuEvent::trkLength1, NuEvent::trkLength2, and NuEvent::trkLength3. 00924 {
00925 if (nu.ntrk<=0) return -1;
00926 else if (nu.ntrk==1) return 1;
00927 else if (nu.ntrk==2) {
00928 Int_t trkIndex=1;
00929 Int_t longestTrack=nu.trkLength1;
00930 if (nu.trkLength2>longestTrack) {
00931 trkIndex=2;
00932 longestTrack=nu.trkLength2;
00933 }
00934 return trkIndex;
00935 }
00936 else {
00937 Int_t trkIndex=1;
00938 Int_t longestTrack=nu.trkLength1;
00939 if (nu.trkLength2>longestTrack) {
00940 trkIndex=2;
00941 longestTrack=nu.trkLength2;
00942 }
00943 if (nu.trkLength3>longestTrack) {
00944 trkIndex=3;
00945 longestTrack=nu.trkLength3;
00946 }
00947 return trkIndex;
00948 }
00949 }
|
|
||||||||||||
|
Definition at line 953 of file NuReco.cxx. References abs(), NtpSRPlane::beg, NtpSRPlane::end, NtpSRTrack::fit, MAXMSG, NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpSRFitTrack::pass, NtpSRTrack::plane, NtpSREvent::trk, and NtpStRecord::trk. Referenced by GetBestTrack(). 00955 {
00956 TClonesArray& trkTca=(*ntp.trk);
00957 //const Int_t numTrks=trkTca.GetEntriesFast();
00958
00959 MAXMSG("NuReco",Msg::kInfo,1)
00960 <<"Setting best track as longest track (end-beg+1)"<<endl;
00961 Int_t trkToUse=0;
00962 Double_t trkLength=0;
00963 if (evt.ntrack>1){
00964 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00965 const NtpSRTrack& trk=
00966 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00967
00968 Int_t absLength=abs(trk.plane.end-trk.plane.beg+1);
00969
00970 //select the longest trk
00971 if (absLength>trkLength) {
00972 trkLength=absLength;
00973 trkToUse=itrk;
00974 }
00975
00976 MAXMSG("NuReco",Msg::kDebug,200)
00977 <<"itrk="<<itrk<<" of "<<evt.ntrack
00978 <<", trkToUse="<<trkToUse
00979 <<", pass="<<trk.fit.pass
00980 <<", ntrklike="<<trk.plane.ntrklike
00981 <<", absLength="<<absLength
00982 <<", longest="<<trkLength
00983 <<endl;
00984 }
00985 MAXMSG("NuReco",Msg::kDebug,200)
00986 <<"Finished loop, selected trkToUse="<<trkToUse<<endl;
00987 const NtpSRTrack& trk=
00988 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00989 return &trk;
00990 }
00991 else if (evt.ntrack==1){
00992 const NtpSRTrack& trk=
00993 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00994 return &trk;
00995 }
00996 else return 0;
00997 }
|
|
||||||||||||
|
Definition at line 881 of file NuReco.cxx. References NtpSRTrack::fit, MAXMSG, NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpSRFitTrack::pass, NtpSRTrack::plane, NtpSREvent::trk, and NtpStRecord::trk. 00882 {
00883 TClonesArray& trkTca=(*ntp.trk);
00884 //const Int_t numTrks=trkTca.GetEntriesFast();
00885
00886 MAXMSG("NuReco",Msg::kInfo,1)
00887 <<"Setting best track as longest track (track-like planes)"<<endl;
00888 Int_t trkToUse=0;
00889 Double_t trkLength=0;
00890 if (evt.ntrack>1){
00891 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00892 const NtpSRTrack& trk=
00893 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00894
00895 MAXMSG("NuReco",Msg::kDebug,200)
00896 <<"trk="<<itrk<<"/"<<evt.ntrack
00897 <<", pass="<<trk.fit.pass
00898 <<", ntrklike="<<trk.plane.ntrklike
00899 <<endl;
00900
00901 //select the longest trk (in track-like planes)
00902 if (trk.plane.ntrklike>trkLength) {
00903 trkLength=trk.plane.ntrklike;
00904 trkToUse=itrk;
00905 }
00906 }
00907 MAXMSG("NuReco",Msg::kDebug,200)
00908 <<"selected trkToUse="<<trkToUse<<endl;
00909 const NtpSRTrack& trk=
00910 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00911 return &trk;
00912 }
00913 else if (evt.ntrack==1){
00914 const NtpSRTrack& trk=
00915 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00916 return &trk;
00917 }
00918 else return 0;
00919 }
|
|
||||||||||||
|
Experimental: Attempts to calculate the 'containment truth'.
Definition at line 5032 of file NuReco.cxx. References NtpMCStdHep::child, GetPrimaryMuonStdHepIndex(), NtpMCStdHep::IdHEP, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 05033 {
05034 Int_t muonindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
05035 // If there is no muon, then we have no containment....
05036 if (muonindex == -1) return false;
05037
05038 // Pull out the TCA from the sntp record
05039 TClonesArray& hepTca=(*ntp.stdhep);
05040
05041 // Grab ahold of this muon
05042 NtpMCStdHep* muon = dynamic_cast<NtpMCStdHep*>(hepTca[muonindex]);
05043
05044 // Trace down the children tree. Since it appears that exiting muons do not
05045 // have their decays simulated in the MC, check for decay products and if
05046 // there are none, then assume that it decayed outside the detector i.e.
05047 // it exited
05048
05049 // What is our calculated containment?
05050 Bool_t containment = false;
05051
05052 // Loop until we have a result
05053 while (true) {
05054 // Pointer to the muon 'child'
05055 NtpMCStdHep* child = 0;
05056
05057 // If we have children....
05058 if (muon->child[0] != -1) {
05059
05060 // Loop over all children (may just be one)
05061 for (Int_t i = muon->child[0]; i <= muon->child[1]; ++i) {
05062 // Grab the child object
05063 NtpMCStdHep* childcand = dynamic_cast<NtpMCStdHep*>(hepTca[i]);
05064
05065 // If this child is a muon....
05066 if (childcand->IdHEP == 13 || childcand->IdHEP == -13) {
05067 // there is a muon child!
05068 child = childcand;
05069 }
05070 } // Loop over children
05071
05072 // We have now looped over all children. If none of them were muons,
05073 // then assume that this is because the muon has decayed in the detector!
05074 if (child == 0) {
05075 // No muon children found. Must have decayed
05076 containment = true;
05077 break;
05078 } else {
05079 // We have a muon child. This is now the muon for the next loop.
05080 muon = child;
05081 }
05082 } else {
05083 // We have no children. This must mean we have exited the detector!
05084 containment = false;
05085 break;
05086 }
05087 } // End of infinite loop
05088
05089 // We should have now resolved the containment.
05090 /* this->PrintStdhep(ntp, mc);
05091 if (containment) {
05092 MSG("NuReco", Msg::kInfo) << "==== CONTAINED ====" << endl;
05093 } else {
05094 MSG("NuReco", Msg::kInfo) << "~~~~~ ESCAPED ~~~~~" << endl;
05095 }*/
05096 // Finally, return.
05097 return containment;
05098 }
|
|
||||||||||||
|
Retrieves the energy of muon on the first detector hit.
Definition at line 4993 of file NuReco.cxx. References NtpMCStdHep::dethit, GetPrimaryMuonStdHepIndex(), NtpMCStdHepHit::mom, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 04994 {
04995 Int_t stdindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
04996 if (stdindex == -1) return -1;
04997
04998 // Get the muon hit object
04999 TClonesArray& hepTca=(*ntp.stdhep);
05000 const NtpMCStdHep& muon = *dynamic_cast<NtpMCStdHep*>(hepTca[stdindex]);
05001
05002 // return the energy
05003 return muon.dethit[0].mom[3];
05004 }
|
|
||||||||||||
|
Retrieves the energy of muon on the last detector hit.
Definition at line 5008 of file NuReco.cxx. References NtpMCStdHep::dethit, GetPrimaryMuonStdHepIndex(), NtpMCStdHepHit::mom, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 05009 {
05010 Int_t stdindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
05011 if (stdindex == -1) return -1;
05012
05013 // Get the muon hit object
05014 TClonesArray& hepTca=(*ntp.stdhep);
05015 const NtpMCStdHep& muon = *dynamic_cast<NtpMCStdHep*>(hepTca[stdindex]);
05016
05017 // return the energy
05018 return muon.dethit[1].mom[3];
05019 }
|
|
||||||||||||
|
This is a rewrite of the code in Mad/MadQuantities Definition at line 3992 of file NuReco.cxx. References NuEvent::aMC, and NuEvent::zMC. Referenced by GetTruthInfo(). 03994 {
03996
03997
03998 Int_t z=nu.zMC;//int(ntpTruth->z);
03999 Int_t a=nu.aMC;//int(ntpTruth->a);
04000 Int_t nucleus=1;
04001
04002 switch (z) {
04003 // case 1:
04004 //nucleus=0; // free nucleon
04005 //break;
04006 case 1:
04007 switch (a) {
04008 case 1:
04009 nucleus=256; // hydrogen1
04010 break;
04011 case 2:
04012 nucleus=257; // hydrogen2
04013 break;
04014 case 3:
04015 nucleus=258; // hydrogen2
04016 break;
04017 default:
04018 nucleus=256; // hydrogen1
04019 break;
04020 }
04021 break;
04022 case 6:
04023 switch (a) {
04024 case 11:
04025 nucleus=274; // carbon11
04026 break;
04027 case 12:
04028 nucleus=275; // carbon12
04029 break;
04030 case 13:
04031 nucleus=276; // carbon13
04032 break;
04033 case 14:
04034 nucleus=277; // carbon14
04035 break;
04036 default:
04037 nucleus=275; // carbon12
04038 break;
04039 }
04040 break;
04041 case 7:
04042 switch (a) {
04043 case 13:
04044 nucleus=278; // nitrogen13
04045 break;
04046 case 14:
04047 nucleus=279; // nitrogen14
04048 break;
04049 case 15:
04050 nucleus=280; // nitrogen15
04051 break;
04052 case 16:
04053 nucleus=281; // nitrogen16
04054 break;
04055 case 17:
04056 nucleus=282; // nitrogen17
04057 break;
04058 default:
04059 nucleus=279; // nitrogen14
04060 break;
04061 }
04062 break;
04063 case 8:
04064 switch (a) {
04065 case 15:
04066 nucleus=283; // oxygen15
04067 break;
04068 case 16:
04069 nucleus=284; // oxygen16
04070 break;
04071 case 17:
04072 nucleus=285; // oxygen17
04073 break;
04074 case 18:
04075 nucleus=286; // oxygen18
04076 break;
04077 default:
04078 nucleus=284; // oxygen16
04079 break;
04080 }
04081 break;
04082 case 13:
04083 switch (a) {
04084 case 26:
04085 nucleus=303; // aluminium26
04086 break;
04087 case 27:
04088 nucleus=304; // aluminium27
04089 break;
04090 case 28:
04091 nucleus=305; // aluminium28
04092 break;
04093 case 29:
04094 nucleus=306; // aluminium29
04095 break;
04096 default:
04097 nucleus=304; // aluminium27
04098 break;
04099 }
04100 break;
04101 case 14:
04102 switch (a) {
04103 case 27:
04104 nucleus=307; // silicon27
04105 break;
04106 case 28:
04107 nucleus=308; // silicon28
04108 break;
04109 case 29:
04110 nucleus=309; // silicon29
04111 break;
04112 case 30:
04113 nucleus=310; // silicon30
04114 break;
04115 default:
04116 nucleus=308; // silicon28
04117 break;
04118 }
04119 break;
04120 case 15:
04121 switch (a) {
04122 case 30:
04123 nucleus=311; //phosphorus30
04124 break;
04125 case 31:
04126 nucleus=312; //phosphorus31
04127 break;
04128 case 32:
04129 nucleus=313; //phosphorus32
04130 break;
04131 case 33:
04132 nucleus=314; //phosphorus33
04133 break;
04134 default:
04135 nucleus=312; //phosphorus31
04136 break;
04137 }
04138 break;
04139 case 16:
04140 switch (a) {
04141 case 31:
04142 nucleus=315; //sulphur31
04143 break;
04144 case 32:
04145 nucleus=316; //sulphur32
04146 break;
04147 case 33:
04148 nucleus=317; //sulphur33
04149 break;
04150 case 34:
04151 nucleus=318; //sulphur34
04152 break;
04153 case 35:
04154 nucleus=319; //sulphur35
04155 break;
04156 case 36:
04157 nucleus=320; //sulphur36
04158 break;
04159 default:
04160 nucleus=316; //sulphur32
04161 break;
04162 }
04163 break;
04164 case 22:
04165 switch (a) {
04166 case 45:
04167 nucleus=347; //titanium45
04168 break;
04169 case 46:
04170 nucleus=348; //titanium46
04171 break;
04172 case 47:
04173 nucleus=349; //titanium47
04174 break;
04175 case 48:
04176 nucleus=350; //titanium48
04177 break;
04178 case 49:
04179 nucleus=351; //titanium49
04180 break;
04181 case 50:
04182 nucleus=352; //titanium50
04183 break;
04184 default:
04185 nucleus=350; //titanium48
04186 break;
04187 }
04188 break;
04189 case 23:
04190 switch (a) {
04191 case 49:
04192 nucleus=353; //vanadium49
04193 break;
04194 case 50:
04195 nucleus=354; //vanadium50
04196 break;
04197 case 51:
04198 nucleus=355; //vanadium51
04199 break;
04200 case 52:
04201 nucleus=356; //vanadium52
04202 break;
04203 case 53:
04204 nucleus=357; //vanadium53
04205 break;
04206 default:
04207 nucleus=355; //vanadium51
04208 break;
04209 }
04210 break;
04211 case 24:
04212 switch (a) {
04213 case 49:
04214 nucleus=358; //chromium49
04215 break;
04216 case 50:
04217 nucleus=359; //chromium50
04218 break;
04219 case 51:
04220 nucleus=360; //chromium51
04221 break;
04222 case 52:
04223 nucleus=361; //chromium52
04224 break;
04225 case 53:
04226 nucleus=362; //chromium53
04227 break;
04228 case 54:
04229 nucleus=363; //chromium54
04230 break;
04231 default:
04232 nucleus=361; //chromium52
04233 break;
04234 }
04235 break;
04236 case 25:
04237 switch (a) {
04238 case 53:
04239 nucleus=364; //manganese53
04240 break;
04241 case 54:
04242 nucleus=365; //manganese54
04243 break;
04244 case 55:
04245 nucleus=366; //manganese55
04246 break;
04247 case 56:
04248 nucleus=367; //manganese56
04249 break;
04250 case 57:
04251 nucleus=368; //manganese57
04252 break;
04253 default:
04254 nucleus=366; //manganese55
04255 break;
04256 }
04257 break;
04258 case 26:
04259 switch (a) {
04260 case 53:
04261 nucleus=369; //iron53
04262 break;
04263 case 54:
04264 nucleus=370; //iron54
04265 break;
04266 case 55:
04267 nucleus=371; //iron55
04268 break;
04269 case 56:
04270 nucleus=372; //iron56
04271 break;
04272 case 57:
04273 nucleus=373; //iron57
04274 break;
04275 case 58:
04276 nucleus=374; //iron58
04277 break;
04278 default:
04279 nucleus=372; //iron56
04280 break;
04281 }
04282 break;
04283 case 28:
04284 switch (a) {
04285 case 57:
04286 nucleus=382; //nickel57
04287 break;
04288 case 58:
04289 nucleus=383; //nickel58
04290 break;
04291 case 59:
04292 nucleus=384; //nickel59
04293 break;
04294 case 60:
04295 nucleus=385; //nickel60
04296 break;
04297 case 61:
04298 nucleus=386; //nickel61
04299 break;
04300 case 62:
04301 nucleus=387; //nickel62
04302 break;
04303 case 63:
04304 nucleus=388; //nickel63
04305 break;
04306 case 64:
04307 nucleus=389; //nickel64
04308 break;
04309 default:
04310 nucleus=383; //nickel58
04311 break;
04312 }
04313 break;
04314 case 29:
04315 switch (a) {
04316 case 62:
04317 nucleus=390; //copper62
04318 break;
04319 case 63:
04320 nucleus=391; //copper63
04321 break;
04322 case 64:
04323 nucleus=392; //copper64
04324 break;
04325 case 65:
04326 nucleus=393; //copper65
04327 break;
04328 case 66:
04329 nucleus=394; //copper66
04330 break;
04331 case 67:
04332 nucleus=395; //copper67
04333 break;
04334 default:
04335 nucleus=392; //copper64
04336 break;
04337 }
04338 break;
04339 default:
04340 nucleus=1; // unknown
04341 break;
04342 }
04343
04344 return nucleus;
04345 }
|
|
||||||||||||||||||||
|
How many strips in, on a diagonal edge, is a hit?
Definition at line 5182 of file NuReco.cxx. References det. Referenced by CalculateParallelStripInset(). 05183 {
05184 Short_t parallelStripVtx = -10;
05185
05186 // Don't calculate this value for the near detector, for now
05187 if (det == Detector::kNear) return parallelStripVtx;
05188
05189 if (edgeRegion == 1) parallelStripVtx = 191-stripBegU;
05190 else if(edgeRegion == 3) parallelStripVtx = 191-stripBegV;
05191 else if(edgeRegion == 5) parallelStripVtx = stripBegU;
05192 else if(edgeRegion == 7) parallelStripVtx = stripBegV;
05193 else parallelStripVtx = -10;
05194
05195 return parallelStripVtx;
05196 }
|
|
||||||||||||||||
|
Is the hit both on a diagonal edge and on the overhanging (perpendicular) strip?
Definition at line 5224 of file NuReco.cxx. References det. Referenced by CalculatePerpendicularStripFlag(). 05225 {
05226 // Convenience variables
05227 const Bool_t Isu = IsHitu;
05228 const Bool_t Isv = !Isu;
05229
05230 Char_t flag = -10;
05231
05232 // We don't have this defined for the near detector
05233 if (det == Detector::kNear) return flag;
05234
05235 if(edgeRegion == 1 || edgeRegion == 5) flag = Isv;
05236 else if(edgeRegion == 3 || edgeRegion == 7) flag = Isu;
05237
05238 return flag;
05239 }
|
|
|
Definition at line 5102 of file NuReco.cxx. References NtpSRVertex::x, and NtpSRVertex::y. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), NuExtraction::ExtractThirdTrkInfo(), and GetEdgeRegion().
|
|
||||||||||||
|
Definition at line 3867 of file NuReco.cxx. References NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, NtpMCStdHep::mass, NuEvent::maxMomPreInukeFSneutMC, NuEvent::maxMomPreInukeFSprotMC, NtpMCStdHep::mc, NuEvent::mc, NuEvent::numPreInukeFSneutMC, NuEvent::numPreInukeFSprotMC, NtpMCStdHep::p4, and NtpStRecord::stdhep. Referenced by GetTruthInfo(). 03869 {
03870 Int_t itr=nu.mc;
03871
03872 TClonesArray* pointStdhepArray = NULL;
03873 pointStdhepArray=ntp.stdhep;
03874 TClonesArray& stdhepArray = *pointStdhepArray;
03875 Int_t nStdHep = stdhepArray.GetEntries();
03876
03877 vector<Int_t> vChildIndices;
03878
03879 for(int i=0;i<nStdHep;i++){
03880 //LoadStdHep(i);
03881 const NtpMCStdHep* ntpStdHep=
03882 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03883 if(ntpStdHep->mc==itr){
03884
03885 if (ntpStdHep->IstHEP==3){
03886 Int_t childStart = ntpStdHep->child[0];
03887 Int_t childEnd = ntpStdHep->child[1];
03888 for (Int_t goodChild = childStart; goodChild<=childEnd; ++goodChild){
03889 vChildIndices.push_back(goodChild);
03890 }
03891 }
03892 }
03893 }
03894
03895 Int_t numPreInukeFSprot = 0;
03896 Int_t numPreInukeFSneut = 0;
03897 Float_t maxMomPreInukeFSprot = 0.;
03898 Float_t maxMomPreInukeFSneut = 0.;
03899
03900 for (UInt_t i=0;i<vChildIndices.size();i++){
03901 //LoadStdHep(i);
03902 const NtpMCStdHep* ntpStdHep=
03903 dynamic_cast<NtpMCStdHep*>(stdhepArray[vChildIndices[i]]);
03904 if(ntpStdHep->mc==itr){
03905 Float_t energy = ntpStdHep->p4[3];
03906 Float_t mass = ntpStdHep->mass;
03907 Float_t mom = 0.0;
03908 if (energy*energy - mass*mass > 0.0){
03909 mom = TMath::Sqrt(energy*energy - mass*mass);
03910 }
03911 if (ntpStdHep->IdHEP==2212){
03912 //I am a proton
03913 numPreInukeFSprot++;
03914 if (mom>maxMomPreInukeFSprot){
03915 maxMomPreInukeFSprot = mom;
03916 }
03917 }
03918 else if (ntpStdHep->IdHEP==2112){
03919 //I am a neutron
03920 numPreInukeFSneut++;
03921 if (mom>maxMomPreInukeFSneut){
03922 maxMomPreInukeFSneut = mom;
03923 }
03924 }
03925 }
03926 }
03927
03928 nu.numPreInukeFSprotMC = numPreInukeFSprot;
03929 nu.numPreInukeFSneutMC = numPreInukeFSneut;
03930 nu.maxMomPreInukeFSprotMC = maxMomPreInukeFSprot;
03931 nu.maxMomPreInukeFSneutMC = maxMomPreInukeFSneut;
03932
03933 }
|
|
||||||||||||
|
Find the stdhep array entry that corresponds to the muon.
Definition at line 4917 of file NuReco.cxx. References NtpMCStdHep::child, NtpMCTruth::iaction, NtpMCStdHep::IdHEP, NtpMCTruth::inu, MSG, PrintStdhep(), NtpMCTruth::stdhep, and NtpStRecord::stdhep. Referenced by GetMuonContainmentMC(), GetMuonFirstHitEnergy(), and GetMuonLastHitEnergy(). 04918 {
04919 // Ignore NC interactions
04920 if (mc.iaction == 0) return -1;
04921
04922 // Ignore non-muon interactions
04923 if (mc.inu != 14 && mc.inu != -14) return -1;
04924
04925 // Pull out the TCA from the sntp record
04926 TClonesArray& hepTca=(*ntp.stdhep);
04927
04928 // Read the first item in the array to find the incident neutrino
04929 const NtpMCStdHep& neutrino = *dynamic_cast<NtpMCStdHep*>(hepTca[mc.stdhep[0]]);
04930 // Read the first child into the array. All relevant numus should have only one
04931 const NtpMCStdHep& child = *dynamic_cast<NtpMCStdHep*>(hepTca[neutrino.child[0]]);
04932
04933 // Ensure we only had one child
04934 Int_t childcount = neutrino.child[1] - neutrino.child[0];
04935 if (childcount > 0) {
04936 // Display debugging information and die
04937 this->PrintStdhep(ntp, mc);
04938 MSG("NuReco",Msg::kFatal) << "More than one primary child in interaction" << endl;
04939 return 0;
04940 }
04941
04942 // // // Print out the array
04943 // this->PrintStdhep(ntp, mc);
04944 // //
04945 // // // Print out the incident particle information
04946 // MSG("NuReco",Msg::kInfo) << "Incident particle: " << neutrino.IdHEP << endl;
04947 // MSG("NuReco",Msg::kInfo) << "Primary Child: " << child.IdHEP << endl;
04948
04949 // If the child is not a muon
04950 if (child.IdHEP != 13 && child.IdHEP != -13) {
04951 MSG("NuReco", Msg::kFatal) << "Not a muon child" << endl;
04952 }
04953
04954 // We believe this child is the muon
04955 return neutrino.child[0];
04956
04957 // // At this point, we know that that child is an only child, and a muon.
04958 // // This is good enough to decide we have found our primary muon
04959 //
04960 // // Print out the momentum 4-vector, to ensure we know which is the energy
04961 // // ->dethit[1].mom[3];
04962 // MSG("NuReco", Msg::kInfo) << "Momentum at start: "
04963 // << child.dethit[0].mom[0] << ", "
04964 // << child.dethit[0].mom[1] << ", "
04965 // << child.dethit[0].mom[2] << ", "
04966 // << child.dethit[0].mom[3] << endl;
04967 // MSG("NuReco", Msg::kInfo) << "Momentum at end: "
04968 // << child.dethit[1].mom[0] << ", "
04969 // << child.dethit[1].mom[1] << ", "
04970 // << child.dethit[1].mom[2] << ", "
04971 // << child.dethit[1].mom[3] << endl;
04972 //
04973 // MSG("NuReco", Msg::kInfo) << "Muon Vertex: "
04974 // << child.vtx[0] << ", "
04975 // << child.vtx[1] << ", "
04976 // << child.vtx[2] << ", "
04977 // << child.vtx[3] << endl;
04978 //
04979 // const NtpMCStdHep& childer = *dynamic_cast<NtpMCStdHep*>(hepTca[child.child[0]]);
04980 // MSG("NuReco", Msg::kInfo) << "Muon's Childs Vertex: "
04981 // << childer.vtx[0] << ", "
04982 // << childer.vtx[1] << ", "
04983 // << childer.vtx[2] << ", "
04984 // << childer.vtx[3] << endl;
04985 //
04986 //
04987 // cout << endl << endl;
04988 // return 0;
04989 }
|
|
|
if a track exists this algorithm relies on the best track having been set - so must set best track first Definition at line 1760 of file NuReco.cxx. References MAXMSG, NuEvent::nshw, NuEvent::primshw, NuEvent::shwEnCor1, NuEvent::trkExists, NuEvent::zShwVtx1, and NuEvent::zTrkVtx. Referenced by GetBestShower(). 01761 {
01764 if (nu.nshw<=0) {
01765 MAXMSG("NuReco",Msg::kDebug,3)
01766 <<"GetPrimaryShower: No showers so return 0"<<endl;
01767 return 0;
01768 }
01769
01770 if (nu.primshw<0) {
01771 // CC only uses a shower if primshw is set, so return here
01772 return 0;
01773 }
01774 //if no track then just return the shw
01775 //doesn't matter about proximity to a non-existant track
01776 if (!nu.trkExists) return 1;
01777
01778 MAXMSG("NuReco",Msg::kInfo,1)
01779 <<"GetPrimaryShw: require trk vtx within 0.5 m || >2 GeV"<<endl;
01780
01781 //in the case of 1 shower you don't know if it
01782 //is actually declared as the primary shower in Birch
01783 //use vtx1 for shower since don't know if primary yet
01784 if (fabs(nu.zTrkVtx-nu.zShwVtx1)<(0.5*Munits::m) ||
01785 (fabs(nu.zTrkVtx-nu.zShwVtx1)>=(0.5*Munits::m) &&
01786 nu.shwEnCor1>(2*Munits::GeV))) {
01787 return 1;
01788 }
01789 else {
01790 return 0;
01791 }
01792 }
|
|
|
Definition at line 1796 of file NuReco.cxx. References MAXMSG, NuEvent::nshw, NuEvent::primshw, NuEvent::shwExists1, NuEvent::shwExists2, and NuEvent::shwExists3. Referenced by GetBestShower(). 01797 {
01798 if (nu.nshw<=0) {
01799 MAXMSG("NuReco",Msg::kInfo,3)
01800 <<"GetPrimaryShower: No showers so return 0"<<endl;
01801 return 0;
01802 }
01803
01804 MAXMSG("NuReco",Msg::kInfo,1)
01805 <<"GetPrimaryShw: using evt.primshw"<<endl;
01806 if (nu.primshw<0) return 0;
01807 else {
01808 if (nu.primshw==0 && nu.shwExists1) return 1;
01809 else if (nu.primshw==1 && nu.shwExists2) {
01810 MAXMSG("NuReco",Msg::kWarning,300)
01811 <<"nu.primshw="<<nu.primshw<<endl;
01812 return 2;
01813 }
01814 else if (nu.primshw==2 && nu.shwExists3) {
01815 MAXMSG("NuReco",Msg::kWarning,300)
01816 <<"nu.primshw="<<nu.primshw<<endl;
01817 return 3;
01818 }
01819 else {
01820 cout<<"ahhhhhhhh nu.primshw="<<nu.primshw<<endl;
01821 return 0;
01822 }
01823 }
01824 }
|
|
|
Definition at line 1752 of file NuReco.cxx. References NuEvent::trkExists1. Referenced by GetBestTrack(). 01753 {
01754 if (nu.trkExists1) return 1;
01755 else return 0;
01756 }
|
|
||||||||||||||||
|
Definition at line 2950 of file NuReco.cxx. References ChooseTruthIndexToUse(), GetShowerWithIndexX(), GetTrackWithIndexX(), NtpSRTrack::index, NtpSRShower::index, NtpSREvent::index, MAXMSG, NuEvent::mcEvt, NuEvent::mcShw, NuEvent::mcTrk, NtpTHTrack::neumc, NtpTHShower::neumc, NtpTHEvent::neumc, NtpStRecord::thevt, NtpStRecord::thshw, and NtpStRecord::thtrk. Referenced by GetTruthInfo(). 02952 {
02953 TClonesArray& thevtTca=(*ntp.thevt);//TruthHelper Event TCA
02954 const Int_t numthevts=thevtTca.GetEntriesFast();
02955 if (numthevts<=0) {
02956 MAXMSG("NuReco",Msg::kDebug,1)
02957 <<"No THEvents, so can't GetPrimaryTruthIndex..."<<endl;
02958 return;
02959 }
02960 MAXMSG("NuReco",Msg::kDebug,1)
02961 <<"Found THEvent, GetPrimaryTruthIndex..."<<endl;
02962
02963 Int_t mcTrk=-1;//track mc index
02964 Int_t mcShw=-1;//shower mc index
02965 Int_t mcEvt=-1;//event mc index
02966
02967 const NtpTHEvent& thevt=
02968 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02969 mcEvt=thevt.neumc;
02970
02971 //get an instance of the code library
02972 //const NuLibrary& lib=NuLibrary::Instance();
02973
02974 //get the shower mc
02975 TClonesArray& thshwTca=(*ntp.thshw);//TruthHelper Shower TCA
02976 const Int_t numthshws=thshwTca.GetEntriesFast();
02977 if (numthshws>0) {
02978 Int_t bestShower=1;//primary shower
02979 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02980 bestShower-1);
02981 if (pshw) {
02982 const NtpSRShower& shw=*pshw;
02983 const NtpTHShower& thshw=
02984 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02985 mcShw=thshw.neumc;
02986 }
02987 }
02988
02989 //get the track mc
02990 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Track TCA
02991 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02992 if (numthtrks>0) {
02993 Int_t bestTrack=1;//primary trk
02994 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02995 bestTrack-1);
02996 if (ptrk) {
02997 const NtpSRTrack& trk=*ptrk;
02998 const NtpTHTrack& thtrk=
02999 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03000 mcTrk=thtrk.neumc;
03001 }
03002 }
03003
03004 //set the variables in the object
03005 nu.mcTrk=mcTrk;
03006 nu.mcShw=mcShw;
03007 nu.mcEvt=mcEvt;
03008 this->ChooseTruthIndexToUse(nu);
03009 }
|
|
||||||||||||||||
|
Definition at line 4460 of file NuReco.cxx. References NuEvent::anaVersion, and NuEvent::detector. Referenced by CalcEvtVariables(), CalcExtraTruthVariables(), CalcTrkVariables(), and NuPlots::FillContainmentHistos(). 04461 {
04462 Float_t zerox=0;
04463 Float_t zeroy=0;
04464 if (nu.detector==Detector::kNear) {
04465 //use beam spot centre
04466 if (nu.anaVersion==NuCuts::kCC0093Std) {
04467 zerox=1.4885*(Munits::m);
04468 zeroy=0.1397*(Munits::m);
04469 }
04470 else {
04471 zerox=1.4828*(Munits::m);//new for 2.5 analysis
04472 zeroy=0.2384*(Munits::m);//new for 2.5 analysis
04473 }
04474 }
04475 //return the radius
04476 return sqrt(pow((x-zerox),2)+pow((y-zeroy),2));
04477 }
|
|
||||||||||||||||||||||||||||
|
Calculates the detector region from direct variables. Currently needs plane - shouldn't this be calculable from the z? Definition at line 4793 of file NuReco.cxx. References det, FidVol::getEvtVtxZOffset(), FidVol::getFarRinner(), FidVol::getFarRouter(), FidVol::getFarZData(), FidVol::getFarZMC(), and FidVol::getTrkVtxZOffset(). 04794 {
04795 // Region of the detector where this track vertex is.
04796 // (Do not call this for shower or event vertices)
04797 //
04798 // -1 if there is no track
04799 // -2 if we don't know how to answer (i.e. this isn't the far detector)
04800 //
04801 // If there is a track, then for the highest energy one's vertex, returns
04802 // the sum of:
04803 //
04804 // 1 if it's near the edge
04805 // 2 if it's near the coil hole
04806 //
04807 // 4 if it's near the front
04808 // 8 if it's near the south side of the supermodule gap
04809 // 12 if it's near the north side of the supermodule gap
04810 // 16 if it's near the back
04811 //
04812 // Or to put it another way, it uses the bit format:
04813 //
04814 // 0 1 2 3
04815 // 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
04816 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
04817 // | reserved | L | T |
04818 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
04819 //
04820 // T: Transverse position:
04821 // 0 if transversely fiducial
04822 // 1 if near the edge
04823 // 2 if near the coil hole
04824 //
04825 // L: Longitudinal position:
04826 // 0 if longitudinally fiducial
04827 // 1 if near the front
04828 // 2 if near the south side of the supermodule gap
04829 // 3 if near the north side of the supermodule gap
04830 // 4 if near the back
04831 //
04832 // Exhaustively listing the possibilities:
04833 //
04834 // 0: Fiducial
04835 // 1: Near the edge and not near any corner
04836 // 2: Near the coil hole and not near any corner
04837 //
04838 // 4: Near the front and not near any corner
04839 // 5: Near the front and the edge
04840 // 6: Near the front and the coil hole
04841 //
04842 // 8: Near the south side of the supermodule gap and not near any corner
04843 // 9: Near the south side of the supermodule gap and the edge
04844 // 10: Near the south side of the supermodule gap and the coil hole
04845 //
04846 // 12: Near the north side of the supermodule gap and not near any corner
04847 // 13: Near the north side of the supermodule gap and the edge
04848 // 14: Near the north side of the supermodule gap and the coil hole
04849 //
04850 // 16: Near the back and not near any corner
04851 // 17: Near the back and the edge
04852 // 18: Near the back and the coil hole
04853 // An enum for the rock detector region - this probably shouldn't go here,
04854 // but for now it will have to do (would rather use an enum than plain numbers
04855 enum Rock_DetectorRegion {
04856 kUnknown = -2,
04857 kNoTrack = -1,
04858 kFiducial = 0,
04859 kEdge = 1, // radial
04860 kCoilHole = 2, // radial
04861 kFront = 4, // longidudinal
04862 kGapSouth = 8, // longidudinal
04863 kGapNorth = 12, // longidudinal
04864 kBack = 16 // longidudinal
04865 };
04866
04867 // Placeholder to avoid nesting if's in this function
04868 // Bool_t inFid = false; // Don't use this variable now
04869 Bool_t nearFront = false, nearEdge = false;
04870 Bool_t nearGapSouth = false, nearGapNorth = false;
04871 Bool_t nearCoil = false, nearBack = false;
04872 Double_t radius2 = vtxX*vtxX + vtxY*vtxY;
04873 Double_t z = vtxZ - FidVol::getTrkVtxZOffset();
04874
04875 // If we are the near detector, we can't calculate the region - return unknown
04876 if (det != Detector::kFar) return kUnknown;
04877
04878 // Get the fiducial volume Z boundaries correctly.
04879 // Read them out into a placeholder array, as this is easier than
04880 // switching every time we call it, and clearer than using a function
04881 // pointer to do the same thing.
04882 Double_t getFarZ[4] = {};
04883 if (simFlag == SimFlag::kData) {
04884 getFarZ[0] = FidVol::getFarZData(0);
04885 getFarZ[1] = FidVol::getFarZData(1);
04886 getFarZ[2] = FidVol::getFarZData(2);
04887 getFarZ[3] = FidVol::getFarZData(3);
04888 } else {
04889 getFarZ[0] = FidVol::getFarZMC(0);
04890 getFarZ[1] = FidVol::getFarZMC(1);
04891 getFarZ[2] = FidVol::getFarZMC(2);
04892 getFarZ[3] = FidVol::getFarZMC(3);
04893 }
04894
04895 // Now - examine each of the regions and determine if the vertex is in there
04896
04897 const int GAPPLANE = 249; // plane at the SM gap, without scintillator
04898
04899 nearEdge = radius2 >= (FidVol::getFarRouter()*FidVol::getFarRouter());
04900 nearCoil = radius2 < FidVol::getFarRinner()*FidVol::getFarRinner();
04901 nearFront = z < getFarZ[0];
04902 nearGapSouth = z > getFarZ[1] && vtxPlane < GAPPLANE;
04903 nearGapNorth = vtxPlane > GAPPLANE && z < getFarZ[2];
04904 nearBack = vtxZ - FidVol::getEvtVtxZOffset() > getFarZ[3];
04905
04906 return kEdge * nearEdge +
04907 kCoilHole * nearCoil +
04908 kFront * nearFront +
04909 kGapSouth * nearGapSouth +
04910 kGapNorth * nearGapNorth +
04911 kBack * nearBack;
04912 }
|
|
||||||||||||
|
Retrieves the (rock muon) detector region from SNTP.
Definition at line 4780 of file NuReco.cxx. References det, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), VldContext::GetSimFlag(), RecHeader::GetVldContext(), NtpSRVertex::plane, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04781 {
04782 // Get the detector and simFlag
04783 Detector::Detector_t det = rec.GetHeader().GetVldContext().GetDetector();
04784 SimFlag::SimFlag_t simFlag = rec.GetHeader().GetVldContext().GetSimFlag();
04785
04786 return GetRegion(det, simFlag, vtx.x, vtx.y, vtx.z, vtx.plane);
04787 }
|
|
||||||||||||||||
|
Definition at line 3012 of file NuReco.cxx. References GetShowerWithIndexX(), GetTrackWithIndexX(), NtpSRTrack::index, NtpSRShower::index, NuEvent::mcShw1, NuEvent::mcShw2, NuEvent::mcShw3, NuEvent::mcShw4, NuEvent::mcShw5, NuEvent::mcTrk1, NuEvent::mcTrk2, NuEvent::mcTrk3, NtpTHTrack::neumc, NtpTHShower::neumc, NuEvent::shwExists1, NuEvent::shwExists2, NuEvent::shwExists3, NuEvent::shwExists4, NuEvent::shwExists5, NtpStRecord::thshw, NtpStRecord::thtrk, NuEvent::trkExists1, NuEvent::trkExists2, and NuEvent::trkExists3. Referenced by GetTruthInfo(). 03015 {
03016 Int_t mcShw1 = -1;
03017 Int_t mcShw2 = -1;
03018 Int_t mcShw3 = -1;
03019 Int_t mcShw4 = -1;
03020 Int_t mcShw5 = -1;
03021
03022 //get the shower mc
03023 TClonesArray& thshwTca=(*ntp.thshw);//TruthHelper Shower TCA
03024 const Int_t numthshws=thshwTca.GetEntriesFast();
03025 if (numthshws>0) {
03026 if (nu.shwExists1){
03027 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,0);
03028 if (pshw) {
03029 const NtpSRShower& shw=*pshw;
03030 const NtpTHShower& thshw=
03031 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03032 mcShw1=thshw.neumc;
03033 }
03034 }
03035 if (nu.shwExists2){
03036 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,1);
03037 if (pshw) {
03038 const NtpSRShower& shw=*pshw;
03039 const NtpTHShower& thshw=
03040 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03041 mcShw2=thshw.neumc;
03042 }
03043 }
03044 if (nu.shwExists3){
03045 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,2);
03046 if (pshw) {
03047 const NtpSRShower& shw=*pshw;
03048 const NtpTHShower& thshw=
03049 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03050 mcShw3=thshw.neumc;
03051 }
03052 }
03053 if (nu.shwExists4){
03054 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,3);
03055 if (pshw) {
03056 const NtpSRShower& shw=*pshw;
03057 const NtpTHShower& thshw=
03058 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03059 mcShw4=thshw.neumc;
03060 }
03061 }
03062 if (nu.shwExists5){
03063 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,4);
03064 if (pshw) {
03065 const NtpSRShower& shw=*pshw;
03066 const NtpTHShower& thshw=
03067 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03068 mcShw5=thshw.neumc;
03069 }
03070 }
03071 }
03072
03073 nu.mcShw1 = mcShw1;
03074 nu.mcShw2 = mcShw2;
03075 nu.mcShw3 = mcShw3;
03076 nu.mcShw4 = mcShw4;
03077 nu.mcShw5 = mcShw5;
03078
03079
03080 Int_t mcTrk1 = -1;
03081 Int_t mcTrk2 = -1;
03082 Int_t mcTrk3 = -1;
03083
03084 //get the track mc
03085 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Track TCA
03086 const Int_t numthtrks=thtrkTca.GetEntriesFast();
03087 if (numthtrks>0) {
03088 if (nu.trkExists1){
03089 const NtpSRTrack* ptrk =
03090 this->GetTrackWithIndexX(ntp,evt,0);
03091 if (ptrk){
03092 const NtpSRTrack& trk=*ptrk;
03093 const NtpTHTrack& thtrk=
03094 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03095 mcTrk1=thtrk.neumc;
03096 }
03097 }
03098 if (nu.trkExists2){
03099 const NtpSRTrack* ptrk =
03100 this->GetTrackWithIndexX(ntp,evt,1);
03101 if (ptrk){
03102 const NtpSRTrack& trk=*ptrk;
03103 const NtpTHTrack& thtrk=
03104 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03105 mcTrk2=thtrk.neumc;
03106 }
03107 }
03108 if (nu.trkExists3){
03109 const NtpSRTrack* ptrk =
03110 this->GetTrackWithIndexX(ntp,evt,2);
03111 if (ptrk){
03112 const NtpSRTrack& trk=*ptrk;
03113 const NtpTHTrack& thtrk=
03114 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03115 mcTrk3=thtrk.neumc;
03116 }
03117 }
03118 }
03119
03120 nu.mcTrk1 = mcTrk1;
03121 nu.mcTrk2 = mcTrk2;
03122 nu.mcTrk3 = mcTrk3;
03123 }
|
|
|
Definition at line 642 of file NuReco.cxx. References GetShowerEnergyCor(), MAXMSG, NuEvent::shwEnLinCCCor, NuEvent::shwEnLinCCNoCor, and NuEvent::shwEnNoCor. Referenced by CalcShwVariables(), and NuDSTAna::FDTestAna(). 00643 {
00644 //this function adds a level of indirection
00645 //to allow flexibility
00646 //could sum up all the shower energies here for example
00647 //or just return the primary shower
00648
00649 //if (nu.redoShwEnCor) {
00650 if (true) {
00651 // Recover old behavior as needed
00652 double shwEn = nu.shwEnLinCCNoCor;
00653 if (shwEn == -1) {
00654 MAXMSG("NuReco",Msg::kInfo,1)
00655 <<"GetShowerEnergy: This appears to be an old DST -- using shwEnNoCor instead of shwEnLinCCNoCor"<<endl;
00656 shwEn = nu.shwEnNoCor;
00657 }
00658 MAXMSG("NuReco",Msg::kInfo,1)
00659 <<"GetShowerEnergy: (re)doing EnergyCorrection"<<endl;
00660 Double_t shwEnCor=this->GetShowerEnergyCor(shwEn, CandShowerHandle::kCC,nu);
00661 MAXMSG("NuReco",Msg::kDebug,300)
00662 <<"Redo energy cor: shwEnNoCor="<<nu.shwEnNoCor
00663 <<", shwEnCor="<<shwEnCor<<endl;
00664 if (shwEnCor == -1) {
00665 MAXMSG("NuReco",Msg::kInfo,1)
00666 <<"GetShowerEnergy: Shower Energy Correction faulty -- using uncorrected"<<endl;
00667 return nu.shwEnNoCor;
00668 }
00669 return shwEnCor;
00670 }
00671 else {
00672 MAXMSG("NuReco",Msg::kInfo,1)
00673 <<"GetShowerEnergy: NOT (re)doing EnergyCorrection"<<endl;
00674 return nu.shwEnLinCCCor;
00675 }
00676 }
|
|
||||||||||||||||
|
Definition at line 704 of file NuReco.cxx. References EnergyCorrections::FullyCorrectShowerEnergy(), GetVldContext(), and NuEvent::releaseType. Referenced by NuExtraction::ExtractShwInfo(), NuDSTAna::FDTestAna(), NuPlots::FillShwHistos(), GetShowerEnergyCC(), GetShowerEnergyNC(), and NuDSTAna::NDTestAna(). 00707 {
00708 //check if shower exists
00709 if (shwEn<0) return -1;
00710
00711 //calc energy correction
00712 VldContext vc=this->GetVldContext(nu);
00713 return EnergyCorrections::FullyCorrectShowerEnergy(shwEn,st,vc,nu.releaseType);
00714 }
|
|
||||||||||||
|
res will be filled with an energy resolution estimate
Definition at line 793 of file NuReco.cxx. References C, NuKDTree< T, D >::Committed(), det, NuEvent::detector, Detector::Detector_t, NuEvent::energyMC, NuKDTree< T, D >::FindNearestPts(), InitializeShowerEnergykNN(), NuEvent::nplaneShw, NuEvent::nshw, NuEvent::shwEn, NuEvent::shwEnCor2, and NuEvent::trkShwEnNearDW. Referenced by CalcShwVariables(). 00794 {
00795 // If something went wrong with the kNN initialization, then bail out
00796 // early on subsequent attempts
00797 static bool treeOK = true;
00798 if(!treeOK){
00799 if(res) *res = -1;
00800 return -1;
00801 }
00802
00803 // Initialize the tree the first time we are called
00804 static NuKDTree<double, 3> kdTree;
00805 // Energy correction constants
00806 static vector<double> C;
00807 static double cutoff_lo;
00808 static double cutoff_hi;
00809
00810 static int num_neighbours;
00811
00812 if(!kdTree.Committed()){
00813 const Detector::Detector_t det = Detector::Detector_t(nu.detector);
00814 const bool success = InitializeShowerEnergykNN(&kdTree, det,
00815 C, cutoff_lo, cutoff_hi,
00816 num_neighbours);
00817
00818 if(!success){
00819 treeOK = false;
00820 if(res) *res = -1;
00821 return -1;
00822 }
00823 } // end if not committed
00824
00825 double sum2Shw = nu.shwEn;
00826 if(nu.nshw > 1) sum2Shw += nu.shwEnCor2;
00827
00828 const NuKDPt<3> pt(nu.trkShwEnNearDW, double(nu.nplaneShw), sum2Shw);
00829
00830 // Sometimes need extra neighbours for MC because results will be this
00831 // same event. This can happen more than once because the same true event
00832 // gets split by the reconstruction. Also can happen because rock overlays
00833 // appear to be reused - although these should fail fiducial and PID cuts.
00834 vector<NuKDPtWithPayload<double, 3> > nearest =
00835 kdTree.FindNearestPts(pt, num_neighbours + 10);
00836
00837 // Calculate the average true energy of all the nearby events
00838 double tot = 0;
00839 double sqrTot = 0;
00840 unsigned int N = num_neighbours;
00841 for(unsigned int n = 0; n < N; ++n){
00842 const double pl = nearest[n].payload;
00843
00844 // We found ourself in the list of neighbours, should ignore this point
00845 // as it can't happen for data
00846 if(nu.energyMC > 0 && fabs(pl-nu.energyMC) < 1e-5){
00847 ++N; // Do one more neighbour to make up for missing this one
00848 // If this assertion fails it means we found more than 10 events with
00849 // the same energyMC as this one. That's probably worth investigating.
00850 assert(N <= nearest.size());
00851 continue;
00852 }
00853
00854 tot += pl;
00855 sqrTot += pl*pl;
00856 }
00857
00858 const double sd = sqrt(N*sqrTot-tot*tot)/num_neighbours;
00859 if(res) *res = sd;
00860
00861 double ret = tot/num_neighbours;
00862
00863 // Apply energy correction - polynomial in log(E)
00864 // NB - the resolution estimate isn't scaled
00865 double clampedE = ret;
00866 if(clampedE < cutoff_lo) clampedE = cutoff_lo;
00867 if(clampedE > cutoff_hi) clampedE = cutoff_hi;
00868 const double y = TMath::Log(clampedE);
00869 const double y2 = y2*y; const double y3 = y3*y;
00870 const double y4 = y4*y; const double y5 = y5*y;
00871 const double y6 = y6*y; const double y7 = y7*y;
00872 const double corrFactor = C[0]+C[1]*y+C[2]*y2+C[3]*y3+
00873 C[4]*y4+C[5]*y5+C[6]*y6+C[7]*y7;
00874
00875 ret *= corrFactor;
00876
00877 return ret;
00878 }
|
|
|
Definition at line 679 of file NuReco.cxx. References GetShowerEnergyCor(), MAXMSG, NuEvent::shwEnLinNCCor, NuEvent::shwEnLinNCNoCor, and NuEvent::shwEnNoCor. Referenced by CalcShwVariables(). 00680 {
00681 //this function adds a level of indirection
00682 //to allow flexibility
00683 //could sum up all the shower energies here for example
00684 //or just return the primary shower
00685
00686 //if (nu.redoShwEnCor) {
00687 if (true) {
00688 MAXMSG("NuReco",Msg::kInfo,1)
00689 <<"GetShowerEnergy: (re)doing EnergyCorrection"<<endl;
00690 Double_t shwEnCor=this->GetShowerEnergyCor(nu.shwEnLinNCNoCor, CandShowerHandle::kNC,nu);
00691 MAXMSG("NuReco",Msg::kDebug,300)
00692 <<"Redo energy cor: shwEnNoCor="<<nu.shwEnNoCor
00693 <<", shwEnCor="<<shwEnCor<<endl;
00694 return shwEnCor;
00695 }
00696 else {
00697 MAXMSG("NuReco",Msg::kInfo,1)
00698 <<"GetShowerEnergy: NOT (re)doing EnergyCorrection"<<endl;
00699 return nu.shwEnLinNCCor;
00700 }
00701 }
|
|
||||||||||||||||
|
Calculates the shower energy 'near' the track vertex. nearShowerEnergyDW is the old version of this value, which used the (broken) deweighted shower energy. Returns zero if there are no showers near enough Definition at line 4617 of file NuReco.cxx. References det, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), GetShwMedianTime(), RecHeader::GetVldContext(), NtpSRShowerPulseHeight::linCCgev, NtpStRecord::shw, NtpSRShower::shwph, NtpSRVertex::t, NtpSRShower::vtx, NtpSRTrack::vtx, NtpSRShowerPulseHeight::wtCCgev, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04620 {
04621 Double_t nearshowere = 0;
04622 // Sum of deweighted shower energies. This is probably not a good idea - but
04623 // we used to define this this way before so someone might want the old value
04624 if(nearshowereDW) *nearshowereDW = 0;
04625
04626 bool hasgoodshower = false;
04627
04628 const double trkTime = trk.vtx.t;
04629
04630 const Detector::Detector_t det = ntp.GetHeader().GetVldContext().GetDetector();
04631
04632 // Loop over all shower entries
04633 for(Int_t it = 0; it < ntp.shw->GetEntries(); it++){
04634 // Grab the shower
04635 const NtpSRShower *shw = dynamic_cast<NtpSRShower *>(ntp.shw->At(it));
04636
04637 // Apply +75ns, -25ns cut in near detector only - see docdb 6667-v1 p4
04638 if(det == Detector::kNear){
04639 const double shwTime = GetShwMedianTime(ntp, *shw);
04640 const double dt = (shwTime-trkTime)*1e9; // nanoseconds
04641 if(dt > +75 || dt < -25) continue;
04642 }
04643
04644 // Calculate the vertex seperation
04645 Double_t vertex_sep_sqr = (trk.vtx.x-shw->vtx.x)*(trk.vtx.x-shw->vtx.x) +
04646 (trk.vtx.y-shw->vtx.y)*(trk.vtx.y-shw->vtx.y) +
04647 (trk.vtx.z-shw->vtx.z)*(trk.vtx.z-shw->vtx.z);
04648
04649 // If the shower vertex is 'near' the track vertex (<1m away)
04650 if(vertex_sep_sqr < 1.0){
04651 nearshowere += shw->shwph.linCCgev;
04652 if(nearshowereDW) *nearshowereDW += shw->shwph.wtCCgev;
04653 hasgoodshower = true;
04654 }
04655 }
04656
04657 // Match old definition
04658 if(!hasgoodshower && nearshowereDW) *nearshowereDW = -1;
04659
04660 return nearshowere;
04661 }
|
|
||||||||||||||||
|
Definition at line 1882 of file NuReco.cxx. References NtpSREvent::nshower, NtpSREvent::shw, and NtpStRecord::shw. Referenced by NuCounter::CountTrkStdhepId(), NuAnalysis::DemoInfidSRInterface(), NuExtraction::ExtractShwInfo(), NuPlots::FillShwHistos(), GetBestTruthIndex(), GetPrimaryTruthIndex(), and GetSecondaryTruthIndex(). 01885 {
01886 TClonesArray& shwTca=(*ntp.shw);
01887
01888 if (evt.nshower>=index+1 && index>=0){
01889 const NtpSRShower* shw=
01890 dynamic_cast<NtpSRShower*>(shwTca[evt.shw[index]]);
01891 return shw;
01892 }
01893 else {
01894 return 0;
01895 }
01896 }
|
|
||||||||||||
|
Which plane does shw have the most ph deposited in?
Definition at line 2816 of file NuReco.cxx. References MAXMSG, NtpSRShower::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRShower::stp, and NtpStRecord::stp. Referenced by NuExtraction::ExtractShwInfo(). 02818 {
02819 typedef map<UShort_t, Float_t> planePH_t;
02820 planePH_t planePH;
02821
02822 for(int i = 0; i < shw.nstrip; ++i){
02823 const NtpSRStrip* stp = (NtpSRStrip*)ntp.stp->At(shw.stp[i]);
02824
02825 planePH[stp->plane] += stp->ph0.pe+stp->ph1.pe;
02826 }
02827
02828 Float_t maxPH = 0;
02829 Int_t maxPlane = -1;
02830
02831 for(planePH_t::iterator it = planePH.begin(); it != planePH.end(); ++it){
02832 if(it->second > maxPH){
02833 maxPH = it->second;
02834 maxPlane = it->first;
02835 }
02836 }
02837
02838 MAXMSG("NuReco", Msg::kDebug, 100) << "Max plane=" << maxPlane << endl;
02839
02840 return maxPlane;
02841 }
|
|
||||||||||||
|
Definition at line 2774 of file NuReco.cxx. References NtpSRStrip::index, MAXMSG, NtpSRShower::nstrip, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSRShower::stp, NtpStRecord::stp, NtpSRStrip::strip, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSRStrip::tpos. Referenced by GetShowerEnergyNearTrack(). 02776 {
02777 //Msg::LogLevel_t logLevel=Msg::kDebug;
02778
02779 multiset<Double_t> times;
02780 const TClonesArray& stpTca=(*ntp.stp);
02781 //const Int_t numStps=stpTca.GetEntriesFast();
02782
02783 MAXMSG("NuReco",Msg::kDebug,200)
02784 <<"shw.nstrip="<<shw.nstrip<<endl;
02785 for (Int_t i=0;i<shw.nstrip;i++) {
02786 const NtpSRStrip& stp=
02787 *(dynamic_cast<NtpSRStrip*>(stpTca[shw.stp[i]]));
02788
02789 Double_t time=stp.time0;
02790 if (time<0 || time>1) time=stp.time1;
02791
02792 if (time>0 && time<=1) {
02793 times.insert(time);
02794 }
02795 //else just don't put the time in the map - clearly rubbish
02796
02797 MAXMSG("NuReco",Msg::kDebug,500)
02798 <<" Strip index="<<stp.index<<", stp="<<stp.strip
02799 <<", t="<<time
02800 <<", pl="<<stp.plane
02801 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02802 }
02803
02804 //get the median time from the map
02805 multiset<Double_t>::const_iterator it=times.begin();
02806 advance(it,times.size()/2);
02807 Double_t medianTime=*it;
02808 MAXMSG("NuReco",Msg::kDebug,100)
02809 <<"Median time="<<medianTime<<endl;
02810
02811 return medianTime;
02812 }
|
|
|
Definition at line 4481 of file NuReco.cxx. References PlaneOutline::DistanceToEdge(), PlaneOutline::IsInside(), MAXMSG, NuEvent::xEvtVtx, and NuEvent::yEvtVtx. Referenced by CalcEvtVariables(), and TestGetSmallestDeepDistToEdge(). 04482 {
04483 //Deep containment is defined as being within the overlap of
04484 //all different plane coverages
04485
04486 Float_t distUPI=0;
04487 Float_t distVPI=0;
04488 Float_t distUFI=0;
04489 Float_t distVFI=0;
04490 Float_t xedge=0;
04491 Float_t yedge=0;
04492
04493 PlaneOutline po;
04494
04495 //calc for U-view partial
04496 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04497 PlaneCoverage::kNearPartial,
04498 distUPI,xedge,yedge);
04499 Bool_t isInsideUPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04500 PlaneCoverage::kNearPartial);
04501
04502 //calc for U-view full
04503 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04504 PlaneCoverage::kNearFull,
04505 distUFI,xedge,yedge);
04506 Bool_t isInsideUFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04507 PlaneCoverage::kNearFull);
04508
04509 //calc for V-view partial
04510 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04511 PlaneCoverage::kNearPartial,
04512 distVPI,xedge,yedge);
04513 Bool_t isInsideVPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04514 PlaneCoverage::kNearPartial);
04515
04516 //calc for V-view full
04517 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04518 PlaneCoverage::kNearFull,
04519 distVFI,xedge,yedge);
04520 Bool_t isInsideVFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04521 PlaneCoverage::kNearFull);
04522
04523 //check that hit has "deep" containment and find
04524 //the closest distance to the edge of the combined overlap region
04525 Float_t smallestDist=999;
04526 if (isInsideUPI && isInsideUFI && isInsideVPI && isInsideVFI){
04527 if (smallestDist>distUPI) smallestDist=distUPI;
04528 if (smallestDist>distUFI) smallestDist=distUFI;
04529 if (smallestDist>distVPI) smallestDist=distVPI;
04530 if (smallestDist>distVFI) smallestDist=distVFI;
04531 }
04532
04533 //if not deeply contained then find the closest distance
04534 //to the combined overlap region
04535 //this is actually the furthest distance!
04536 if (smallestDist==999) {
04537 smallestDist=0;
04538 if (smallestDist<distUPI && !isInsideUPI) smallestDist=distUPI;
04539 if (smallestDist<distUFI && !isInsideUFI) smallestDist=distUFI;
04540 if (smallestDist<distVPI && !isInsideVPI) smallestDist=distVPI;
04541 if (smallestDist<distVFI && !isInsideVFI) smallestDist=distVFI;
04542
04543 //make the distance negative for hits outside deep containment
04544 smallestDist*=-1;
04545 }
04546
04547 if (smallestDist==999) cout<<"Ahhhhhhh"<<endl;
04548
04549 MAXMSG("NuReco",Msg::kDebug,500)
04550 <<"smallestDist="<<smallestDist
04551 <<", UPI="<<distUPI<<" (in="<<isInsideUPI<<")"
04552 <<", UFI="<<distUFI<<" (in="<<isInsideUFI<<")"
04553 <<", VPI="<<distVPI<<" (in="<<isInsideVPI<<")"
04554 <<", VFI="<<distVFI<<" (in="<<isInsideVFI<<")"<<endl;
04555
04556 return smallestDist;
04557 }
|
|
||||||||||||
|
Definition at line 1739 of file NuReco.cxx. 01740 {
01741 Int_t charge=qp<0 ? -1 : 1;
01742 if (qp==0) {
01743 if (isReverse) charge=+1;//assume numubar!!!
01744 else charge=-1;//assume numu!!!
01745 }
01746
01747 return charge;
01748 }
|
|
||||||||||||
|
Definition at line 1717 of file NuReco.cxx. References NtpSRTrack::fit, MAXMSG, NtpSRTrack::momentum, NtpSRFitTrack::pass, and NtpSRMomentum::qp. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSAFitInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 01718 {
01719 Int_t charge=this->GetTrackCharge(trk.momentum.qp, isReverse);
01720 if (trk.momentum.qp==0) {
01721 if (isReverse) {
01722 MAXMSG("NuReco",Msg::kInfo,3)
01723 <<"Note: assigning charge to be +1 since trk.momentum.qp="
01724 <<trk.momentum.qp<<", trkfitpass="<<trk.fit.pass
01725 <<" and this is RHC" << endl;
01726 }
01727 else {
01728 MAXMSG("NuReco",Msg::kInfo,3)
01729 <<"Note: assigning charge to be -1 since trk.momentum.qp="
01730 <<trk.momentum.qp<<", trkfitpass="<<trk.fit.pass
01731 <<" and this is FHC" << endl;
01732 }
01733 }
01734 return charge;
01735 }
|
|
|
Definition at line 532 of file NuReco.cxx. References GetTrackEnergyFromCurv(), GetTrackEnergyFromRange(), UseCurvature(), and UseRange(). Referenced by CalcTrkVariables(). 00533 {
00534 //get the track energy
00535 //FC
00536 if (this->UseRange(nu)) {
00537 return this->GetTrackEnergyFromRange(nu);
00538 }
00539 //PC
00540 else if (this->UseCurvature(nu)){
00541 return this->GetTrackEnergyFromCurv(nu);
00542 }
00543 else {
00544 cout<<"Ahhhh, can't use range or curvature"<<endl;
00545 return -1;
00546 }
00547 }
|
|
|
Definition at line 600 of file NuReco.cxx. References GetTrackEnergyFromCurvCor(), MAXMSG, NuEvent::qp, NuEvent::trkEnCorCurv, and NuEvent::trkEnCurv. Referenced by CalcTrkVariables(), and GetTrackEnergy(). 00601 {
00602 //if (nu.redoTrkEnCor) {
00603 if (true) {
00604 MAXMSG("NuReco",Msg::kInfo,1)
00605 <<"GetTrackEnergyFromCurv: (re)doing EnergyCorrection"<<endl;
00606 Double_t trkEn=this->GetTrackEnergyFromCurvCor(nu.qp,nu);
00607 Double_t p=-1;
00608 if (nu.qp) p=1./nu.qp;
00609 MAXMSG("NuReco",Msg::kDebug,3000)
00610 <<"Curv: trkEn="<<trkEn<<", trkEnCorCurv="<<nu.trkEnCorCurv
00611 <<", 1/qp="<<p<<endl;
00612 return trkEn;
00613 }
00614 else {
00615 MAXMSG("NuReco",Msg::kInfo,1)
00616 <<"GetTrackEnergyFromCurv: NOT (re)doing EnergyCorrection"<<endl;
00617 return nu.trkEnCurv;
00618 }
00619 }
|
|
||||||||||||
|
this function should only ever be called if the track exists there is no protection against value being initialised to -1 since it is a perfectly valid value for qp Definition at line 623 of file NuReco.cxx. References CapTrackMomentum(), GetVldContext(), and NuEvent::releaseType. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), NuExtraction::ExtractThirdTrkInfo(), and GetTrackEnergyFromCurv(). 00625 {
00629 if (qp!=0){
00630 Double_t p=1./qp;
00631 this->CapTrackMomentum(p);
00632 VldContext vc=this->GetVldContext(nu);
00633 p=EnergyCorrections::FullyCorrectSignedMomentumFromCurvature
00634 (p,vc,nu.releaseType);
00635 return sqrt(pow(p,2)+pow(0.105658357*(Munits::GeV),2));
00636 }
00637 else return 0.;
00638 }
|
|
|
Definition at line 551 of file NuReco.cxx. References GetTrackEnergyFromRangeCor(), MAXMSG, NuEvent::trkEnCorRange, NuEvent::trkEnRange, and NuEvent::trkMomentumRange. Referenced by CalcTrkVariables(), and GetTrackEnergy(). 00552 {
00553 //if (nu.redoTrkEnCor) {
00554 if (true) {
00555 MAXMSG("NuReco",Msg::kInfo,1)
00556 <<"GetTrackEnergyFromRange: (re)doing EnergyCorrection"<<endl;
00557 Double_t trkEn=this->GetTrackEnergyFromRangeCor
00558 (nu.trkMomentumRange,nu);
00559 MAXMSG("NuReco",Msg::kDebug,3000)
00560 <<"Range: trkEn="<<trkEn<<", trkEnCorRange="<<nu.trkEnCorRange
00561 <<", trkMomentumRange="<<nu.trkMomentumRange<<endl;
00562 return trkEn;
00563 }
00564 else {
00565 MAXMSG("NuReco",Msg::kInfo,1)
00566 <<"GetTrackEnergyFromRange: NOT (re)doing EnergyCorrection"<<endl;
00567 return nu.trkEnRange;
00568 }
00569 }
|
|
||||||||||||
|
Definition at line 582 of file NuReco.cxx. References CapTrackMomentum(), GetVldContext(), and NuEvent::releaseType. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), NuExtraction::ExtractThirdTrkInfo(), and GetTrackEnergyFromRange(). 00584 {
00585 //check if the track exists
00586 if (trkMomentumRange<0) return -1;
00587
00588 this->CapTrackMomentum(trkMomentumRange);
00589 //calc the corrected momentum from range
00590 Double_t mr=0;
00591 VldContext vc=this->GetVldContext(nu);
00592 mr=EnergyCorrections::FullyCorrectMomentumFromRange
00593 (trkMomentumRange,vc,nu.releaseType);
00594
00595 return sqrt(pow(mr,2)+pow(0.105658357*(Munits::GeV),2));
00596 }
|
|
||||||||||||
|
Definition at line 4684 of file NuReco.cxx. References NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04685 {
04686 NtpSRStrip * firststrip =
04687 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[0]));
04688 return firststrip->strip;
04689 }
|
|
||||||||||||
|
Definition at line 4721 of file NuReco.cxx. References NtpSRStrip::planeview, NtpSRTrack::stp, and NtpStRecord::stp. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04722 {
04723 NtpSRStrip * firststrip =
04724 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[0]));
04725 if(firststrip->planeview == PlaneView::kU) return true;
04726 else return false;
04727 }
|
|
||||||||||||
|
Definition at line 4693 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRStrip::planeview, NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04694 {
04695 // Look for the first U strip in the track
04696 for(int i = 0; i < trk.nstrip; i++) {
04697 NtpSRStrip * firststrip =
04698 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04699 if(firststrip->planeview == PlaneView::kU) return firststrip->strip;
04700 }
04701 // There is no U strip (unlikely)
04702 return -1;
04703 }
|
|
||||||||||||
|
Definition at line 4707 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRStrip::planeview, NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04708 {
04709 // Look for the first U strip in the track
04710 for(int i = 0; i < trk.nstrip; i++) {
04711 NtpSRStrip * firststrip =
04712 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04713 if(firststrip->planeview == PlaneView::kV) return firststrip->strip;
04714 }
04715 // There is no V strip (unlikely)
04716 return -1;
04717 }
|
|
||||||||||||
|
Definition at line 4731 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04732 {
04733 NtpSRStrip * laststrip =
04734 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[trk.nstrip-1]));
04735 return laststrip->strip;
04736 }
|
|
||||||||||||
|
Definition at line 4768 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRStrip::planeview, NtpSRTrack::stp, and NtpStRecord::stp. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04769 {
04770 NtpSRStrip * laststrip =
04771 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[trk.nstrip-1]));
04772 if(laststrip->planeview == PlaneView::kU) return true;
04773 else return false;
04774 }
|
|
||||||||||||
|
Definition at line 4740 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRStrip::planeview, NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04741 {
04742 // Loop backwards for the first u strip
04743 for(Int_t i = trk.nstrip-1; i >= 0; i--) {
04744 NtpSRStrip * laststrip =
04745 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04746 if(laststrip->planeview == PlaneView::kU) return laststrip->strip;
04747 }
04748 // There is no U strip (unlikely)
04749 return -1;
04750 }
|
|
||||||||||||
|
Definition at line 4754 of file NuReco.cxx. References NtpSRTrack::nstrip, NtpSRStrip::planeview, NtpSRTrack::stp, NtpStRecord::stp, and NtpSRStrip::strip. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 04755 {
04756 // Loop backwards for the first v strip
04757 for(Int_t i = trk.nstrip-1; i >= 0; i--) {
04758 NtpSRStrip * laststrip =
04759 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04760 if(laststrip->planeview == PlaneView::kV) return laststrip->strip;
04761 }
04762 // There is no U strip (unlikely)
04763 return -1;
04764 }
|
|
|
Definition at line 5023 of file NuReco.cxx. References NuEvent::dirCosNu. Referenced by CalcTrkVariables(). 05024 {
05025 Double_t bestmomentum = -1;
05026
05027 return bestmomentum*TMath::Sin(TMath::ACos(nu.dirCosNu));
05028 }
|
|
||||||||||||||||
|
|
Definition at line 2065 of file NuReco.cxx. References MAXMSG, NtpSRTrack::nstrip, and NtpSRTrack::stpfitqp. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 02066 {
02067 Int_t nPositive=0;
02068 Int_t nNegative=0;
02069 Int_t nZero=0;
02070
02071 for (Int_t i=0;i<trk.nstrip;i++){
02072 //Double_t x=trk.stpx[i];
02073 Double_t stpfitqp=trk.stpfitqp[i];
02074
02075 //count number of strips with positive and negative qp
02076 if (stpfitqp>0) nPositive++;
02077 else if (stpfitqp<0) nNegative++;
02078 else nZero++;
02079
02080 MAXMSG("NuReco",Msg::kDebug,1000)
02081 <<"stpfitqp="<<stpfitqp<<endl;
02082 }
02083
02084 Double_t qpFraction=-1;
02085 if (trk.nstrip>0) qpFraction=1.*nPositive/trk.nstrip;
02086
02087 if (trk.nstrip!=nPositive+nNegative) {
02088 MAXMSG("NuReco",Msg::kDebug,100)
02089 <<"??? stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
02090 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
02091 <<", trkn="<<trk.nstrip
02092 <<", qpFraction="<<qpFraction
02093 <<endl;
02094 }
02095 else {
02096 MAXMSG("NuReco",Msg::kDebug,100)
02097 <<"stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
02098 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
02099 <<", trkn="<<trk.nstrip
02100 <<", qpFraction="<<qpFraction
02101 <<endl;
02102 }
02103
02104 return qpFraction;
02105 }
|
|
||||||||||||
|
Get the true track identity.
Definition at line 5308 of file NuReco.cxx. References NtpMCStdHep::IdHEP, NtpSRTrack::index, MSG, NtpStRecord::stdhep, NtpStRecord::thtrk, and NtpTHTrack::trkstdhep. Referenced by NuExtraction::ExtractPrimaryTrkInfo(), NuExtraction::ExtractSecondTrkInfo(), and NuExtraction::ExtractThirdTrkInfo(). 05309 {
05310 // Firstly, find the THTrack object by searching the MC info
05311 //get the track mc
05312 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Track TCA
05313 // Are there any entries?
05314 if (thtrkTca.GetEntriesFast() <= 0) return 0;
05315
05316 // Get the THTrack object
05317 const NtpTHTrack* thtrk= dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
05318
05319 if(!thtrk) {
05320 MSG("NuReco",Msg::kError) << "Truth track information available but track "
05321 << trk.index << " not found!" << endl;
05322 return 0;
05323 }
05324
05325 // Grab the Stdhep entry
05326 NtpMCStdHep * stdhep =
05327 dynamic_cast<NtpMCStdHep *>(ntp.stdhep->At(thtrk->trkstdhep));
05328
05329 if(!stdhep) {
05330 MSG("NuReco",Msg::kError) << "Truth Helper, but no stdhep!\n";
05331 return 0;
05332 }
05333
05334 return stdhep->IdHEP;
05335 }
|
|
||||||||||||||||
|
||||||||||||||||
|
Definition at line 3146 of file NuReco.cxx. References GetPrimaryTruthIndex(), GetSecondaryTruthIndex(), NtpMCTruth::iaction, MAXMSG, NtpStRecord::mc, NuEvent::mc, NtpMCTruth::p4mu1, NtpMCTruth::p4neu, NtpMCTruth::p4shw, NtpMCTruth::p4tgt, NuEvent::simFlag, and NtpMCTruth::y. Referenced by NuDemoModule::Ana(), NuAnalysis::ChargeSignCut(), NuAnalysis::DoExtractions(), NuAnalysis::Efficiencies(), NuAnalysis::EnergySpect(), NuAnalysis::LIRejectionTest(), and NuAnalysis::LoopOverTruthInfo(). 03148 {
03149 //get the index of the mc neutrino interaction
03150 //Mad uses the primary track, so use that here
03151 this->GetPrimaryTruthIndex(ntp,evt,nu);
03152 this->GetSecondaryTruthIndex(ntp,evt,nu);
03153
03154 //MadBase::LoadTruthForRecoTH does not use the "best" trk/shw
03155 //this->GetBestTruthIndex(ntp,evt,nu);
03156
03157 //use the index
03158 if (nu.mc<0) {
03159 MAXMSG("NuReco",Msg::kInfo,1)
03160 <<"Not Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
03161 return;
03162 }
03163 MAXMSG("NuReco",Msg::kInfo,1)
03164 <<"Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
03165
03166 //now get the mc object (neutrino interaction) that
03167 //corresponds to the trk using the nu.mc index
03168 TClonesArray& mcTca=(*ntp.mc);//TruthHelper Track TCA
03169 const NtpMCTruth& mc=
03170 *(dynamic_cast<NtpMCTruth*>(mcTca[nu.mc]));
03171
03172 //extract the info from this mc object
03173 this->GetTruthInfo(ntp,mc,nu);
03174
03175 if (mc.iaction==1){
03176 MAXMSG("NuReco",Msg::kDebug,300)
03177 <<"y="<<mc.y<<endl
03178 <<"p4neu[3]="<<mc.p4neu[3]<<", p4tgt[3]="<<mc.p4tgt[3]
03179 <<", sum="<<mc.p4neu[3]+mc.p4tgt[3]
03180 <<endl
03181 <<"p4mu1[3]="<<mc.p4mu1[3]<<", mc.p4shw[3]="<<mc.p4shw[3]
03182 <<", sum="<<fabs(mc.p4mu1[3])+mc.p4shw[3]
03183 <<", labY="<<mc.p4shw[3]/(fabs(mc.p4mu1[3])+mc.p4shw[3])
03184 <<endl;
03185 }
03186 else if (mc.iaction==0){
03187 MAXMSG("NuReco",Msg::kDebug,300)
03188 <<"NC: y="<<mc.y<<endl
03189 <<"p4neu[3]="<<mc.p4neu[3]
03190 <<", p4tgt[3]="<<mc.p4tgt[3]
03191 <<endl
03192 <<"p4mu1[3]="<<mc.p4mu1[3]<<", p4shw[3]="<<mc.p4shw[3]
03193 <<endl;
03194 }
03195 }
|
|
|
Definition at line 1954 of file NuReco.cxx. References NuEvent::detector, NuEvent::simFlag, NuEvent::timeNanoSec, and NuEvent::timeSec. Referenced by CalcExtraTruthVariables(), NuExtraction::ExtractBeamInfoDB(), NuExtraction::ExtractCoilCurrent(), NuExtraction::ExtractCoilInfo(), NuExtraction::ExtractDataQuality(), NuExtraction::ExtractSAFitInfo(), NuExtraction::ExtractTimeToNearestSpill(), NuPlots::FillDetectorEdge(), GetShowerEnergyCor(), GetTrackEnergyFromCurvCor(), and GetTrackEnergyFromRangeCor(). 01955 {
01956 //build the vc from the NuEvent
01957 VldTimeStamp time(nu.timeSec,nu.timeNanoSec);
01958 //VldContext (const Detector::Detector_t &detector, const SimFlag::SimFlag_t mcflag, const VldTimeStamp &time)
01959
01960 VldContext vc(static_cast<Detector::Detector_t>(nu.detector),
01961 static_cast<SimFlag::SimFlag_t>(nu.simFlag),time);
01962 return vc;
01963 }
|
|
||||||||||||||||||||||||||||
|
Helper function for GetShowerEnergykNN. C, cutoff_lo and cutoff_hi are the constants to use in kNN energy corrections. Will be called automatically on first call to GetShowerEnergykNN Returns false if we fail to find the training file Definition at line 717 of file NuReco.cxx. References NuKDTree< T, D >::AddPt(), C, NuKDTree< T, D >::Commit(), det, NuLibrary::general, Registry::Get(), Registry::GetDouble(), NuGeneral::GetReleaseDirToUse(), NuLibrary::Instance(), and MSG. Referenced by GetShowerEnergykNN(). 00723 {
00724 MSG("NuReco", Msg::kInfo) << "Initializing shower energy kNN..." << endl;
00725
00726 // Protect people's gDirectory from being unexpectedly stomped on
00727 TDirectory* restore = gDirectory;
00728
00729 const TString fname = (det == Detector::kFar) ? "kNNShwEnTrainFar.root"
00730 : "kNNShwEnTrainNear.root";
00731
00732 const NuGeneral& gen = NuLibrary::Instance().general;
00733 const TString releaseDir = gen.GetReleaseDirToUse("NtupleUtils");
00734
00735 const TString fullPath = releaseDir+"/NtupleUtils/data/"+fname;
00736
00737 TFile trainFile(fullPath, "READ");
00738 // Restore the global directory the instant we have changed it
00739 restore->cd();
00740
00741 if(trainFile.IsZombie()){
00742 MSG("NuReco", Msg::kError) << "No kNN shower training file at "
00743 << fullPath
00744 << " variables will not be filled" << endl;
00745 return false;
00746 }
00747
00748 TTree* tr = (TTree*)trainFile.Get("train");
00749
00750 double trkShwEnNearDW, nplaneShw, sum2Shw, shwEnMC;
00751
00752 const int status = tr->SetBranchAddress("trkShwEnNearDW", &trkShwEnNearDW);
00753 if(status < 0) tr->SetBranchAddress("trkShwEnNear", &trkShwEnNearDW);
00754
00755 tr->SetBranchAddress("nplaneShw", &nplaneShw);
00756 tr->SetBranchAddress("sum2Shw", &sum2Shw);
00757 tr->SetBranchAddress("shwEnMC", &shwEnMC);
00758
00759 TRandom3 r;
00760 const int N = tr->GetEntries();
00761 for(int n = 0; n < N; ++n){
00762 tr->GetEntry(n);
00763 const NuKDPt<3> pt(trkShwEnNearDW,
00764 nplaneShw + r.Gaus(0, 1e-5), // avoid all the same
00765 sum2Shw);
00766 kdTree->AddPt(shwEnMC, pt);
00767 }
00768
00769 kdTree->Commit();
00770
00771 // Extract the energy corrections constants from the training file
00772 Registry* rEnCorr = (Registry*)trainFile.Get("energy_corrections");
00773 // Polynomial goes up to 7th order
00774 C.resize(8);
00775 for(unsigned int n = 0; n < C.size(); ++n)
00776 C[n] = rEnCorr->GetDouble(TString::Format("C%d", n));
00777 cutoff_lo = rEnCorr->GetDouble("cutoff_lo");
00778 cutoff_hi = rEnCorr->GetDouble("cutoff_hi");
00779
00780 const bool success = rEnCorr->Get("num_neighbours", neighbours);
00781 if(!success){
00782 MSG("NuReco", Msg::kWarning)
00783 << "num_neighbours not found. Defaulting to 400" << endl;
00784 neighbours = 400;
00785 }
00786
00787 MSG("NuReco", Msg::kInfo) << "Done initializing shower energy kNN" << endl;
00788
00789 return true;
00790 }
|
|
||||||||||||
|
Definition at line 2568 of file NuReco.cxx. References abs(), NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpSRTrack::index, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, NtpStRecord::mc, NtpTHTrack::neumc, NtpMCStdHep::p4, NtpMCStdHep::parent, NtpMCTruth::stdhep, NtpStRecord::stdhep, and NtpStRecord::thtrk. Referenced by NuCounter::CountTrkStdhepId(). 02570 {
02571 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Event TCA
02572 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02573 if (numthtrks<=0) {
02574 MAXMSG("NuReco",Msg::kDebug,1)
02575 <<"No THTracks, so can't do IsDimuon..."<<endl;
02576 return false;
02577 }
02578 MAXMSG("NuReco",Msg::kDebug,1)
02579 <<"Found THTrack, IsDimuon..."<<endl;
02580 const NtpTHTrack& thtrk=
02581 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02582
02583 //now get the mc object (neutrino interaction) that
02584 //corresponds to the trk using the thtrk.neumc index
02585 TClonesArray& mcTca=(*ntp.mc);
02586 const NtpMCTruth& mc=
02587 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02588
02589 TClonesArray& hepTca=(*ntp.stdhep);
02590 //const Int_t numHeps=hepTca.GetEntriesFast();
02591
02592 Int_t charmEvent=-1;
02593 Bool_t isDimuon=false;
02594
02595 Int_t child0=-999;
02596 Int_t child1=-999;
02597
02598 //loop over stdhep
02599 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02600 const NtpMCStdHep& stdhep=
02601 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02602
02603 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02604 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02605 charmEvent=stdhep.mc;
02606 MAXMSG("NuReco",Msg::kDebug,3000)
02607 <<"Found Charm particle, mc="<<charmEvent
02608 <<", ihep="<<ihep
02609 <<", id="<<stdhep.IdHEP
02610 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02611 <<endl;
02612
02613 child0=stdhep.child[0];
02614 child1=stdhep.child[1];
02615 }
02616 else if ((abs(stdhep.IdHEP)>4000 && abs(stdhep.IdHEP)<5000) ||
02617 (abs(stdhep.IdHEP)>400 && abs(stdhep.IdHEP)<500)) {
02618 MAXMSG("NuReco",Msg::kInfo,3000)
02619 <<endl<<endl
02620 <<"What is this particle??? Charm?, mc="<<charmEvent
02621 <<", ihep="<<ihep
02622 <<", id="<<stdhep.IdHEP
02623 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02624 <<endl<<endl<<endl;
02625 }
02626
02627 if ((Int_t) stdhep.index>=child0 && (Int_t) stdhep.index<=child1){
02628 MAXMSG("NuReco",Msg::kDebug,3000)
02629 <<"Found Charm child..."
02630 <<", i="<<stdhep.index
02631 <<", child0="<<child0
02632 <<", child1="<<child1
02633 <<", id="<<stdhep.IdHEP
02634 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02635 <<endl;
02636 if (abs(stdhep.IdHEP)==13){
02637 isDimuon=true;
02638 MAXMSG("NuReco",Msg::kDebug,3000)
02639 <<"It's a muon!!!"
02640 <<" m="<<stdhep.mass
02641 <<", E="<<stdhep.p4[3]<<endl<<endl;
02642 break;
02643 }
02644 }
02645 }
02646
02647 if (isDimuon) return true;
02648 else return false;
02649 }
|
|
||||||||||||
|
this doesn't actually check that the tracked hits are from a muon it just checks that there was a muon in the THTrack neumc event Definition at line 2427 of file NuReco.cxx. References abs(), NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpSRTrack::index, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, NtpStRecord::mc, MSG, NtpTHTrack::neumc, NtpMCStdHep::p4, NtpMCStdHep::parent, PrintStdhep(), NtpMCTruth::stdhep, NtpStRecord::stdhep, and NtpStRecord::thtrk. Referenced by NuCounter::CountTrkStdhepId(). 02429 {
02432
02433 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Event TCA
02434 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02435 if (numthtrks<=0) {
02436 MAXMSG("NuReco",Msg::kDebug,1)
02437 <<"No THTracks, so can't do IsTrkWithMuonFromKaonDecay..."<<endl;
02438 return false;
02439 }
02440 MAXMSG("NuReco",Msg::kDebug,1)
02441 <<"Found THTrack, IsTrkWithMuonFromKaonDecay..."<<endl;
02442 const NtpTHTrack& thtrk=
02443 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02444
02445 //now get the mc object (neutrino interaction) that
02446 //corresponds to the trk using the thtrk.neumc index
02447 TClonesArray& mcTca=(*ntp.mc);
02448 const NtpMCTruth& mc=
02449 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02450
02451 TClonesArray& hepTca=(*ntp.stdhep);
02452 //const Int_t numHeps=hepTca.GetEntriesFast();
02453
02454 Int_t muonEvent=-1;
02455 Bool_t isFromKaon=false;
02456
02457 Int_t parent0=-999;
02458 Int_t parent1=-999;
02459
02460 Int_t parent0b=-999;
02461 Int_t parent1b=-999;
02462
02463 Int_t parent0c=-999;
02464 Int_t parent1c=-999;
02465
02466 //loop over stdhep
02467 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02468 const NtpMCStdHep& stdhep=
02469 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02470
02471 if (abs(stdhep.IdHEP)==13){
02472 muonEvent=stdhep.mc;
02473 MAXMSG("NuReco",Msg::kDebug,3000)
02474 <<"Found muon, mc="<<muonEvent
02475 <<", ihep="<<ihep
02476 <<", id="<<stdhep.IdHEP
02477 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02478 <<", parent=["<<stdhep.parent[0]
02479 <<","<<stdhep.parent[1]<<"]"
02480 <<", E="<<stdhep.p4[3]
02481 <<endl;
02482
02483 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02484 MSG("NuReco",Msg::kWarning)
02485 <<"Ahhhhhhhh, multiple parents of a muon"
02486 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02487 this->PrintStdhep(ntp,mc);
02488 }
02489
02490 if (parent0==-999) {
02491 parent0=stdhep.parent[0];
02492 parent1=stdhep.parent[1];
02493 }
02494 else if (parent0!=-999 && parent0b==-999) {//check for second muon
02495 MAXMSG("NuReco",Msg::kDebug,3000)
02496 <<"Found second muon, 1st parent=["
02497 <<parent0<<","<<parent1<<"]"
02498 <<", 2nd parent=["<<stdhep.parent[0]
02499 <<","<<stdhep.parent[1]<<"]"
02500 <<endl;
02501 parent0b=stdhep.parent[0];
02502 parent1b=stdhep.parent[1];
02503 }
02504 else if (parent0b!=-999 && parent0c==-999) {//check for third muon
02505 MAXMSG("NuReco",Msg::kDebug,3000)
02506 <<endl
02507 <<"Found third muon, parentb=["
02508 <<parent0b<<","<<parent1b<<"]"
02509 <<", new parent=["<<stdhep.parent[0]
02510 <<","<<stdhep.parent[1]<<"]"
02511 <<endl;
02512 parent0c=stdhep.parent[0];
02513 parent1c=stdhep.parent[1];
02514 }
02515 else if (parent0c!=-999) {//check for fourth muon
02516 MAXMSG("NuReco",Msg::kDebug,3000)
02517 <<endl
02518 <<"Ahhh, already found >=3 muons, parentb=["
02519 <<parent0b<<","<<parent1b<<"]"
02520 <<", new parent=["<<stdhep.parent[0]
02521 <<","<<stdhep.parent[1]<<"]"
02522 <<endl;
02523 //this->PrintStdhep(ntp,mc);
02524 }
02525 }
02526 }
02527
02528 //have to have a second loop since the parents will be before muon!
02529 //loop over stdhep (again)
02530 //don't have to do this - could jump about within stdhep in above loop
02531 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02532 const NtpMCStdHep& stdhep=
02533 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02534
02535 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02536 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02537
02538 if (abs(stdhep.IdHEP)!=14) {
02539 MAXMSG("NuReco",Msg::kDebug,3000)
02540 <<endl
02541 <<"**** Found non-nu muon parent..."
02542 <<", i="<<stdhep.index
02543 <<", id="<<stdhep.IdHEP
02544 <<", parent0="<<parent0
02545 <<", parent1="<<parent1
02546 <<", parent0b="<<parent0b
02547 <<", parent1b="<<parent1b
02548 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02549 <<endl;
02550 if (abs(stdhep.IdHEP)>300 && abs(stdhep.IdHEP)<400){
02551 isFromKaon=true;
02552 MAXMSG("NuReco",Msg::kDebug,3000)
02553 <<"From a kaon!!!"
02554 <<" m="<<stdhep.mass
02555 <<", E="<<stdhep.p4[3]<<endl<<endl;
02556 break;
02557 }
02558 }
02559 }
02560 }
02561
02562 if (isFromKaon) return true;
02563 else return false;
02564 }
|
|
||||||||||||
|
NOT IMPLEMENTED YET!!! Definition at line 2110 of file NuReco.cxx. References abs(), NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, MSG, NtpTHTrack::neumc, NtpMCStdHep::p4, NtpMCStdHep::parent, and NtpMCTruth::stdhep. Referenced by NuCounter::CountTrkStdhepId(). 02111 {
02113 return false;
02114
02115 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Event TCA
02116 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02117 if (numthtrks<=0) {
02118 MAXMSG("NuReco",Msg::kDebug,1)
02119 <<"No THTracks, so can't do IsMuonFromNotNuNotPiNotCharm..."
02120 <<endl;
02121 return false;
02122 }
02123 MAXMSG("NuReco",Msg::kDebug,1)
02124 <<"Found THTrack, IsMuonFromNotNuNotPiNotCharm..."<<endl;
02125 const NtpTHTrack& thtrk=
02126 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02127
02128 //now get the mc object (neutrino interaction) that
02129 //corresponds to the trk using the thtrk.neumc index
02130 TClonesArray& mcTca=(*ntp.mc);
02131 const NtpMCTruth& mc=
02132 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02133
02134 TClonesArray& hepTca=(*ntp.stdhep);
02135 //const Int_t numHeps=hepTca.GetEntriesFast();
02136
02137 Int_t muonEvent=-1;
02138 Bool_t isFromOther=false;
02139
02140 Int_t parent0=-999;
02141 Int_t parent1=-999;
02142
02143 Int_t parent0b=-999;
02144 Int_t parent1b=-999;
02145
02146 Int_t parent0c=-999;
02147 Int_t parent1c=-999;
02148
02149 //loop over stdhep
02150 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02151 const NtpMCStdHep& stdhep=
02152 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02153
02154 if (abs(stdhep.IdHEP)==13){
02155 muonEvent=stdhep.mc;
02156 MAXMSG("NuReco",Msg::kDebug,3000)
02157 <<"Found muon, mc="<<muonEvent
02158 <<", ihep="<<ihep
02159 <<", id="<<stdhep.IdHEP
02160 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02161 <<", parent=["<<stdhep.parent[0]
02162 <<","<<stdhep.parent[1]<<"]"
02163 <<", E="<<stdhep.p4[3]
02164 <<endl;
02165
02166 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02167 MSG("NuReco",Msg::kWarning)
02168 <<"Ahhhhhhhh, multiple parents of a muon"
02169 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02170 this->PrintStdhep(ntp,mc);
02171 }
02172
02173 if (parent0==-999) {
02174 parent0=stdhep.parent[0];
02175 parent1=stdhep.parent[1];
02176 }
02177 else if (parent0!=-999 && parent0b==-999) {//check for second muon
02178 MAXMSG("NuReco",Msg::kDebug,3000)
02179 <<"Found second muon, 1st parent=["
02180 <<parent0<<","<<parent1<<"]"
02181 <<", 2nd parent=["<<stdhep.parent[0]
02182 <<","<<stdhep.parent[1]<<"]"
02183 <<endl;
02184 parent0b=stdhep.parent[0];
02185 parent1b=stdhep.parent[1];
02186 }
02187 else if (parent0b!=-999 && parent0c==-999) {//check for third muon
02188 MAXMSG("NuReco",Msg::kDebug,3000)
02189 <<endl
02190 <<"Found third muon, parentb=["
02191 <<parent0b<<","<<parent1b<<"]"
02192 <<", new parent=["<<stdhep.parent[0]
02193 <<","<<stdhep.parent[1]<<"]"
02194 <<endl;
02195 parent0c=stdhep.parent[0];
02196 parent1c=stdhep.parent[1];
02197 //this->PrintStdhep(ntp,mc);
02198 }
02199 else if (parent0c!=-999) {//check for fourth muon
02200 MAXMSG("NuReco",Msg::kDebug,3000)
02201 <<endl
02202 <<"Ahhh, already found >=3 muons, parentb=["
02203 <<parent0b<<","<<parent1b<<"]"
02204 <<", new parent=["<<stdhep.parent[0]
02205 <<","<<stdhep.parent[1]<<"]"
02206 <<endl;
02207 //this->PrintStdhep(ntp,mc);
02208 }
02209 }
02210 }
02211
02212 //have to have a second loop since the parents will be before muon!
02213 //loop over stdhep (again)
02214 //don't have to do this - could jump about within stdhep in above loop
02215 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02216 const NtpMCStdHep& stdhep=
02217 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02218
02219 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02220 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02221
02222 if (abs(stdhep.IdHEP)!=14) {
02223 MAXMSG("NuReco",Msg::kDebug,3000)
02224 <<endl<<endl<<endl
02225 <<"**** Found non-nu muon parent..."
02226 <<", i="<<stdhep.index
02227 <<", id="<<stdhep.IdHEP
02228 <<", parent0="<<parent0
02229 <<", parent1="<<parent1
02230 <<", parent0b="<<parent0b
02231 <<", parent1b="<<parent1b
02232 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02233 <<endl<<endl<<endl;
02234 if (!((abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<400) ||
02235 (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02236 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122))){
02237 isFromOther=true;
02238 MAXMSG("NuReco",Msg::kInfo,3000)
02239 <<"From a ?? ?? !!!"
02240 <<" m="<<stdhep.mass
02241 <<", E="<<stdhep.p4[3]
02242 <<", id="<<stdhep.IdHEP
02243 <<endl<<endl;
02244
02245 this->PrintStdhep(ntp,mc);
02246 break;
02247 }
02248 }
02249 }
02250 }
02251
02252 if (isFromOther) return true;
02253 else return false;
02254 }
|
|
||||||||||||
|
this doesn't actually check that the tracked hits are from a muon it just checks that there was a muon in the THTrack neumc event Definition at line 2258 of file NuReco.cxx. References abs(), NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpSRTrack::index, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, NtpStRecord::mc, MSG, NtpTHTrack::neumc, NtpMCStdHep::p4, NtpMCStdHep::parent, PrintStdhep(), NtpMCTruth::stdhep, NtpStRecord::stdhep, and NtpStRecord::thtrk. Referenced by NuCounter::CountTrkStdhepId(). 02260 {
02263
02264 TClonesArray& thtrkTca=(*ntp.thtrk);//TruthHelper Event TCA
02265 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02266 if (numthtrks<=0) {
02267 MAXMSG("NuReco",Msg::kDebug,1)
02268 <<"No THTracks, so can't do IsTrkWithMuonFromPionDecay..."<<endl;
02269 return false;
02270 }
02271 MAXMSG("NuReco",Msg::kDebug,1)
02272 <<"Found THTrack, IsTrkWithMuonFromPionDecay..."<<endl;
02273 const NtpTHTrack& thtrk=
02274 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02275
02276 //now get the mc object (neutrino interaction) that
02277 //corresponds to the trk using the thtrk.neumc index
02278 TClonesArray& mcTca=(*ntp.mc);
02279 const NtpMCTruth& mc=
02280 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02281
02282 TClonesArray& hepTca=(*ntp.stdhep);
02283 //const Int_t numHeps=hepTca.GetEntriesFast();
02284
02285 Int_t muonEvent=-1;
02286 Bool_t isFromPion=false;
02287
02288 Int_t parent0=-999;
02289 Int_t parent1=-999;
02290
02291 Int_t parent0b=-999;
02292 Int_t parent1b=-999;
02293
02294 Int_t parent0c=-999;
02295 Int_t parent1c=-999;
02296
02297 //loop over stdhep
02298 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02299 const NtpMCStdHep& stdhep=
02300 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02301
02302 if (abs(stdhep.IdHEP)==13){
02303 muonEvent=stdhep.mc;
02304 MAXMSG("NuReco",Msg::kDebug,3000)
02305 <<"Found muon, mc="<<muonEvent
02306 <<", ihep="<<ihep
02307 <<", id="<<stdhep.IdHEP
02308 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02309 <<", parent=["<<stdhep.parent[0]
02310 <<","<<stdhep.parent[1]<<"]"
02311 <<", E="<<stdhep.p4[3]
02312 <<endl;
02313
02314 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02315 MSG("NuReco",Msg::kWarning)
02316 <<"Ahhhhhhhh, multiple parents of a muon"
02317 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02318 this->PrintStdhep(ntp,mc);
02319 }
02320
02321 if (parent0==-999) {
02322 parent0=stdhep.parent[0];
02323 parent1=stdhep.parent[1];
02324 }
02325 else if (parent0!=-999 && parent0b==-999) {//check for second muon
02326 MAXMSG("NuReco",Msg::kDebug,3000)
02327 <<"Found second muon, 1st parent=["
02328 <<parent0<<","<<parent1<<"]"
02329 <<", 2nd parent=["<<stdhep.parent[0]
02330 <<","<<stdhep.parent[1]<<"]"
02331 <<endl;
02332 parent0b=stdhep.parent[0];
02333 parent1b=stdhep.parent[1];
02334 }
02335 else if (parent0b!=-999 && parent0c==-999) {//check for third muon
02336 MAXMSG("NuReco",Msg::kDebug,3000)
02337 <<endl
02338 <<"Found third muon, parentb=["
02339 <<parent0b<<","<<parent1b<<"]"
02340 <<", new parent=["<<stdhep.parent[0]
02341 <<","<<stdhep.parent[1]<<"]"
02342 <<endl;
02343 parent0c=stdhep.parent[0];
02344 parent1c=stdhep.parent[1];
02345 }
02346 else if (parent0c!=-999) {//check for fourth muon
02347 MAXMSG("NuReco",Msg::kDebug,3000)
02348 <<endl
02349 <<"Ahhh, already found >=3 muons, parentb=["
02350 <<parent0b<<","<<parent1b<<"]"
02351 <<", new parent=["<<stdhep.parent[0]
02352 <<","<<stdhep.parent[1]<<"]"
02353 <<endl;
02354 //this->PrintStdhep(ntp,mc);
02355 }
02356 }
02357 }
02358
02359 //have to have a second loop since the parents will be before muon!
02360 //loop over stdhep (again)
02361 //don't have to do this - could jump about within stdhep in above loop
02362 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02363 const NtpMCStdHep& stdhep=
02364 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02365
02366 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02367 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02368
02369 if (abs(stdhep.IdHEP)!=14) {
02370 MAXMSG("NuReco",Msg::kDebug,3000)
02371 <<endl
02372 <<"**** Found non-nu muon parent..."
02373 <<", i="<<stdhep.index
02374 <<", id="<<stdhep.IdHEP
02375 <<", parent0="<<parent0
02376 <<", parent1="<<parent1
02377 <<", parent0b="<<parent0b
02378 <<", parent1b="<<parent1b
02379 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02380 <<endl;
02381 if (abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<300){
02382 isFromPion=true;
02383 MAXMSG("NuReco",Msg::kDebug,3000)
02384 <<"From a pion!!!"
02385 <<" m="<<stdhep.mass
02386 <<", E="<<stdhep.p4[3]<<endl<<endl;
02387 break;
02388 }
02389 }
02390 }
02391
02392 /* //this doesn't work since the decay child of the pion is not set
02393 //need to check the other way round too - event may be associated
02394 //with a pion which actually decayed to a muon fairly quickly
02395 if (abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<400){
02396
02397 if (stdhep.child[0]>0 && stdhep.child[1]>0) {
02398 const NtpMCStdHep& stdhepChild0=
02399 *dynamic_cast<NtpMCStdHep*>(hepTca[stdhep.child[0]]);
02400 const NtpMCStdHep& stdhepChild1=
02401 *dynamic_cast<NtpMCStdHep*>(hepTca[stdhep.child[1]]);
02402
02403 if (abs(stdhepChild0.IdHEP)==13 ||
02404 abs(stdhepChild1.IdHEP)==13) {
02405 MAXMSG("NuReco",Msg::kDebug,3000)
02406 <<endl
02407 <<"Found pion/kaon with muon child..."
02408 <<", i="<<stdhep.index
02409 <<", id="<<stdhep.IdHEP
02410 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02411 <<", child.IdHEP=["<<stdhepChild0.IdHEP
02412 <<","<<stdhepChild1.IdHEP<<"]"
02413 <<endl;
02414 }
02415 }
02416 }
02417 */
02418
02419 }
02420
02421 if (isFromPion) return true;
02422 else return false;
02423 }
|
|
||||||||||||
|
Definition at line 2653 of file NuReco.cxx. References abs(), NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpSREvent::index, NtpMCStdHep::IstHEP, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, NtpStRecord::mc, NtpTHEvent::neumc, NtpMCStdHep::p4, NtpMCStdHep::parent, NtpMCTruth::stdhep, NtpStRecord::stdhep, and NtpStRecord::thevt. 02655 {
02656 TClonesArray& thevtTca=(*ntp.thevt);//TruthHelper Event TCA
02657 const Int_t numthevts=thevtTca.GetEntriesFast();
02658 if (numthevts<=0) {
02659 MAXMSG("NuReco",Msg::kDebug,1)
02660 <<"No THEvents, so can't GetTruthInfo..."<<endl;
02661 return false;
02662 }
02663 MAXMSG("NuReco",Msg::kDebug,1)
02664 <<"Found THEvent, GetTruthInfo..."<<endl;
02665 const NtpTHEvent& thevt=
02666 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02667
02668 //now get the mc object (neutrino interaction) that
02669 //corresponds to the trk using the thtrk.neumc index
02670 TClonesArray& mcTca=(*ntp.mc);
02671 const NtpMCTruth& mc=
02672 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
02673
02674 TClonesArray& hepTca=(*ntp.stdhep);
02675 //const Int_t numHeps=hepTca.GetEntriesFast();
02676
02677 Int_t charmEvent=-1;
02678
02679 Int_t child0=-999;
02680 Int_t child1=-999;
02681
02682 //loop over stdhep
02683 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02684 const NtpMCStdHep& stdhep=
02685 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02686
02687 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02688 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02689 charmEvent=stdhep.mc;
02690 MAXMSG("NuReco",Msg::kDebug,3000)
02691 <<"Found Charm particle, mc="<<charmEvent
02692 <<", ihep="<<ihep
02693 <<", id="<<stdhep.IdHEP
02694 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02695 <<endl;
02696
02697 child0=stdhep.child[0];
02698 child1=stdhep.child[1];
02699 }
02700
02701 if (child0==(Int_t)stdhep.index || child1==(Int_t)stdhep.index) {
02702 MAXMSG("NuReco",Msg::kDebug,3000)
02703 <<"Found Charm child..."
02704 <<", i="<<stdhep.index
02705 <<", child0="<<child0
02706 <<", child1="<<child1
02707 <<", id="<<stdhep.IdHEP
02708 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02709 <<endl;
02710 if (abs(stdhep.IdHEP)==13){
02711 MAXMSG("NuReco",Msg::kDebug,3000)
02712 <<"It's a muon!!!"
02713 <<" m="<<stdhep.mass
02714 <<", E="<<stdhep.p4[3]<<endl<<endl;
02715 }
02716 }
02717
02718 }
02719
02720 //loop over stdhep
02721 if (charmEvent!=-1){
02722 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02723 const NtpMCStdHep& stdhep=
02724 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02725
02726 MAXMSG("NuReco",Msg::kInfo,3000)
02727 <<"ihep="<<ihep
02728 <<", mc="<<stdhep.mc
02729 <<", id="<<stdhep.IdHEP
02730 <<", Ist="<<stdhep.IstHEP
02731 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02732 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02733 <<", m="<<stdhep.mass
02734 <<", E="<<stdhep.p4[3]
02735 <<endl;
02736 }
02737 return true;
02738 }
02739 else return false;
02740 }
|
|
||||||||||||
|
Definition at line 2744 of file NuReco.cxx. References NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, NtpMCStdHep::mass, MAXMSG, NtpMCStdHep::mc, NtpMCStdHep::p4, NtpMCStdHep::parent, NtpMCTruth::stdhep, and NtpStRecord::stdhep. Referenced by NuCounter::CountTrkStdhepId(), GetPrimaryMuonStdHepIndex(), IsTrkWithMuonFromKaonDecay(), and IsTrkWithMuonFromPionDecay(). 02746 {
02747 TClonesArray& hepTca=(*ntp.stdhep);
02748 //const Int_t numHeps=hepTca.GetEntriesFast();
02749
02750 MAXMSG("NuReco",Msg::kInfo,3000)
02751 <<"PrintStdhep: mc.stdhep=["<<mc.stdhep[0]
02752 <<","<<mc.stdhep[1]<<"]"<<endl;
02753
02754 //loop over stdhep
02755 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02756 const NtpMCStdHep& stdhep=
02757 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02758
02759 MAXMSG("NuReco",Msg::kInfo,3000)
02760 <<"ihep="<<ihep
02761 <<", mc="<<stdhep.mc
02762 <<", id="<<stdhep.IdHEP
02763 <<", Ist="<<stdhep.IstHEP//Interaction status code, from Geant+200
02764 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02765 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02766 <<", m="<<stdhep.mass
02767 <<", E="<<stdhep.p4[3]
02768 <<endl;
02769 }
02770 }
|
|
||||||||||||||||
|
Definition at line 3331 of file NuReco.cxx. References NtpMCTruth::iaction, NtpTHEvent::index, NtpMCTruth::index, NtpSREvent::index, MAXMSG, NtpStRecord::mc, NtpTHEvent::neumc, NtpMCTruth::p4neu, NtpStRecord::thevt, and NtpMCTruth::y. 03333 {
03334 TClonesArray& thevtTca=(*ntp.thevt);//TruthHelper Event TCA
03335 const Int_t numthevts=thevtTca.GetEntriesFast();
03336 if (numthevts<=0) {
03337 MAXMSG("NuReco",Msg::kInfo,1)
03338 <<"No THEvents, so can't print truth information..."<<endl;
03339 MAXMSG("NuReco",Msg::kInfo,20)
03340 <<"recoEn="<<recoEn<<endl;
03341 return;
03342 }
03343 MAXMSG("NuReco",Msg::kInfo,1)
03344 <<"Found THEvent, printing info..."<<endl;
03345 const NtpTHEvent& thevt=
03346 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
03347
03348 //now get the mc object (neutrino interaction)
03349 TClonesArray& mcTca=(*ntp.mc);
03350 const NtpMCTruth& mc=
03351 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
03352 MAXMSG("NuReco",Msg::kDebug,100)
03353 <<"Found mc="<<mc.index<<", thevt.index="<<thevt.index<<endl;
03354
03355 Double_t nuEn=mc.p4neu[3];
03356
03357 Double_t trueEn=mc.y*mc.p4neu[3];//NC
03358 if (mc.iaction==1) trueEn=mc.p4neu[3];//CC
03359
03360 Double_t muEn=(1-mc.y)*mc.p4neu[3];
03361 Double_t hadEn=mc.y*mc.p4neu[3];
03362
03363 string siaction="NC";
03364 if (mc.iaction==1) siaction="CC";
03365
03366 if (mc.iaction==1) {
03367 MAXMSG("NuReco",Msg::kInfo,100)
03368 <<"RecoEn="<<recoEn
03369 <<", "<<siaction
03370 <<" truth: En="<<trueEn<<", mu="<<muEn<<", had="<<hadEn
03371 <<", (y="<<mc.y<<")"
03372 <<endl;
03373 }
03374 else {
03375 MAXMSG("NuReco",Msg::kInfo,20)
03376 <<"RecoEn="<<recoEn
03377 <<", "<<siaction
03378 <<" truth: had="<<hadEn
03379 <<", (Nu="<<nuEn<<", y="<<mc.y<<")"
03380 <<endl;
03381 }
03382 }
|
|
|
|
|
|
Definition at line 1551 of file NuReco.cxx. References GetBestTrack(), NuLibrary::Instance(), NuEvent::jitter, NuEvent::jitter1, NuEvent::jitter2, NuEvent::jitter3, NuEvent::jPID, NuEvent::jPID1, NuEvent::jPID2, NuEvent::jPID3, NuEvent::majC, NuEvent::majC1, NuEvent::majC2, NuEvent::majC3, MAXMSG, NuEvent::ntrk, NuEvent::primtrk, NuLibrary::reco, NuEvent::smoothMajC, NuEvent::smoothMajC1, NuEvent::smoothMajC2, NuEvent::smoothMajC3, NuEvent::trkExists1, NuEvent::trkExists2, and NuEvent::trkExists3. Referenced by NuExtraction::ExtractMajorityCurvature(), and SetBestTrk(). 01552 {
01553 //get an instance of the code library
01554 const NuLibrary& lib=NuLibrary::Instance();
01555
01556 //get the best track
01557 Int_t bestTrack=lib.reco.GetBestTrack(nu);
01558
01559 //copy the best track across
01560 if (bestTrack<1) {
01561 MAXMSG("NuReco",Msg::kInfo,3)
01562 <<"SetBestTrkMajorityCurvature: No best trk, primtrk="
01563 <<nu.primtrk
01564 <<", ntrk="<<nu.ntrk
01565 <<", trkExists1,2,3="<<nu.trkExists1<<","
01566 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01567 return;
01568 }
01569 else if (bestTrack==1) {
01570 MAXMSG("NuReco",Msg::kDebug,3)
01571 <<"SetBestTrkMajorityCurvature: Best trk=first, primtrk="
01572 <<nu.primtrk
01573 <<", ntrk="<<nu.ntrk
01574 <<", trkExists1,2,3="<<nu.trkExists1<<","
01575 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01576
01577 nu.jitter=nu.jitter1;
01578 nu.jPID=nu.jPID1;
01579 nu.majC=nu.majC1;
01580 //nu.majCRatio=nu.majCRatio1;
01581 //nu.rms=nu.rms1;
01582 //nu.simpleMajC=nu.simpleMajC1;
01583 nu.smoothMajC=nu.smoothMajC1;
01584 //nu.sqJitter=nu.sqJitter1;
01585 //nu.totWidth=nu.totWidth1;
01586 }
01587 else if (bestTrack==2) {
01588 MAXMSG("NuReco",Msg::kDebug,3)
01589 <<"SetBestTrkMajorityCurvature: Best trk=second, primtrk="
01590 <<nu.primtrk
01591 <<", ntrk="<<nu.ntrk
01592 <<", trkExists1,2,3="<<nu.trkExists1<<","
01593 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01594
01595 nu.jitter=nu.jitter2;
01596 nu.jPID=nu.jPID2;
01597 nu.majC=nu.majC2;
01598 //nu.majCRatio=nu.majCRatio2;
01599 //nu.rms=nu.rms2;
01600 //nu.simpleMajC=nu.simpleMajC2;
01601 nu.smoothMajC=nu.smoothMajC2;
01602 //nu.sqJitter=nu.sqJitter2;
01603 //nu.totWidth=nu.totWidth2;
01604 }
01605 else if (bestTrack==3) {
01606 MAXMSG("NuReco",Msg::kDebug,3)
01607 <<"SetBestTrkMajorityCurvature: Best trk=third, primtrk="
01608 <<nu.primtrk
01609 <<", ntrk="<<nu.ntrk
01610 <<", trkExists1,2,3="<<nu.trkExists1<<","
01611 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01612
01613 nu.jitter=nu.jitter3;
01614 nu.jPID=nu.jPID3;
01615 nu.majC=nu.majC3;
01616 //nu.majCRatio=nu.majCRatio3;
01617 //nu.rms=nu.rms3;
01618 //nu.simpleMajC=nu.simpleMajC3;
01619 nu.smoothMajC=nu.smoothMajC3;
01620 //nu.sqJitter=nu.sqJitter3;
01621 //nu.totWidth=nu.totWidth3;
01622 }
01623 else cout<<"Ahhhhhhhhhhhh"<<endl;
01624 }
|
|
|
|
Definition at line 4561 of file NuReco.cxx. References GetSmallestDeepDistToEdge(), NuEvent::xEvtVtx, and NuEvent::yEvtVtx. 04562 {
04563 NuEvent nu;
04564
04565 TH2F* hYvsXDistToEdge=new TH2F
04566 ("hYvsXDistToEdge","hYvsXDistToEdge",
04567 600,-3,3, 600,-3,3);
04568 hYvsXDistToEdge->SetTitle("SmallestDeepDistToEdge");
04569 hYvsXDistToEdge->GetXaxis()->SetTitle("X (m)");
04570 hYvsXDistToEdge->GetXaxis()->CenterTitle();
04571 hYvsXDistToEdge->GetYaxis()->SetTitle("Y (m)");
04572 hYvsXDistToEdge->GetYaxis()->CenterTitle();
04573
04574 for (Float_t x=-3;x<3;x+=0.01) {
04575 for (Float_t y=-3;y<3;y+=0.01) {
04576
04577 nu.xEvtVtx=x;
04578 nu.yEvtVtx=y;
04579 Float_t dist=this->GetSmallestDeepDistToEdge(nu);
04580 //cout<<"x,y="<<x<<","<<y<<", dist="<<dist<<endl;
04581 hYvsXDistToEdge->Fill(x,y,dist);
04582 }
04583 }
04584 }
|
|
|
Definition at line 1945 of file NuReco.cxx. References NuEvent::containmentFlag. Referenced by CalcTrkVariables(), and GetTrackEnergy(). 01946 {
01947 if (nu.containmentFlag==2 ||
01948 nu.containmentFlag==4) return true;
01949 else return false;
01950 }
|
|
|
Definition at line 1936 of file NuReco.cxx. References NuEvent::containmentFlag. Referenced by CalcTrkVariables(), and GetTrackEnergy(). 01937 {
01938 if (nu.containmentFlag==1 ||
01939 nu.containmentFlag==3) return true;
01940 else return false;
01941 }
|
|
||||||||||||||||||||
|
Definition at line 4602 of file NuReco.cxx. References UgliGeomHandle::uvz2xyz(). Referenced by NuExtraction::ExtractSAFitInfo(), and NuPlots::FillDetectorEdge(). 04604 {
04605 //get an ugh
04606 UgliGeomHandle ugh(vc);
04607
04608 //calculate the positions in XYZ system
04609 TVector3 uvz(u,v,z);
04610 TVector3 xyz=ugh.uvz2xyz(uvz);
04611 return xyz;
04612 }
|
|
||||||||||||||||||||
|
Definition at line 4588 of file NuReco.cxx. References UgliGeomHandle::xyz2uvz(). Referenced by CalcExtraTruthVariables(). 04590 {
04591 //get an ugh
04592 UgliGeomHandle ugh(vc);
04593
04594 //calculate the positions in UVZ system
04595 TVector3 xyz(x,y,z);
04596 TVector3 uvz=ugh.xyz2uvz(xyz);
04597 return uvz;
04598 }
|
1.3.9.1