diff --git a/fcl/CrvExpert.fcl b/fcl/CrvExpert.fcl index 0a10def..76f3cd9 100644 --- a/fcl/CrvExpert.fcl +++ b/fcl/CrvExpert.fcl @@ -1,7 +1,3 @@ #include "TrkAna/fcl/TrkAnaReco.fcl" - -physics.analyzers.TrkAnaNeg.FillCRVHits : true physics.analyzers.TrkAnaNeg.FillCRVPulses : false - -physics.analyzers.TrkAnaPos.FillCRVHits : true -physics.analyzers.TrkAnaPos.FillCRVPulses : false \ No newline at end of file +physics.analyzers.TrkAnaPos.FillCRVPulses : false diff --git a/fcl/TrkAnaExtracted.fcl b/fcl/TrkAnaExtracted.fcl index e532fba..e8d5d62 100644 --- a/fcl/TrkAnaExtracted.fcl +++ b/fcl/TrkAnaExtracted.fcl @@ -16,7 +16,7 @@ physics : } analyzers : { @table::TrkAnaReco.analyzers } } -physics.TrkAnaTrigPath : [ BestCrvLine ] +physics.TrkAnaTrigPath : [ ] physics.TrkAnaEndPath : [ TrkAnaExt ] # Include more information (MC, full TrkQual and TrkPID branches) @@ -25,12 +25,10 @@ physics.analyzers.TrkAnaExt.candidate.options.fillTrkQual : false physics.analyzers.TrkAnaExt.candidate.options.fillTrkPID : false physics.analyzers.TrkAnaExt.FillTriggerInfo : false physics.analyzers.TrkAnaExt.FitType : KinematicLine -physics.analyzers.TrkAnaExt.FillCRVHits : true -physics.analyzers.TrkAnaExt.FillCRVPulses : true # for hit level diagnostics, set diagLevel to 2 physics.analyzers.TrkAnaExt.diagLevel : 2 physics.analyzers.TrkAnaExt.FillMCInfo : true -physics.analyzers.TrkAnaExt.FillCRV : true +physics.analyzers.TrkAnaExt.FillCRVHits : true services.TFileService.fileName: "nts.owner.trkanaextracted-reco.version.sequencer.root" diff --git a/fcl/TrkAnaReco.fcl b/fcl/TrkAnaReco.fcl index 07977c1..9722b6a 100644 --- a/fcl/TrkAnaReco.fcl +++ b/fcl/TrkAnaReco.fcl @@ -14,13 +14,11 @@ physics : producers : { # @table::TrkAnaReco.producers PBIWeight : @local::PBIWeight - @table::BestCrvProducers } analyzers : { @table::TrkAnaReco.analyzers } } -#physics.TrkAnaTrigPath : [ @sequence::TrkAnaReco.MCTrigSequence, @sequence::TrkAnaReco.TrigSequence ] -physics.TrkAnaTrigPath : [ PBIWeight, @sequence::BestCrvProducersPath ] +physics.TrkAnaTrigPath : [ PBIWeight ] physics.TrkAnaEndPath : [ @sequence::TrkAnaReco.EndSequence ] # Include more information (MC, full TrkQual and TrkPID branches) @@ -30,13 +28,13 @@ physics.analyzers.TrkAnaPos.candidate.options : @local::AllOpt # for hit level diagnostics, set diagLevel to 2 physics.analyzers.TrkAnaNeg.diagLevel : 1 physics.analyzers.TrkAnaNeg.FillMCInfo : true -physics.analyzers.TrkAnaNeg.FillCRV : true +physics.analyzers.TrkAnaNeg.FillCRVHits : true physics.analyzers.TrkAnaNeg.FillTriggerInfo : false physics.analyzers.TrkAnaNeg.FitType : LoopHelix physics.analyzers.TrkAnaPos.diagLevel : 1 physics.analyzers.TrkAnaPos.FillMCInfo : true -physics.analyzers.TrkAnaPos.FillCRV : true +physics.analyzers.TrkAnaPos.FillCRVHits : true physics.analyzers.TrkAnaPos.FillTriggerInfo : false physics.analyzers.TrkAnaPos.FitType : LoopHelix diff --git a/fcl/TrkAnaReco_wTrkQualFilter.fcl b/fcl/TrkAnaReco_wTrkQualFilter.fcl index 518c334..b095c19 100644 --- a/fcl/TrkAnaReco_wTrkQualFilter.fcl +++ b/fcl/TrkAnaReco_wTrkQualFilter.fcl @@ -13,7 +13,7 @@ physics : { producers : { @table::TrkAnaReco.producers } analyzers : { @table::TrkAnaReco.analyzers } - filters : { + filters : { trkQualFilter : { module_type : TrkQualFilter trainName : "TrkQual" trkHypo : "DeM" @@ -53,11 +53,11 @@ physics.analyzers.TrkAnaPos.candidate.options : @local::AllOpt # for hit level diagnostics, set diagLevel to 2 physics.analyzers.TrkAnaNeg.diagLevel : 1 physics.analyzers.TrkAnaNeg.FillMCInfo : true -physics.analyzers.TrkAnaNeg.FillCRV : true +physics.analyzers.TrkAnaNeg.FillCRVHits : true physics.analyzers.TrkAnaNeg.FillTriggerInfo : true physics.analyzers.TrkAnaPos.diagLevel : 1 physics.analyzers.TrkAnaPos.FillMCInfo : true -physics.analyzers.TrkAnaPos.FillCRV : true +physics.analyzers.TrkAnaPos.FillCRVHits : true physics.analyzers.TrkAnaPos.FillTriggerInfo : true services.TFileService.fileName: "nts.owner.trkana-reco.version.sequencer.root" diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index e66fd8a..9d072e8 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -67,43 +67,6 @@ TrkPIDProducers : { } TrkPIDProducersPath : [ TrkPIDDeM, TrkPIDDeP ] -BestCrv : { - module_type : BestCrvHitDeltaT - crvCoincidenceTag : "SelectRecoMC:CrvCoincidenceClusterFinder" -} - -BestCrvDeM : @local::BestCrv -BestCrvDeM.kalSeedTag : "KKDeM" -BestCrvUeM : @local::BestCrv -BestCrvUeM.kalSeedTag : "KKUeM" -BestCrvDmuM : @local::BestCrv -BestCrvDmuM.kalSeedTag : "KKDmuM" -BestCrvUmuM : @local::BestCrv -BestCrvUmuM.kalSeedTag : "KKUmuM" -BestCrvDeP : @local::BestCrv -BestCrvDeP.kalSeedTag : "KKDeP" -BestCrvUeP : @local::BestCrv -BestCrvUeP.kalSeedTag : "KKUeP" -BestCrvDmuP : @local::BestCrv -BestCrvDmuP.kalSeedTag : "KKDmuP" -BestCrvUmuP : @local::BestCrv -BestCrvUmuP.kalSeedTag : "KKUmuP" -BestCrvLine : @local::BestCrv -BestCrvLine.kalSeedTag : "KKLine" - -BestCrvProducers : { - BestCrvDeM : @local::BestCrvDeM - BestCrvUeM : @local::BestCrvUeM - BestCrvDeP : @local::BestCrvDeP - BestCrvUeP : @local::BestCrvUeP - BestCrvDmuM : @local::BestCrvDmuM - BestCrvDmuP : @local::BestCrvDmuP - BestCrvUmuM : @local::BestCrvUmuM - BestCrvUmuP : @local::BestCrvUmuP - BestCrvLine : @local::BestCrvLine -} -BestCrvProducersPath : [ BestCrvDeM, BestCrvUeM, BestCrvDmuM, BestCrvDeP, BestCrvUeP, BestCrvDmuP, BestCrvUmuM, BestCrvUmuP ] - # DIO weighting for flat spectrum electrons DIOWeight: { module_type: DecayInOrbitWeight @@ -141,11 +104,6 @@ AllOpt : { fillMC : true fillTrkPID : true fillHits : true genealogyDepth : -1 - fillBestCrv : true - - bestCrvModules : [ "BestCrv" ] - bestCrvInstances : [ "first" ] - bestCrvBranches : [ "bestcrv" ] } DeM : { input : "KK" @@ -192,14 +150,14 @@ TrkAnaTreeMaker : { PBITag : "PBISim" PBTTag : "EWMProducer" PBTMCTag : "EWMProducer" - CrvCoincidenceModuleLabel : "SelectRecoMC:CrvCoincidenceClusterFinder" - CrvCoincidenceMCModuleLabel : "compressRecoMCs:CrvCoincidenceClusterMatchMC" - CrvRecoPulseLabel : "SelectRecoMC" - CrvStepLabel : "compressRecoMCs" - SimParticleLabel : "compressRecoMCs" - MCTrajectoryLabel : "compressRecoMCs" - CrvWaveformsModuleLabel : "compressRecoMCs" - CrvDigiModuleLabel : "SelectRecoMC" + CrvCoincidencesTag : "SelectRecoMC:CrvCoincidenceClusterFinder" + CrvCoincidenceMCsTag : "compressRecoMCs:CrvCoincidenceClusterMatchMC" + CrvRecoPulsesTag : "SelectRecoMC" + CrvStepsTag : "compressRecoMCs" + SimParticlesTag : "compressRecoMCs" + MCTrajectoriesTag : "compressRecoMCs" + CrvDigiMCsTag : "compressRecoMCs" + CrvDigisTag : "SelectRecoMC" CrvCoincidenceClusterMCAssnsTag : "CrvCoincidenceClusterMCAssns" CrvPlaneY : 2653 FillMCInfo : true @@ -209,8 +167,7 @@ TrkAnaTreeMaker : { FillTriggerInfo : true TriggerProcessName : "Mix" ProcessEmptyEvents : false - FillCRV : true - FillCRVHits : false + FillCRVHits : true FillCRVPulses : false PrimaryParticleTag : "compressRecoMCs" KalSeedMCAssns : "SelectRecoMC" @@ -228,7 +185,6 @@ TrkAnaReco : { PBIWeight : @local::PBIWeight @table::TrkQualProducers @table::TrkPIDProducers - @table::BestCrvProducers } analyzers : { @@ -251,7 +207,7 @@ TrkAnaReco : { genCountLogger : @local::genCountLogger } - TrigSequence : [ PBIWeight, @sequence::TrkQualProducersPath, @sequence::TrkPIDProducersPath, @sequence::BestCrvProducersPath ] + TrigSequence : [ PBIWeight, @sequence::TrkQualProducersPath, @sequence::TrkPIDProducersPath ] EndSequenceNoMC : [ TrkAnaNeg, TrkAnaPos ] EndSequence : [ TrkAnaNeg, TrkAnaPos, genCountLogger ] diff --git a/fcl/prolog_trigger.fcl b/fcl/prolog_trigger.fcl index fab20e9..cf1ae1d 100644 --- a/fcl/prolog_trigger.fcl +++ b/fcl/prolog_trigger.fcl @@ -10,14 +10,13 @@ TrkAnaTrigger : { @table::DeM options : { @table::AllOpt - fillBestCrv : false fillHits : true } } supplements : [] ExtraMCStepCollectionTags : [] diagLevel : 2 - FillCRV : false + FillCRVHits : false FillMCInfo : true FillCaloMC : false RecoCountTag : "" @@ -28,8 +27,8 @@ TrkAnaTrigger : { FillTriggerInfo : false supplements : [] PrimaryParticleTag: "compressDigiMCs" - SimParticleLabel: "compressDigiMCs" - MCTrajectoryLabel: "compressDigiMCs" + SimParticlesTag: "compressDigiMCs" + MCTrajectoriesTag: "compressDigiMCs" ExtraMCStepCollectionTags: [ ] InfoMCStructHelper : { SimParticleCollectionTag : "compressDigiMCs" diff --git a/inc/CrvHitInfoMC.hh b/inc/CrvHitInfoMC.hh new file mode 100644 index 0000000..50c04c8 --- /dev/null +++ b/inc/CrvHitInfoMC.hh @@ -0,0 +1,53 @@ +#ifndef CrvHitInfoMC_hh +#define CrvHitInfoMC_hh + +#include "CLHEP/Vector/ThreeVector.h" +#include +#include "Rtypes.h" +#include "Offline/DataProducts/inc/GenVector.hh" + +namespace mu2e +{ + struct CrvHitInfoMC //information about the MC track which most likely caused the CRV coincidence triplets + { + Bool_t valid = false; //was an MC particle found that matches the coincidence triplets? + Int_t pdgId = -1; //PDG ID of this MC particle + Int_t primaryPdgId = -1; //PDG ID of the primary particle of this MC particle (helps to determine whether it was a cosmic ray, etc.) + Float_t primaryE = -1; //energy of the primary particle of this MC particle + XYZVectorF primary; //starting point of the primary particle of this MC particle + Int_t parentPdgId = -1; //PDG ID of the parent particle of this MC particle (helps to determine whether it was a cosmic ray, etc.) + Float_t parentE = -1; //energy of the parent particle of this MC particle + XYZVectorF parent; //starting point of the parent particle of this MC particle + Int_t gparentPdgId = -1; //PDG ID of the gparent particle of this MC particle (helps to determine whether it was a cosmic ray, etc.) + Float_t gparentE = -1; //energy of the gparent particle of this MC particle + XYZVectorF gparent; //starting point of the gparent particle of this MC particle + XYZVectorF pos; //position of the MC particle when it "created" the first StepPointMC + Float_t time = -1; //time of the MC particle when it "created" the first StepPointMC + Float_t depositedEnergy = -1; //total energy deposited for this cluster (not just for this track) + CrvHitInfoMC(){} + CrvHitInfoMC(bool valid, int pdgId, + int primaryPdgId, float primaryE, CLHEP::Hep3Vector primaryPos, + int parentPdgId, float parentE, CLHEP::Hep3Vector parentPos, + int gparentPdgId, float gparentE, CLHEP::Hep3Vector gparentPos, + CLHEP::Hep3Vector pos, float time, float depositedEnergy) : + valid(valid), + pdgId(pdgId), + primaryPdgId(primaryPdgId), + primaryE(primaryE), + primary(primaryPos), + parentPdgId(parentPdgId), + parentE(parentE), + parent(parentPos), + gparentPdgId(gparentPdgId), + gparentE(gparentE), + gparent(gparentPos), + pos(pos), + time(time), + depositedEnergy(depositedEnergy) + {} + }; + + typedef std::vector CrvHitInfoMCCollection; //this is the MC vector which will be stored in the main TTree + +} +#endif diff --git a/inc/CrvHitInfoReco.hh b/inc/CrvHitInfoReco.hh new file mode 100644 index 0000000..2d4789d --- /dev/null +++ b/inc/CrvHitInfoReco.hh @@ -0,0 +1,41 @@ +#ifndef CrvHitInfoReco_hh +#define CrvHitInfoReco_hh + +#include "Offline/DataProducts/inc/GenVector.hh" +#include "CLHEP/Vector/ThreeVector.h" +#include +#include "Rtypes.h" + +namespace mu2e +{ + + struct CrvHitInfoReco //information about a cluster of CRV coincidence triplets + { + Int_t sectorType =-1; //CRV sector type + XYZVectorF pos; //average position of counters + Float_t timeStart = -1; //first hit time + Float_t timeEnd = -1; //last hit time + Float_t time = -1; // average hit time + Int_t PEs = -1; //total number of PEs for this cluster + Int_t nHits = -1; //number of coincidence hits in this cluster + Int_t nLayers = -1; //number of coincidence layers in this cluster + Float_t angle = -999; //coincidence direction + + CrvHitInfoReco(){} + CrvHitInfoReco(int sectorType, CLHEP::Hep3Vector hpos, float timeWindowStart, float timeWindowEnd, float timeAvg, int PEs, int nCoincidenceHits, int nCoincidenceLayers, float coincidenceAngle) : + sectorType(sectorType), + pos(hpos), + timeStart(timeWindowStart), + timeEnd(timeWindowEnd), + time(timeAvg), + PEs(PEs), + nHits(nCoincidenceHits), + nLayers(nCoincidenceLayers), + angle(coincidenceAngle) + {} + }; + + typedef std::vector CrvHitInfoRecoCollection; //this is the reco vector which will be stored in the main TTree + +} +#endif diff --git a/inc/CrvInfoHelper.hh b/inc/CrvInfoHelper.hh new file mode 100644 index 0000000..3b0d469 --- /dev/null +++ b/inc/CrvInfoHelper.hh @@ -0,0 +1,67 @@ +// +// Replacement for CRVAnalysis +// +#include "TrkAna/inc/CrvHitInfoReco.hh" +#include "TrkAna/inc/CrvHitInfoMC.hh" +#include "TrkAna/inc/CrvWaveformInfo.hh" +#include "TrkAna/inc/CrvSummaryReco.hh" +#include "TrkAna/inc/CrvSummaryMC.hh" +#include "TrkAna/inc/CrvPlaneInfoMC.hh" +#include "TrkAna/inc/CrvPulseInfoReco.hh" +#include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" +#include "Offline/RecoDataProducts/inc/CrvRecoPulse.hh" +#include "Offline/RecoDataProducts/inc/CrvDigi.hh" +#include "Offline/MCDataProducts/inc/CrvDigiMC.hh" +#include "Offline/MCDataProducts/inc/CrvStep.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/CrvCoincidenceClusterMC.hh" +#include "Offline/MCDataProducts/inc/MCTrajectory.hh" +#include "art/Framework/Principal/Handle.h" + +namespace mu2e +{ + class CrvInfoHelper + { + public: + + CrvInfoHelper() {} + + void FillCrvHitInfoCollections( + art::Handle const& crvCoincidences, + art::Handle const& crvCoincidencesMC, + art::Handle const& crvRecoPulses, + art::Handle const& crvSteps, + art::Handle const& mcTrajectories, + CrvHitInfoRecoCollection &recoInfo, CrvHitInfoMCCollection &MCInfo, + CrvSummaryReco &recoSummary, CrvSummaryMC &MCSummary, + CrvPlaneInfoMCCollection &MCInfoPlane, double crvPlaneY); + + void FillCrvPulseInfoCollections( + art::Handle const& crvRecoPulses, + art::Handle const& crvDigiMCs, + art::Handle const& crvDigis, + CrvPulseInfoRecoCollection &recoInfo, CrvHitInfoMCCollection &MCInfo, CrvWaveformInfoCollection &waveformInfo); + + private: + static const art::Ptr &FindPrimaryParticle(const art::Ptr &simParticle) + { + return simParticle->hasParent() ? FindPrimaryParticle(simParticle->parent()) : simParticle; + } + static const SimParticle &FindPrimaryParticle(const SimParticle &simParticle) + { + return simParticle.hasParent() ? *FindPrimaryParticle(simParticle.parent()) : simParticle; + } + static const art::Ptr &FindParentParticle(const art::Ptr &simParticle) + { + return simParticle->hasParent() ? simParticle->parent() : simParticle; + } + static const art::Ptr &FindGParentParticle(const art::Ptr &simParticle) + { + const art::Ptr &parentParticle = simParticle->hasParent() ? simParticle->parent() : simParticle; + return parentParticle->hasParent() ? parentParticle->parent() : parentParticle; + } + + static const int _trajectorySimParticleId = 300001; //only temporarily here for some tests + }; + +} diff --git a/inc/CrvPlaneInfoMC.hh b/inc/CrvPlaneInfoMC.hh new file mode 100644 index 0000000..ad78841 --- /dev/null +++ b/inc/CrvPlaneInfoMC.hh @@ -0,0 +1,41 @@ +#ifndef CrvPlaneInfoMC_hh +#define CrvPlaneInfoMC_hh + +#include "CLHEP/Vector/ThreeVector.h" +#include +#include "Rtypes.h" + +namespace mu2e +{ + struct CrvPlaneInfoMC //information about the point where the MC trajectory crosses the xz plane of CRV-T + { + Int_t pdgId = -1; //PDG ID of this MC particle + Int_t primaryPdgId = -1; //PDG ID of the primary particle of this MC particle + Float_t primaryE = -1; //energy of the primary particle of this MC particle + XYZVectorF primary; //starting point of the primary particle of this MC particle + XYZVectorF pos; //position of the MC particle when it crosses the xz plane of CRV-T + XYZVectorF dir; //direction of the MC particle when it crosses the xz plane of CRV-T + Float_t time = -1; //time of the MC particle when it crosses the xz plane of CRV-T + Float_t kineticEnergy = -1; //time of the MC particle when it crosses the xz plane of CRV-T + Int_t dataSource = -1; //temporary variable; will be removed (1...data from stepPointMCs, 2...data from trajectory extrapolation) + CrvPlaneInfoMC(){} + CrvPlaneInfoMC(int pdgId, int primaryPdgId, float primaryE, CLHEP::Hep3Vector primaryPos, + CLHEP::Hep3Vector ppos, CLHEP::Hep3Vector pdir, float time, float kineticEnergy, int dataSource) : + pdgId(pdgId), + primaryPdgId(primaryPdgId), + primaryE(primaryE), + primary(primaryPos), + pos(ppos), + dir(pdir), + time(time), + kineticEnergy(kineticEnergy), + dataSource(dataSource) + {} + }; + + typedef std::vector CrvPlaneInfoMCCollection; //this is the MC vector which will be stored in the main TTree + +} +#endif + + diff --git a/inc/CrvPulseInfoReco.hh b/inc/CrvPulseInfoReco.hh new file mode 100644 index 0000000..7b73a24 --- /dev/null +++ b/inc/CrvPulseInfoReco.hh @@ -0,0 +1,44 @@ +#ifndef CrvPulseInfoReco_hh +#define CrvPulseInfoReco_hh + +#include "CLHEP/Vector/ThreeVector.h" +#include +#include +#include "Rtypes.h" + +namespace mu2e +{ + + struct CrvPulseInfoReco //information about CRV reco pulses + { + + XYZVectorF pos; //average position of counters + Int_t barId = -1; //CRV counter ID + Int_t sectorId = -1; //CRV sector ID + Int_t SiPMId = -1; //SiPMId number + Int_t PEs = -1; //PEs using pulse integral + Int_t PEsPulseHeight = -1; //PEs using pulse height + Float_t pulseHeight = -1; //Pulse height + Float_t pulseBeta = -1; //Pulse beta + Float_t pulseFitChi2 = -1; //Pulse Fit chi2 + Float_t time = -1; //Time + + CrvPulseInfoReco(){} + CrvPulseInfoReco(CLHEP::Hep3Vector ppos, int barId, int sectorId, int SiPMId, int PEs, int PEsPulseHeight, float pulseHeight, float pulseBeta, float pulseFitChi2, float time) : + pos(ppos), + barId(barId), + sectorId(sectorId), + SiPMId(SiPMId), + PEs(PEs), + PEsPulseHeight(PEsPulseHeight), + pulseHeight(pulseHeight), + pulseBeta(pulseBeta), + pulseFitChi2(pulseFitChi2), + time(time) + {} + }; + + typedef std::vector CrvPulseInfoRecoCollection; //this is the reco vector which will be stored in the main TTree + +} +#endif diff --git a/inc/CrvSummaryMC.hh b/inc/CrvSummaryMC.hh new file mode 100644 index 0000000..0eaa330 --- /dev/null +++ b/inc/CrvSummaryMC.hh @@ -0,0 +1,32 @@ +#ifndef CrvSummaryMC_hh +#define CrvSummaryMC_hh + +#include "Rtypes.h" +#include "Offline/DataProducts/inc/GenVector.hh" +namespace mu2e +{ + struct CrvSummaryMC + { + Double_t totalEnergyDeposited = -1; + Double_t minPathLayer = -1; + Double_t maxPathLayer = -1; + Int_t nHitCounters = -1; + XYZVectorF pos; //position of the first MC particle in CRV + Int_t sectorNumber = -1; //CRV sector number of the first MC particle in CR + Int_t sectorType = -1; //CRV sector type of the first MC particle in CRV + Int_t pdgId = -1; //PDG ID of the first MC particle in CRV + + CrvSummaryMC(){} + CrvSummaryMC(double totalEnergyDeposited, double minPathLayer, double maxPathLayer, int nHitCounters, CLHEP::Hep3Vector ppos, int sectorNumber, int sectorType, int pdgId) : + totalEnergyDeposited(totalEnergyDeposited), + minPathLayer(minPathLayer), + maxPathLayer(maxPathLayer), + nHitCounters(nHitCounters), + pos(ppos), + sectorNumber(sectorNumber), + sectorType(sectorType), + pdgId(pdgId) + {} + }; +} +#endif diff --git a/inc/CrvSummaryReco.hh b/inc/CrvSummaryReco.hh new file mode 100644 index 0000000..73ebed8 --- /dev/null +++ b/inc/CrvSummaryReco.hh @@ -0,0 +1,21 @@ +#ifndef CrvSummaryReco_hh +#define CrvSummaryReco_hh + +#include "Rtypes.h" + +namespace mu2e +{ + struct CrvSummaryReco + { + Int_t totalPEs = -1; + Int_t nHitCounters = -1; + CrvSummaryReco(){} + CrvSummaryReco(int totalPEs, int nHitCounters) : + totalPEs(totalPEs), + nHitCounters(nHitCounters) + {} + }; +} +#endif + + diff --git a/inc/CrvWaveformInfo.hh b/inc/CrvWaveformInfo.hh new file mode 100644 index 0000000..24573de --- /dev/null +++ b/inc/CrvWaveformInfo.hh @@ -0,0 +1,25 @@ +#ifndef CrvWaveformInfo_hh +#define CrvWaveformInfo_hh + +#include +#include "Rtypes.h" + +namespace mu2e +{ + struct CrvWaveformInfo //information about CRV waveforms + { + Float_t adc = -1; + Float_t time = -1; + Int_t SiPMId = -1; + CrvWaveformInfo(){} + CrvWaveformInfo(float adc, float time, int SiPMId) : + adc(adc), + time(time), + SiPMId(SiPMId) + {} + }; + + typedef std::vector CrvWaveformInfoCollection; //this is the reco vector which will be stored in the main TTree + +} +#endif diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index dc225ad..e0afdc0 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -16,11 +16,11 @@ #include "TrkAna/inc/TrkStrawHitInfoMC.hh" #include "TrkAna/inc/CaloClusterInfoMC.hh" #include "TrkAna/inc/MCStepInfo.hh" +#include "TrkAna/inc/CrvHitInfoMC.hh" #include "Offline/RecoDataProducts/inc/KalSeed.hh" #include "Offline/MCDataProducts/inc/KalSeedMC.hh" #include "BTrk/BbrGeom/HepPoint.h" #include "Offline/MCDataProducts/inc/PrimaryParticle.hh" -#include "Offline/CRVAnalysis/inc/CrvHitInfoMC.hh" #include #include @@ -64,7 +64,6 @@ namespace mu2e { void fillVDInfo(KalSeed const& kseed, const KalSeedMC& kseedmc, std::vector& vdinfos); void fillHitInfoMCs(const KalSeedMC& kseedmc, std::vector& tshinfomcs); void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, CaloClusterInfoMC& ccimc); - void fillCrvHitInfoMC(art::Ptr const& crvCoincMC, CrvHitInfoMC& crvHitInfoMC); void fillExtraMCStepInfos(KalSeedMC const& kseedmc, StepPointMCCollection const& mcsteps, MCStepInfos& mcsic, MCStepSummaryInfo& mcssi); }; diff --git a/inc/InfoStructHelper.hh b/inc/InfoStructHelper.hh index c0e544d..733df7f 100644 --- a/inc/InfoStructHelper.hh +++ b/inc/InfoStructHelper.hh @@ -25,7 +25,7 @@ #include "TrkAna/inc/TrkQualInfo.hh" #include "TrkAna/inc/TrkPIDInfo.hh" #include "TrkAna/inc/HelixInfo.hh" -#include "Offline/CRVAnalysis/inc/CrvHitInfoReco.hh" +#include "TrkAna/inc/CrvHitInfoReco.hh" #include #include @@ -62,7 +62,6 @@ namespace mu2e { void fillTrkQualInfo(const TrkQual& tqual, TrkQualInfo& trkqualInfo); void fillTrkPIDInfo(const TrkCaloHitPID& tchp, const KalSeed& kseed, TrkPIDInfo& trkpidInfo); void fillHelixInfo(art::Ptr const& hptr, HelixInfo& hinfo); - void fillCrvHitInfo(art::Ptr const& crvCoinc, CrvHitInfoReco& crvHitInfo); }; } diff --git a/inc/TrkInfo.hh b/inc/TrkInfo.hh index 8030840..fe989e5 100644 --- a/inc/TrkInfo.hh +++ b/inc/TrkInfo.hh @@ -69,8 +69,6 @@ namespace mu2e int nnullambig =0; // number of hits without any ambiguity assigned int nmat =0, nmatactive =0; // number materials (straw) assigned and used (active) to this fit int nseg =0; // number of trajectory segments - float t0 =0; // time this particle was estimated to cross z=0 - float t0err =-1; // error on t0 float chisq =-1; // Kalman fit chisquared float fitcon =-1; // Kalman fit chisqured consistency float radlen; // total radiation length of (active) material crossed by this particle inside the tracker diff --git a/src/.InfoMCStructHelper.cc.swp b/src/.InfoMCStructHelper.cc.swp new file mode 100644 index 0000000..08cde10 Binary files /dev/null and b/src/.InfoMCStructHelper.cc.swp differ diff --git a/src/CrvInfoHelper.cc b/src/CrvInfoHelper.cc new file mode 100644 index 0000000..27b375c --- /dev/null +++ b/src/CrvInfoHelper.cc @@ -0,0 +1,263 @@ +#include "TrkAna/inc/CrvInfoHelper.hh" + +#include "Offline/CosmicRayShieldGeom/inc/CosmicRayShield.hh" +#include "Offline/CRVConditions/inc/CRVDigitizationPeriod.hh" +#include "Offline/GeometryService/inc/DetectorSystem.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/MCDataProducts/inc/CrvCoincidenceClusterMC.hh" +#include "Offline/MCDataProducts/inc/CrvDigiMC.hh" +#include "Offline/MCDataProducts/inc/CrvStep.hh" +#include "Offline/MCDataProducts/inc/MCTrajectory.hh" +#include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" +#include "Offline/RecoDataProducts/inc/CrvDigi.hh" +#include "Offline/RecoDataProducts/inc/CrvRecoPulse.hh" +#include "Offline/DataProducts/inc/PDGCode.hh" +#include "art/Framework/Principal/Handle.h" +#include "Offline/CRVResponse/inc/CrvMCHelper.hh" +#include "Offline/CRVReco/inc/CrvHelper.hh" +#include "Offline/GeometryService/inc/DetectorSystem.hh" + +namespace mu2e +{ + void CrvInfoHelper::FillCrvHitInfoCollections( + art::Handle const& crvCoincidences, + art::Handle const& crvCoincidencesMC, + art::Handle const& crvRecoPulses, + art::Handle const& crvSteps, + art::Handle const& mcTrajectories, + CrvHitInfoRecoCollection &recoInfo, CrvHitInfoMCCollection &MCInfo, + CrvSummaryReco &recoSummary, CrvSummaryMC &MCSummary, + CrvPlaneInfoMCCollection &MCInfoPlane, double crvPlaneY) { + GeomHandle CRS; + GeomHandle tdet; + if(!crvCoincidences.isValid()) return; + size_t nClusters=crvCoincidences->size(); + for(size_t i=0; iat(i); + + //fill the Reco collection + recoInfo.emplace_back( + cluster.GetCrvSectorType(), + tdet->toDetector(cluster.GetAvgHitPos()), + cluster.GetStartTime(), cluster.GetEndTime(), cluster.GetAvgHitTime(), + cluster.GetPEs(), + cluster.GetCrvRecoPulses().size(), + cluster.GetLayers().size(), + cluster.GetSlope()); + } + + if(!crvRecoPulses.isValid()) return; + size_t nRecoPulses=crvRecoPulses->size(); + recoSummary.totalPEs=0; + std::set counters; + for(size_t i=0; iat(i).GetPEs(); + counters.insert(crvRecoPulses->at(i).GetScintillatorBarIndex()); + } + recoSummary.nHitCounters=counters.size(); + + + //fill the MC collection + if(crvCoincidencesMC.isValid()) + { + size_t nClustersMC=crvCoincidencesMC->size(); + if(nClusters!=nClustersMC) std::cout<<"The number of MC and reco CRV coincidence clusters does not match!"<at(i); + if(clusterMC.HasMCInfo()) + { + const art::Ptr &simParticle = clusterMC.GetMostLikelySimParticle(); + const art::Ptr &primaryParticle = FindPrimaryParticle(simParticle); + const art::Ptr &parentParticle = FindParentParticle(simParticle); + const art::Ptr &gparentParticle = FindGParentParticle(simParticle); + MCInfo.emplace_back( + true, + simParticle->pdgId(), + primaryParticle->pdgId(), + primaryParticle->startMomentum().e() - primaryParticle->startMomentum().m(), + tdet->toDetector(primaryParticle->startPosition()), + parentParticle->pdgId(), + parentParticle->startMomentum().e() - parentParticle->startMomentum().m(), + tdet->toDetector(parentParticle->startPosition()), + gparentParticle->pdgId(), + gparentParticle->startMomentum().e() - gparentParticle->startMomentum().m(), + tdet->toDetector(gparentParticle->startPosition()), + tdet->toDetector(clusterMC.GetAvgHitPos()), + clusterMC.GetAvgHitTime(), + clusterMC.GetTotalEnergyDeposited()); + } + else MCInfo.emplace_back(); + } + } + + MCSummary = CrvSummaryMC(); + if(crvSteps.isValid()) + { + size_t nSteps=crvSteps->size(); + MCSummary.totalEnergyDeposited=0; + std::set counters; + double totalStep[] = {0, 0, 0, 0}; + for(size_t i=0; iat(i).visibleEDep(); + counters.insert(crvSteps->at(i).barIndex()); + const CRSScintillatorBarId &CRVCounterId = CRS->getBar(crvSteps->at(i).barIndex()).id(); + int layer = CRVCounterId.getLayerNumber(); + int pdgId = crvSteps->at(i).simParticle()->pdgId(); + if(abs(pdgId)==PDGCode::mu_minus) + totalStep[layer] = totalStep[layer] + crvSteps->at(i).pathLength(); + + // Save info from the first step in the CRV + if(i==0){ + CLHEP::Hep3Vector CrvPos = crvSteps->at(i).startPosition(); + MCSummary.pos = XYZVectorF(tdet->toDetector(CrvPos)); + int sectorNumber = CRVCounterId.getShieldNumber(); + MCSummary.sectorNumber = sectorNumber; + MCSummary.sectorType = CRS->getCRSScintillatorShield(sectorNumber).getSectorType(); + MCSummary.pdgId = pdgId; + } + } + + MCSummary.nHitCounters=counters.size(); + MCSummary.minPathLayer=*std::min_element(totalStep,totalStep+4); + MCSummary.maxPathLayer=*std::max_element(totalStep,totalStep+4); + } + + //locate points where the cosmic MC trajectories cross the xz plane of CRV-T + if(mcTrajectories.isValid()) + { + for(auto trajectoryIter=mcTrajectories->begin(); trajectoryIter!=mcTrajectories->end(); trajectoryIter++) + { + const art::Ptr &trajectorySimParticle = trajectoryIter->first; + if(abs(trajectorySimParticle->pdgId())!=PDGCode::mu_minus) continue; + const art::Ptr &trajectoryPrimaryParticle = FindPrimaryParticle(trajectorySimParticle); + GenId genId = trajectoryPrimaryParticle->genParticle()->generatorId(); + if(genId==GenId::cosmicToy || genId==GenId::cosmicDYB || genId==GenId::cosmic || genId==GenId::cosmicCRY) + { + const std::vector &points = trajectoryIter->second.points(); + if(points.size()<1) continue; + CLHEP::Hep3Vector previousPos=points[0].pos(); + for(size_t i=1; icrvPlaneY && pos.y()<=crvPlaneY) || (previousPos.y()=crvPlaneY)) + { + double fraction=(crvPlaneY-pos.y())/(previousPos.y()-pos.y()); + CLHEP::Hep3Vector planePos=fraction*(previousPos-pos)+pos; + CLHEP::Hep3Vector planeDir=(pos-previousPos).unit(); + double planeTime=fraction*(points[i-1].t()-points[i].t())+points[i].t(); + double planeKineticEnergy=fraction*(points[i-1].kineticEnergy()-points[i].kineticEnergy())+points[i].kineticEnergy(); + MCInfoPlane.emplace_back(trajectorySimParticle->pdgId(), + trajectoryPrimaryParticle->pdgId(), + trajectoryPrimaryParticle->startMomentum().e(), + tdet->toDetector(trajectoryPrimaryParticle->startPosition()), + tdet->toDetector(planePos), + planeDir, + planeTime, + planeKineticEnergy, + 0); //unused + } + previousPos=pos; + } + } + } + } + + }//FillCrvInfoStructure + + void CrvInfoHelper::FillCrvPulseInfoCollections ( + art::Handle const& crvRecoPulses, + art::Handle const& crvDigiMCs, + art::Handle const& crvDigis, + CrvPulseInfoRecoCollection &recoInfo, CrvHitInfoMCCollection &MCInfo, CrvWaveformInfoCollection &waveformInfo){ + GeomHandle tdet; + + if(!crvRecoPulses.isValid()) return; + + GeomHandle CRS; + + // Create SiPM map to extract sequantial SiPM IDs + const std::vector > &counters = CRS->getAllCRSScintillatorBars(); + std::vector >::const_iterator iter; + int iSiPM = 0; + std::map sipm_map; + for(iter=counters.begin(); iter!=counters.end(); iter++) + { + const CRSScintillatorBarIndex &barIndex = (*iter)->index(); + for(int SiPM=0; SiPM<4; SiPM++) + { + if(!(*iter)->getBarDetail().hasCMB(SiPM%2)) continue; + sipm_map[barIndex.asInt()*4 + SiPM] = iSiPM; + iSiPM++; + } + } + + // Loop through all reco pulses + for(size_t recoPulseIndex=0; recoPulseIndexsize(); recoPulseIndex++) + { + const art::Ptr crvRecoPulse(crvRecoPulses, recoPulseIndex); + //get information about the counter + const CRSScintillatorBarIndex &barIndex = crvRecoPulse->GetScintillatorBarIndex(); + int sectorNumber = -1; + int moduleNumber = -1; + int layerNumber = -1; + int counterNumber = -1; + CrvHelper::GetCrvCounterInfo(CRS, barIndex, sectorNumber, moduleNumber, layerNumber, counterNumber); + + //Reco pulses information + int SiPM = crvRecoPulse->GetSiPMNumber(); + int SiPMId = sipm_map.find(barIndex.asInt()*4 + SiPM)->second; + CLHEP::Hep3Vector HitPos = CrvHelper::GetCrvCounterPos(CRS, barIndex); + recoInfo.emplace_back(HitPos, barIndex.asInt(), sectorNumber, SiPMId, + crvRecoPulse->GetPEs(), crvRecoPulse->GetPEsPulseHeight(), crvRecoPulse->GetPulseHeight(), + crvRecoPulse->GetPulseBeta(), crvRecoPulse->GetPulseFitChi2(), crvRecoPulse->GetPulseTime()); + + //MCtruth pulses information + double visibleEnergyDeposited = 0; + double earliestHitTime = 0; + double avgHitTime = 0; + CLHEP::Hep3Vector earliestHitPos; + CLHEP::Hep3Vector avgHitPos; + art::Ptr mostLikelySimParticle; + //for this reco pulse + CrvMCHelper::GetInfoFromCrvRecoPulse(crvRecoPulse, crvDigiMCs, visibleEnergyDeposited,earliestHitTime, earliestHitPos, avgHitTime, avgHitPos, + mostLikelySimParticle); + + bool hasMCInfo = (mostLikelySimParticle.isNonnull()?true:false); //MC + if(hasMCInfo) + { + const art::Ptr &primaryParticle = FindPrimaryParticle(mostLikelySimParticle); + const art::Ptr &parentParticle = FindParentParticle(mostLikelySimParticle); + const art::Ptr &gparentParticle = FindGParentParticle(mostLikelySimParticle); + MCInfo.emplace_back(true, mostLikelySimParticle->pdgId(), + primaryParticle->pdgId(), + primaryParticle->startMomentum().e() - primaryParticle->startMomentum().m(), + tdet->toDetector(primaryParticle->startPosition()), + parentParticle->pdgId(), + parentParticle->startMomentum().e() - parentParticle->startMomentum().m(), + tdet->toDetector(parentParticle->startPosition()), + gparentParticle->pdgId(), + gparentParticle->startMomentum().e() - gparentParticle->startMomentum().m(), + tdet->toDetector(gparentParticle->startPosition()), + tdet->toDetector(avgHitPos), + avgHitTime, + visibleEnergyDeposited); + } + else + MCInfo.emplace_back(); + } + + // Fill waveforms struct + for(size_t j=0; jsize(); j++) + { + mu2e::CrvDigi const& digi(crvDigis->at(j)); + int SiPMId = sipm_map.find(digi.GetScintillatorBarIndex().asInt()*4 + digi.GetSiPMNumber())->second; + for(size_t k=0; k const& crvCoincMC, CrvHitInfoMC& crvHitInfoMC) { - const art::Ptr &simParticle = crvCoincMC->GetMostLikelySimParticle(); - if(crvCoincMC->HasMCInfo() && simParticle.isAvailable()) { - art::Ptr primaryParticle = simParticle; - art::Ptr parentParticle = simParticle; - art::Ptr gparentParticle = simParticle; - int i_gen = 0; - while (primaryParticle->hasParent()) { - primaryParticle = primaryParticle->parent(); - if (i_gen ==0) { - parentParticle = primaryParticle; - gparentParticle = primaryParticle; - } - else if (i_gen == 1) { - gparentParticle = primaryParticle; - } - ++i_gen; - } - crvHitInfoMC._valid = true; - crvHitInfoMC._pdgId = simParticle->pdgId(); - crvHitInfoMC._primaryPdgId = primaryParticle->pdgId(); - crvHitInfoMC._primaryE = primaryParticle->startMomentum().e() - primaryParticle->startMomentum().m(); - crvHitInfoMC._primaryX = primaryParticle->startPosition().x(); - crvHitInfoMC._primaryY = primaryParticle->startPosition().y(); - crvHitInfoMC._primaryZ = primaryParticle->startPosition().z(); - crvHitInfoMC._parentPdgId = parentParticle->pdgId(); - crvHitInfoMC._parentE = parentParticle->startMomentum().e() - parentParticle->startMomentum().m(); - crvHitInfoMC._parentX = parentParticle->startPosition().x(); - crvHitInfoMC._parentY = parentParticle->startPosition().y(); - crvHitInfoMC._parentZ = parentParticle->startPosition().z(); - crvHitInfoMC._gparentPdgId = gparentParticle->pdgId(); - crvHitInfoMC._gparentE = gparentParticle->startMomentum().e() - gparentParticle->startMomentum().m(); - crvHitInfoMC._gparentX = gparentParticle->startPosition().x(); - crvHitInfoMC._gparentY = gparentParticle->startPosition().y(); - crvHitInfoMC._gparentZ = gparentParticle->startPosition().z(); - crvHitInfoMC._x = crvCoincMC->GetEarliestHitPos().x(); - crvHitInfoMC._y = crvCoincMC->GetEarliestHitPos().y(); - crvHitInfoMC._z = crvCoincMC->GetEarliestHitPos().z(); - crvHitInfoMC._time = crvCoincMC->GetEarliestHitTime(); - crvHitInfoMC._depositedEnergy = crvCoincMC->GetTotalEnergyDeposited(); - } - } - void InfoMCStructHelper::fillExtraMCStepInfos(KalSeedMC const& kseedmc, StepPointMCCollection const& mcsteps, MCStepInfos& mcsic, MCStepSummaryInfo& mcssi) { GeomHandle det; diff --git a/src/InfoStructHelper.cc b/src/InfoStructHelper.cc index 3bbb8e3..bd2541e 100644 --- a/src/InfoStructHelper.cc +++ b/src/InfoStructHelper.cc @@ -61,8 +61,6 @@ namespace mu2e { trkinfo.fitalg = 0; trkinfo.pdg = kseed.particle(); - trkinfo.t0 = kseed.t0().t0(); - trkinfo.t0err = kseed.t0().t0Err(); fillTrkInfoHits(kseed, trkinfo); @@ -420,14 +418,4 @@ namespace mu2e { } } - void InfoStructHelper::fillCrvHitInfo(art::Ptr const& crvCoinc, CrvHitInfoReco& crvHitInfo) { - crvHitInfo._crvSectorType = crvCoinc->GetCrvSectorType(); - crvHitInfo._x = crvCoinc->GetAvgHitPos().x(); - crvHitInfo._y = crvCoinc->GetAvgHitPos().y(); - crvHitInfo._z = crvCoinc->GetAvgHitPos().z(); - crvHitInfo._timeWindowStart = crvCoinc->GetStartTime(); - crvHitInfo._timeWindowEnd = crvCoinc->GetEndTime(); - crvHitInfo._PEs = crvCoinc->GetPEs(); - crvHitInfo._nCoincidenceHits = crvCoinc->GetCrvRecoPulses().size(); - } } diff --git a/src/SConscript b/src/SConscript index 2999912..4e683c2 100644 --- a/src/SConscript +++ b/src/SConscript @@ -23,6 +23,7 @@ mainlib = helper.make_mainlib ( [ 'art_Framework_Services_Registry', 'art_Utilities', 'mu2e_GlobalConstantsService', + 'mu2e_GeometryService', 'mu2e_KinKalGeom', 'mu2e_Mu2eInterfaces', 'mu2e_Mu2eUtilities', @@ -32,10 +33,14 @@ mainlib = helper.make_mainlib ( [ 'mu2e_MCDataProducts', 'mu2e_RecoDataProducts', 'mu2e_CalorimeterGeom', + 'mu2e_CosmicRayShieldGeom', + 'mu2e_CRVResponse', + 'mu2e_CRVReco', + 'mu2e_CRVReco', 'KinKal_Geometry', 'KinKal_Trajectory', 'KinKal_General' -] ) + ] ) # Fixme: split into link lists for each module. helper.make_plugins( [ @@ -60,8 +65,8 @@ helper.make_plugins( [ 'mu2e_RecoDataProducts', 'mu2e_Mu2eUtilities', 'mu2e_GeneralUtilities', + 'mu2e_CRVReco', 'mu2e_BFieldGeom', - 'mu2e_CRVAnalysis', 'mu2e_TrkDiag', 'mu2e_KinKalGeom', 'KinKal_Geometry', diff --git a/src/TrkAnaTreeMaker_module.cc b/src/TrkAnaTreeMaker_module.cc index 0a1131a..58a5df3 100644 --- a/src/TrkAnaTreeMaker_module.cc +++ b/src/TrkAnaTreeMaker_module.cc @@ -29,6 +29,7 @@ #include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" #include "Offline/MCDataProducts/inc/CrvCoincidenceClusterMC.hh" #include "Offline/MCDataProducts/inc/CrvCoincidenceClusterMCAssns.hh" +#include "Offline/MCDataProducts/inc/CrvDigiMC.hh" #include "Offline/Mu2eUtilities/inc/fromStrings.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" @@ -77,12 +78,12 @@ #include "TrkAna/inc/TrkPIDInfo.hh" #include "TrkAna/inc/HelixInfo.hh" #include "TrkAna/inc/InfoStructHelper.hh" +#include "TrkAna/inc/CrvInfoHelper.hh" #include "TrkAna/inc/InfoMCStructHelper.hh" #include "Offline/RecoDataProducts/inc/RecoQual.hh" #include "TrkAna/inc/RecoQualInfo.hh" #include "TrkAna/inc/BestCrvAssns.hh" #include "TrkAna/inc/MCStepInfo.hh" -#include "Offline/CRVAnalysis/inc/CRVAnalysis.hh" // C++ includes. #include @@ -112,10 +113,6 @@ namespace mu2e { fhicl::Atom filltrkpid{Name("fillTrkPID"), Comment("Switch to turn on filling of the full TrkPIDInfo for this set of tracks"), false}; fhicl::Atom required{Name("required"), Comment("True/false if you require this type of track in the event"), false}; fhicl::Atom genealogyDepth{Name("genealogyDepth"), Comment("The depth of the genealogy information you want to keep"), 1}; - fhicl::OptionalSequence bestCrvModules{Name("bestCrvModules"), Comment("Sequence of module labels that create the BestCrvAssns you want written out (use prefix if fcl parameter suffix (e.g. DeM) is defined)")}; - fhicl::OptionalSequence bestCrvInstances{Name("bestCrvInstances"), Comment("Sequence of instance names for modules that create multiple BestCrvAssns")}; - fhicl::OptionalSequence bestCrvBranches{Name("bestCrvBranches"), Comment("Sequence of branch names that will be created to store the bestcrv information")}; - fhicl::Atom fillbestcrv{Name("fillBestCrv"), Comment("Switch to turn on filling of the bestcrv branch for this set of tracks"), false}; }; struct BranchConfig { @@ -153,31 +150,31 @@ namespace mu2e { fhicl::Atom fittype{Name("FitType"),Comment("Type of track Fit: LoopHelix, CentralHelix, KinematicLine, or Unknown"),"Unknown"}; fhicl::Atom helices{Name("FillHelixInfo"),false}; // CRV -- input tags - fhicl::Atom crvCoincidenceModuleLabel{Name("CrvCoincidenceModuleLabel"), Comment("CrvCoincidenceModuleLabel")}; - fhicl::Atom crvRecoPulseLabel{Name("CrvRecoPulseLabel"), Comment("CrvRecoPulseLabel")}; - fhicl::Atom crvStepLabel{Name("CrvStepLabel"), Comment("CrvStepLabel")}; - fhicl::Atom crvWaveformsModuleLabel{ Name("CrvWaveformsModuleLabel"), Comment("CrvWaveformsModuleLabel")}; - fhicl::Atom crvDigiModuleLabel{ Name("CrvDigiModuleLabel"), Comment("CrvDigiModuleLabel")}; + fhicl::Atom crvCoincidencesTag{Name("CrvCoincidencesTag"), Comment("Tag for CrvCoincidenceCluster Collection"), art::InputTag()}; + fhicl::Atom crvRecoPulsesTag{Name("CrvRecoPulsesTag"), Comment("Tag for CrvRecopPulse Collection"), art::InputTag()}; + fhicl::Atom crvStepsTag{Name("CrvStepsTag"), Comment("Tag for CrvStep Collection"), art::InputTag()}; + fhicl::Atom crvDigiMCsTag{Name("CrvDigiMCsTag"), Comment("Tag for CrvDigiMC Collection"), art::InputTag()}; + fhicl::Atom crvDigisTag{Name("CrvDigisTag"), Comment("Tag for CrvDigi Collection"), art::InputTag()}; // CRV -- flags - fhicl::Atom crv{Name("FillCRV"),Comment("Flag for turning on bestcrv(mc) branches"), false}; - fhicl::Atom crvhits{Name("FillCRVHits"), Comment("Flag for turning on crvinfo(mc), crvsummary(mc), and crvinfomcplane branches"), false}; + fhicl::Atom crvhits{Name("FillCRVHits"),Comment("Flag for turning on crv CoincidenceClusterbranches"), false}; fhicl::Atom crvpulses{Name("FillCRVPulses"),Comment("Flag for turning on crvpulseinfo(mc), crvwaveforminfo branches"), false}; // CRV -- other - fhicl::Atom crvPlaneY{Name("CrvPlaneY"),2751.485}; //y of center of the top layer of the CRV-T counters + fhicl::Atom crvPlaneY{Name("CrvPlaneY"),2751.485}; //y of center of the top layer of the CRV-T counters. This belongs in KinKalGeom as an intersection plane, together with the rest of the CRV planes FIXME // MC truth fhicl::Atom fillmc{Name("FillMCInfo"),Comment("Global switch to turn on/off MC info"),true}; fhicl::Table infoMCStructHelper{Name("InfoMCStructHelper"), Comment("Configuration for the InfoMCStructHelper")}; fhicl::Atom PBTMCTag{Name("PBTMCTag"), Comment("Tag for ProtonBunchTimeMC object") ,art::InputTag()}; - fhicl::Atom simParticleLabel{Name("SimParticleLabel"), Comment("SimParticleLabel")}; - fhicl::Atom mcTrajectoryLabel{Name("MCTrajectoryLabel"), Comment("MCTrajectoryLabel")}; + fhicl::Atom simParticlesTag{Name("SimParticlesTag"), Comment("SimParticle Collection Tag")}; + fhicl::Atom mcTrajectoriesTag{Name("MCTrajectoriesTag"), Comment("MCTrajectory Collection Tag")}; fhicl::Atom primaryParticleTag{Name("PrimaryParticleTag"), Comment("Tag for PrimaryParticle"), art::InputTag()}; fhicl::Atom kalSeedMCTag{Name("KalSeedMCAssns"), Comment("Tag for KalSeedMCAssn"), art::InputTag()}; + // extra MC fhicl::OptionalSequence extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out")}; // Calo MC fhicl::Atom fillCaloMC{ Name("FillCaloMC"),Comment("Fill CaloMC information"), true}; fhicl::Atom caloClusterMCTag{Name("CaloClusterMCTag"), Comment("Tag for CaloClusterMCCollection") ,art::InputTag()}; // CRV MC - fhicl::Atom crvCoincidenceMCModuleLabel{Name("CrvCoincidenceMCModuleLabel"), Comment("CrvCoincidenceMCModuleLabel")}; + fhicl::Atom crvCoincidenceMCsTag{Name("CrvCoincidenceMCsTag"), Comment("Tag for CrvCoincidenceClusterMC Collection"), art::InputTag()}; fhicl::Atom crvMCAssnsTag{ Name("CrvCoincidenceClusterMCAssnsTag"), Comment("art::InputTag for CrvCoincidenceClusterMCAssns")}; // Pre-processed analysis info; are these redundant with the branch config ? fhicl::Atom filltrkqual{Name("FillTrkQualInfo"),false}; @@ -204,10 +201,6 @@ namespace mu2e { EventInfo _einfo; EventInfoMC _einfomc; art::InputTag _PBITag, _PBTTag, _PBTMCTag; - std::vector _extraMCStepTags; - std::vector> _extraMCStepCollections; - std::vector> _extraMCStepInfos; - std::vector> _extraMCStepSummaryInfos; // hit counting HitCount _hcnt; // track counting @@ -236,9 +229,15 @@ namespace mu2e { // MC truth (fcl parameters) bool _fillmc; // MC truth inputs + std::vector _extraMCStepTags; + std::vector> _extraMCStepCollections; + std::vector> _extraMCStepInfos; + std::vector> _extraMCStepSummaryInfos; + // art::Handle _pph; art::Handle _ksmcah; - art::InputTag _primaryParticleTag; + art::Handle _simParticles; + art::Handle _mcTrajectories; // MC truth branches (outputs) std::vector _allMCTIs; std::map> _allMCSimTIs; @@ -255,39 +254,36 @@ namespace mu2e { // event weights std::vector > _wtHandles; EventWeightInfo _wtinfo; - // CRV - // CRV -- fhicl parameters - bool _crv; - bool _crvhits; - bool _crvpulses; - std::string _crvCoincidenceModuleLabel; - std::string _crvCoincidenceMCModuleLabel; - std::string _crvRecoPulseLabel; - std::string _crvStepLabel; - std::string _crvWaveformsModuleLabel; - std::string _crvDigiModuleLabel; - art::InputTag _crvMCAssnsTag; - double _crvPlaneY; // CRV (inputs) std::map>> _allBestCrvAssns; - art::Handle _crvCoincidenceMCCollHandle; - art::Handle _crvMCAssnsHandle; + art::Handle _crvMCAssns; + art::Handle _crvCoincidences; + art::Handle _crvCoincidenceMCs; + art::Handle _crvRecoPulses; + art::Handle _crvDigiMCs; + art::Handle _crvDigis; + art::Handle _crvSteps; + // CRV -- fhicl parameters + bool _crvhits, _crvpulses; + double _crvPlaneY; // needs to move to KinKalGeom FIXME // CRV (output) - std::vector _crvinfo; + std::vector _crvhit; std::map> _allBestCrvs; // there can be more than one of these per candidate/supplement - std::vector _crvinfomc; + std::vector _crvhitmc; std::map> _allBestCrvMCs; CrvSummaryReco _crvsummary; CrvSummaryMC _crvsummarymc; - std::vector _crvinfomcplane; + std::vector _crvhitmcplane; std::vector _crvpulseinfo; std::vector _crvwaveforminfo; std::vector _crvpulseinfomc; + std::vector _crvrecoinfo; // helices HelixInfo _hinfo; // struct helpers InfoStructHelper _infoStructHelper; InfoMCStructHelper _infoMCStructHelper; + CrvInfoHelper _crvHelper; // branch structure Int_t _buffsize, _splitlevel; enum FType{Unknown=0,LoopHelix,CentralHelix,KinematicLine}; @@ -316,17 +312,7 @@ namespace mu2e { _fillmc(conf().fillmc()), _fillcalomc(conf().fillCaloMC()), // CRV - _crv(conf().crv()), _crvhits(conf().crvhits()), - _crvpulses(conf().crvpulses()), - _crvCoincidenceModuleLabel(conf().crvCoincidenceModuleLabel()), - _crvCoincidenceMCModuleLabel(conf().crvCoincidenceMCModuleLabel()), - _crvRecoPulseLabel(conf().crvRecoPulseLabel()), - _crvStepLabel(conf().crvStepLabel()), - _crvWaveformsModuleLabel(conf().crvWaveformsModuleLabel()), - _crvDigiModuleLabel(conf().crvDigiModuleLabel()), - _crvMCAssnsTag(conf().crvMCAssnsTag()), - _crvPlaneY(conf().crvPlaneY()), _infoMCStructHelper(conf().infoMCStructHelper()), _buffsize(conf().buffsize()), _splitlevel(conf().splitlevel()) @@ -392,19 +378,6 @@ namespace mu2e { std::vector tshimc; _allTSHIMCs.push_back(tshimc); - if(_crv && i_branchConfig.options().fillbestcrv()) { // if we are filling in bestcrv information - std::vector bestcrv; - std::vector bestcrvmc; - std::vector bestCrvBranchNames; // need to know how many bestcrv branches there will be per candidate/supplement - if (i_branchConfig.options().bestCrvBranches(bestCrvBranchNames)) { - for (size_t i_bestCrvBranch = 0; i_bestCrvBranch < bestCrvBranchNames.size(); ++i_bestCrvBranch) { - bestcrv.push_back(CrvHitInfoReco()); // add empty structs to the vector so that ROOT can be given a location to find it - bestcrvmc.push_back(CrvHitInfoMC()); - } - _allBestCrvs[i_branch] = bestcrv; - _allBestCrvMCs[i_branch] = bestcrvmc; - } - } } } @@ -478,19 +451,6 @@ namespace mu2e { _trkana->Branch((branch+"tsh.").c_str(),&_allTSHIs.at(i_branch),_buffsize,_splitlevel); _trkana->Branch((branch+"tsm.").c_str(),&_allTSMIs.at(i_branch),_buffsize,_splitlevel); } - // want to be able to have more than one bestcrv branch for each candidate/supplement - if(_crv && i_branchConfig.options().fillbestcrv()) { - std::vector bestCrvBranchNames; - if (i_branchConfig.options().bestCrvBranches(bestCrvBranchNames)) { - for (size_t i_bestCrvBranch = 0; i_bestCrvBranch < bestCrvBranchNames.size(); ++i_bestCrvBranch) { - std::string bestCrvBranchName = bestCrvBranchNames.at(i_bestCrvBranch); - _trkana->Branch((branch+bestCrvBranchName+".").c_str(),&_allBestCrvs.at(i_branch).at(i_bestCrvBranch),_buffsize,_splitlevel); - if (_fillmc) { - _trkana->Branch((branch+bestCrvBranchName+"mc.").c_str(),&_allBestCrvMCs.at(i_branch).at(i_bestCrvBranch),_buffsize,_splitlevel); - } - } - } - } // optionally add MC branches if(_fillmc && i_branchConfig.options().fillmc()){ @@ -528,24 +488,21 @@ namespace mu2e { } // calorimeter information for the downstream electron track // general CRV info - if(_crv) { - if (_crvhits) { - _trkana->Branch("crvsummary.",&_crvsummary,_buffsize,_splitlevel); - _trkana->Branch("crvinfo.",&_crvinfo,_buffsize,_splitlevel); - if(_crvpulses) { - _trkana->Branch("crvpulseinfo.",&_crvpulseinfo,_buffsize,_splitlevel); - _trkana->Branch("crvwaveforminfo.",&_crvwaveforminfo,_buffsize,_splitlevel); - } + if(_crvhits) { + // coincidence branches should be here FIXME + _trkana->Branch("crvsummary.",&_crvsummary,_buffsize,_splitlevel); + _trkana->Branch("crvhit.",&_crvhit,_buffsize,_splitlevel); + if(_crvpulses) { + _trkana->Branch("crvpulseinfo.",&_crvpulseinfo,_buffsize,_splitlevel); + _trkana->Branch("crvwaveforminfo.",&_crvwaveforminfo,_buffsize,_splitlevel); } if(_fillmc){ - if (_crvhits) { - _trkana->Branch("crvsummarymc.",&_crvsummarymc,_buffsize,_splitlevel); - _trkana->Branch("crvinfomc.",&_crvinfomc,_buffsize,_splitlevel); - _trkana->Branch("crvinfomcplane.",&_crvinfomcplane,_buffsize,_splitlevel); - if(_crvpulses) { - _trkana->Branch("crvpulseinfomc.",&_crvpulseinfomc,_buffsize,_splitlevel); - } + _trkana->Branch("crvsummarymc.",&_crvsummarymc,_buffsize,_splitlevel); + _trkana->Branch("crvhitmc.",&_crvhitmc,_buffsize,_splitlevel); + _trkana->Branch("crvhitmcplane.",&_crvhitmcplane,_buffsize,_splitlevel); + if(_crvpulses) { + _trkana->Branch("crvpulseinfomc.",&_crvpulseinfomc,_buffsize,_splitlevel); } } } @@ -633,26 +590,6 @@ namespace mu2e { } } _allTCHPCHs.push_back(trkpidCollHandle); - - // BestCrv - if (i_branchConfig.options().fillbestcrv() && _crv) { // if we are filling in bestcrv information - std::vector i_bestcrv_tags; - std::vector i_bestcrv_instances; - std::vector> bestCrvAssnsHandles; - if (i_branchConfig.options().bestCrvModules(i_bestcrv_tags) && i_branchConfig.options().bestCrvInstances(i_bestcrv_instances)) { // get the module labels and instances names - // loop htrough the module lables - art::Handle bestCrvAssnsHandle; - for (size_t i_bestCrvBranch = 0; i_bestCrvBranch < i_bestcrv_tags.size(); ++i_bestCrvBranch) { - art::InputTag bestCrvInputTag = i_bestcrv_tags.at(i_bestCrvBranch) + i_branchConfig.suffix() + ":" + i_bestcrv_instances.at(i_bestCrvBranch); - event.getByLabel(bestCrvInputTag,bestCrvAssnsHandle); // get the Assns for this module and this candidate/supplement - bestCrvAssnsHandles.push_back(bestCrvAssnsHandle); // add this handle to the vector of handles for this candidate/supplement - } - } - _allBestCrvAssns[i_branch] = bestCrvAssnsHandles; - if (_fillmc) { - event.getByLabel(_crvMCAssnsTag, _crvMCAssnsHandle); - } - } } // trigger information @@ -663,6 +600,8 @@ namespace mu2e { if(_fillmc) { // get MC product collections event.getByLabel(_conf.primaryParticleTag(),_pph); event.getByLabel(_conf.kalSeedMCTag(),_ksmcah); + event.getByLabel(_conf.simParticlesTag(),_simParticles); + event.getByLabel(_conf.mcTrajectoriesTag(),_mcTrajectories); if(_fillcalomc)event.getByLabel(_conf.caloClusterMCTag(),_ccmcch); } // fill track counts @@ -696,7 +635,6 @@ namespace mu2e { _infoStructHelper.fillHelixInfo(hptr, _hinfo); } - // Now loop through all the branches (both candidate + supplements)... for (BranchIndex i_branch = 0; i_branch < _allBranches.size(); ++i_branch) { if (i_branch == _candidateIndex) { // ...but actually ignore candidate @@ -729,16 +667,22 @@ namespace mu2e { // TODO we want MC information when we don't have a track // fill general CRV info - if(_crv){ - if (_crvhits) { - CRVAnalysis::FillCrvHitInfoCollections(_crvCoincidenceModuleLabel, _crvCoincidenceMCModuleLabel, - _crvRecoPulseLabel, _crvStepLabel, _conf.simParticleLabel(), - _conf.mcTrajectoryLabel(), event, _crvinfo, _crvinfomc, - _crvsummary, _crvsummarymc, _crvinfomcplane, _crvPlaneY); + if(_crvhits){ + event.getByLabel(_conf.crvCoincidencesTag(),_crvCoincidences); + event.getByLabel(_conf.crvRecoPulsesTag(),_crvRecoPulses); + event.getByLabel(_conf.crvStepsTag(),_crvSteps); + event.getByLabel(_conf.crvDigisTag(),_crvDigis); + if(_fillmc){ + event.getByLabel(_conf.crvCoincidenceMCsTag(),_crvCoincidenceMCs); + event.getByLabel(_conf.crvDigiMCsTag(),_crvDigiMCs); } + _crvHelper.FillCrvHitInfoCollections( + _crvCoincidences, _crvCoincidenceMCs, + _crvRecoPulses, _crvSteps, _mcTrajectories,_crvhit, _crvhitmc, + _crvsummary, _crvsummarymc, _crvhitmcplane, _crvPlaneY); if(_crvpulses){ - CRVAnalysis::FillCrvPulseInfoCollections(_crvRecoPulseLabel, _crvWaveformsModuleLabel, _crvDigiModuleLabel, - event, _crvpulseinfo, _crvpulseinfomc, _crvwaveforminfo); + _crvHelper.FillCrvPulseInfoCollections(_crvRecoPulses, _crvDigiMCs, _crvDigis, + _crvpulseinfo, _crvpulseinfomc, _crvwaveforminfo); } } // fill this row in the TTree @@ -866,30 +810,6 @@ namespace mu2e { } } - // BestCrv branches - if (_crv && branchConfig.options().fillbestcrv()) { - // get the list of bestcrv modules - std::vector i_bestcrv_tags; - if (branchConfig.options().bestCrvModules(i_bestcrv_tags)) { // get the module labels - for (size_t i_bestCrvBranch = 0; i_bestCrvBranch < i_bestcrv_tags.size(); ++i_bestCrvBranch) { - auto hBestCrvAssns = _allBestCrvAssns.at(i_branch).at(i_bestCrvBranch); - if (hBestCrvAssns->size()>0) { - auto bestCrvCoinc = hBestCrvAssns->at(i_kseed).second; - _infoStructHelper.fillCrvHitInfo(bestCrvCoinc, _allBestCrvs.at(i_branch).at(i_bestCrvBranch)); - if (_fillmc) { - // loop through data-MC assns and find bestCrvCoinc (which is an art::Ptr){ - for (const auto& i_crvMCAssn : *_crvMCAssnsHandle) { - if (bestCrvCoinc == i_crvMCAssn.first) { - auto bestCrvCoincMC = i_crvMCAssn.second;//art::Ptr(_crvCoincidenceMCCollHandle, bestCrvCoinc.key()); - _infoMCStructHelper.fillCrvHitInfoMC(bestCrvCoincMC, _allBestCrvMCs.at(i_branch).at(i_bestCrvBranch)); - break; - } - } - } - } - } - } - } // all RecoQuals std::vector recoQuals; // for the output value @@ -1036,16 +956,11 @@ namespace mu2e { _allTSMIs.at(i_branch).clear(); _allTSHIMCs.at(i_branch).clear(); - if (_allBranches.at(i_branch).options().fillbestcrv()) { // only clear the vectors if they exist - // we don't want to remove elements so use assign instead of clear - _allBestCrvs.at(i_branch).assign(_allBestCrvs.at(i_branch).size(), CrvHitInfoReco()); - _allBestCrvMCs.at(i_branch).assign(_allBestCrvMCs.at(i_branch).size(), CrvHitInfoMC()); - } } // clear vectors - _crvinfo.clear(); - _crvinfomc.clear(); - _crvinfomcplane.clear(); + _crvhit.clear(); + _crvhitmc.clear(); + _crvhitmcplane.clear(); _crvpulseinfo.clear(); _crvwaveforminfo.clear(); _crvpulseinfomc.clear(); diff --git a/src/classes.h b/src/classes.h index 77a35d4..ef996b3 100644 --- a/src/classes.h +++ b/src/classes.h @@ -14,4 +14,11 @@ #include "TrkAna/inc/EventWeightInfo.hh" #include "TrkAna/inc/RecoQualInfo.hh" #include "TrkAna/inc/CaloClusterInfoMC.hh" -#include "TrkAna/inc/BestCrvAssns.hh" +#include "TrkAna/inc/CrvHitInfoReco.hh" +#include "TrkAna/inc/CrvHitInfoMC.hh" +#include "TrkAna/inc/CrvWaveformInfo.hh" +#include "TrkAna/inc/CrvPlaneInfoMC.hh" +#include "TrkAna/inc/CrvPulseInfoReco.hh" +#include "TrkAna/inc/CrvSummaryReco.hh" +#include "TrkAna/inc/CrvSummaryMC.hh" + diff --git a/src/classes_def.xml b/src/classes_def.xml index 954228c..345b02a 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -33,8 +33,16 @@ - - - - + + + + + + + + + + + + diff --git a/utils/TrkAnaUtils.C b/utils/TrkAnaUtils.C index d6d58b0..1304aa7 100644 --- a/utils/TrkAnaUtils.C +++ b/utils/TrkAnaUtils.C @@ -23,8 +23,8 @@ using namespace std; class TrkAnaUtils { public: - TrkAnaUtils(TTree* tatree); TrkAnaUtils(TFile* myfile,const char* treename="TrkAnaNeg"); + TrkAnaUtils(TTree* mytree); TrkAnaUtils(const char* filename,const char* treename); TrkAnaUtils(const char* treename="TrkAnaNeg"); TFile const& File() const { return *myfile_; } @@ -41,6 +41,7 @@ class TrkAnaUtils { void Draw(const char* lname,const char* cut="", const char* gopt="") const; void Project(const char* pname, const char* lname, const char* cut="") const; void Scan(const char* lname,const char* cut="") const; + void createEventList(const char* filename,TCut const& cut) const; private: void ListBranch(TBranch* branch, int idepth, int maxdepth) const; mutable TFile* myfile_; @@ -211,3 +212,16 @@ void TrkAnaUtils::Scan(const char* lname,const char* cut="") const { std::cout << "No current tree; call UseTree to set current tree" << std::endl; } } + +void TrkAnaUtils::createEventList(const char* filename,TCut const& cut) const { + // enable the output + TTreePlayer* player = (TTreePlayer*)(mytree_->GetPlayer()); + player->SetScanRedirect(true); + player->SetScanFileName(filename); + mytree_->SetScanField(0); + // create the output + mytree_->Scan("evtinfo.runid:evtinfo.subrunid:evtinfo.eventid","","col=#i:#i:#i"); + player->SetScanRedirect(false); + mytree_->SetScanField(50); +} +