diff --git a/CaloCluster/fcl/prolog.fcl b/CaloCluster/fcl/prolog.fcl index e69cf1d307..ce69466169 100644 --- a/CaloCluster/fcl/prolog.fcl +++ b/CaloCluster/fcl/prolog.fcl @@ -30,15 +30,29 @@ CaloCluster : { strategy : 0 diagLevel : 0 } + + CaloClusterFastMaker: { + module_type : "CaloClusterFast" + caloHitCollection : "CaloHitFastMaker:calo" + EminSeed : 10 + EnoiseCut : 1 + ExpandCut : 1 + deltaTime : 10 + diagLevel : 0 + extendSearch: true + minSiPMPerHit: 2 + timeOffset: -61 + } } CaloCluster : { @table::CaloCluster producers : { CaloProtoClusterMaker : { @table::CaloCluster.CaloProtoClusterMaker } CaloClusterMaker : { @table::CaloCluster.CaloClusterMaker } + CaloClusterFastMaker : { @table::CaloCluster.CaloClusterFastMaker } } - Reco : [ CaloProtoClusterMaker, CaloClusterMaker ] + Reco : [ CaloProtoClusterMaker, CaloClusterMaker, CaloClusterFastMaker ] } END_PROLOG diff --git a/CaloCluster/inc/ClusterUtils.hh b/CaloCluster/inc/ClusterUtils.hh index a968cb60bb..03a9ad296f 100644 --- a/CaloCluster/inc/ClusterUtils.hh +++ b/CaloCluster/inc/ClusterUtils.hh @@ -29,7 +29,6 @@ namespace mu2e { double e9 () const; double e25 () const; - private: const Calorimeter& cal_; const CaloHitPtrVector hits_; diff --git a/CaloCluster/src/ClusterUtils.cc b/CaloCluster/src/ClusterUtils.cc index 5f07a6b017..aecbae1b77 100644 --- a/CaloCluster/src/ClusterUtils.cc +++ b/CaloCluster/src/ClusterUtils.cc @@ -43,8 +43,7 @@ namespace mu2e { } - - + //------------------------------------------------------------------------------------------------- CLHEP::Hep3Vector ClusterUtils::cog3Vector() const { double sx(0),sy(0),sz(0),sx2(0),sy2(0),sw(0); @@ -55,13 +54,15 @@ namespace mu2e { return cal_.geomUtil().mu2eToDiskFF(iSection,cogMu2eFrame); } + //------------------------------------------------------------------------------------------------- double ClusterUtils::secondMoment() const { double sx(0),sy(0),sz(0),sx2(0),sy2(0),sw(0); fill(sx,sy,sz,sx2,sy2,sw); - return sx2-sx*sx/sw + sy2-sy*sy/sw; + return (sx2-sx*sx/sw + sy2-sy*sy/sw)/sw; } + //------------------------------------------------------------------------------------------------- void ClusterUtils::fill(double& sx, double& sy, double& sz, double& sx2, double& sy2, double& sw) const { int iSection = cal_.crystal(hits_[0]->crystalID()).diskID(); diff --git a/CaloDiag/CMakeLists.txt b/CaloDiag/CMakeLists.txt index cffccadbb9..b2141785ad 100644 --- a/CaloDiag/CMakeLists.txt +++ b/CaloDiag/CMakeLists.txt @@ -51,4 +51,16 @@ cet_build_plugin(CaloNNDiag art::module Offline::RecoDataProducts ) +cet_build_plugin(CaloNNTrain art::module + REG_SOURCE src/CaloNNTrain_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::CaloCluster + Offline::CalorimeterGeom + Offline::DataProducts + Offline::GeometryService + Offline::MCDataProducts + Offline::RecoDataProducts +) + install_source(SUBDIRS src) diff --git a/CaloDiag/fcl/prolog.fcl b/CaloDiag/fcl/prolog.fcl new file mode 100644 index 0000000000..99af18db3c --- /dev/null +++ b/CaloDiag/fcl/prolog.fcl @@ -0,0 +1,34 @@ +#include "Offline/Compression/fcl/prolog.fcl" +#include "Offline/CaloMC/fcl/prolog.fcl" +#include "Offline/CaloReco/fcl/prolog.fcl" +#include "Offline/CaloCluster/fcl/prolog.fcl" + +BEGIN_PROLOG +CaloNNDiag : { + module_type : CaloNNDiag + vdCollection : "compressDigiMCs:virtualdetector" + caloHitCollection : CaloHitMaker + caloClusterCollection : CaloClusterMaker + caloHitTruth : CaloHitTruthMatch + caloClusterTruth : CaloClusterTruthMatch + minCluEnergy : 0 +} + +CaloNNTrain: { + module_type : CaloNNTrain + caloClusterCollection : CaloClusterMaker + caloClusterMCCollection : CaloClusterTruthMatch + minEtoTest : 40 + MCEdepCut : 3 + diagLevel : 0 +} + + +CaloDiag : { + analyzers : { + CaloNNDiag : @local::CaloNNDiag + CaloNNTrain : @local::CaloNNTrain + } +} + +END_PROLOG diff --git a/CaloDiag/src/CaloNNDiag_module.cc b/CaloDiag/src/CaloNNDiag_module.cc index ba880c3767..e97c366d71 100644 --- a/CaloDiag/src/CaloNNDiag_module.cc +++ b/CaloDiag/src/CaloNNDiag_module.cc @@ -56,6 +56,7 @@ namespace mu2e { fhicl::Atom caloClusterCollection { Name("caloClusterCollection"), Comment("Calo cluster collection name") }; fhicl::Atom caloHitTruth { Name("caloHitTruth"), Comment("CaloHit truth name") }; fhicl::Atom caloClusterTruth { Name("caloClusterTruth"), Comment("caloCluster truth name") }; + fhicl::Atom minCluEnergy { Name("minCluEnergy"), Comment("Min cluster energy to write cluster") }; fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; }; @@ -68,34 +69,35 @@ namespace mu2e { private: - art::InputTag virtualDetectorTag_; - art::InputTag caloHitTag_; - art::InputTag caloClusterTag_; - art::InputTag caloHitTruthTag_; - art::InputTag caloClusterTruthTag_; - int diagLevel_; - int nProcess_; + art::InputTag virtualDetectorTag_; + art::InputTag caloHitTag_; + art::InputTag caloClusterTag_; + art::InputTag caloHitTruthTag_; + art::InputTag caloClusterTruthTag_; + float minCluEnergy_; + int diagLevel_; + int nProcess_; TTree* Ntup_; int _evt,_run; - int nHits_,cryId_[ntupLen],crySectionId_[ntupLen],crySimIdx_[ntupLen],crySimLen_[ntupLen]; - float cryTime_[ntupLen],cryEdep_[ntupLen],cryEdepErr_[ntupLen],cryPosX_[ntupLen],cryPosY_[ntupLen],cryPosZ_[ntupLen],_cryLeak[ntupLen]; + int nHits_,cryId_[ntupLen],crySectionId_[ntupLen],crySimIdx_[ntupLen],crySimLen_[ntupLen],cryConv_[ntupLen]; + float cryTime_[ntupLen],cryEdep_[ntupLen],cryEdepErr_[ntupLen],cryPosX_[ntupLen],cryPosY_[ntupLen],cryPosZ_[ntupLen]; + float cryTimeErr_[ntupLen],cryT1_[ntupLen],cryT2_[ntupLen],cryT1Err_[ntupLen],cryT2Err_[ntupLen]; - int nSimHit_,crySimId_[ntupLen],crySimPdgId_[ntupLen],crySimCrCode_[ntupLen],crySimGenIdx_[ntupLen],cryConv_[ntupLen]; - float crySimMom_[ntupLen],crySimStartX_[ntupLen],crySimStartY_[ntupLen],crySimStartZ_[ntupLen],crySimStartT_[ntupLen]; - float crySimTime_[ntupLen],crySimEdep_[ntupLen],cryTimeErr_[ntupLen],cryT1_[ntupLen],cryT2_[ntupLen],cryT1Err_[ntupLen],cryT2Err_[ntupLen]; + std::vector> crySimId_,crySimPdgId_,crySimCrCode_; + std::vector> crySimMom_,crySimTime_,crySimEdep_,crySimStartX_,crySimStartY_,crySimStartZ_; int nCluster_,nCluSim_,cluNcrys_[ntupLen],cluDisk_[ntupLen]; float cluEnergy_[ntupLen],cluEnergyErr_[ntupLen],cluTime_[ntupLen],cluTimeErr_[ntupLen],cluCogX_[ntupLen],cluCogY_[ntupLen],cluCogR_[ntupLen], cluCogZ_[ntupLen],cluE1_[ntupLen],cluE2_[ntupLen],cluE9_[ntupLen],cluE25_[ntupLen],cluSecMom_[ntupLen],cluDR_[ntupLen]; int cluSplit_[ntupLen],cluConv_[ntupLen],cluSimIdx_[ntupLen],cluSimLen_[ntupLen]; - std::vector > cluList_; - int cluSimId_[ntupLen],cluSimPdgId_[ntupLen],cluSimGenId_[ntupLen],cluSimGenPdg_[ntupLen],cluSimCrCode_[ntupLen]; - float cluSimMom_[ntupLen],cluSimMom2_[ntupLen],cluSimPosX_[ntupLen],cluSimPosY_[ntupLen],cluSimPosZ_[ntupLen],cluSimStartX_[ntupLen], - cluSimStartY_[ntupLen],cluSimStartZ_[ntupLen],cluSimTime_[ntupLen],cluSimEdep_[ntupLen]; + std::vector > cluList_; + std::vector> cluSimId_,cluSimPdgId_,cluSimCrCode_; + std::vector> cluSimMom_,cluSimMom2_,cluSimTime_,cluSimEdep_,cluSimPosX_,cluSimPosY_,cluSimPosZ_; + std::vector> cluSimStartX_,cluSimStartY_,cluSimStartZ_; int nVd_,vdId_[ntupLen],vdPdgId_[ntupLen],vdGenId_[ntupLen],vdGenIdx_[ntupLen]; float vdTime_[ntupLen],vdPosX_[ntupLen],vdPosY_[ntupLen],vdPosZ_[ntupLen],vdMom_[ntupLen],vdMomX_[ntupLen],vdMomY_[ntupLen],vdMomZ_[ntupLen]; @@ -109,6 +111,7 @@ namespace mu2e { caloClusterTag_ (config().caloClusterCollection()), caloHitTruthTag_ (config().caloHitTruth()), caloClusterTruthTag_(config().caloClusterTruth()), + minCluEnergy_ (config().minCluEnergy()), diagLevel_ (config().diagLevel()), nProcess_(0), Ntup_(0) @@ -138,21 +141,15 @@ namespace mu2e { Ntup_->Branch("cryT1Err", &cryT1Err_ , "cryT1Err[nCry]/F"); Ntup_->Branch("cryT2Err", &cryT2Err_ , "cryT2Err[nCry]/F"); Ntup_->Branch("cryConv", &cryConv_ , "cryConv[nCry]/I"); - - Ntup_->Branch("crySimIdx", &crySimIdx_ , "crySimIdx[nCry]/I"); - Ntup_->Branch("crySimLen", &crySimLen_ , "crySimLen[nCry]/I"); - Ntup_->Branch("nSim", &nSimHit_ , "nSim/I"); - Ntup_->Branch("simId", &crySimId_ , "simId[nSim]/I"); - Ntup_->Branch("simPdgId", &crySimPdgId_ , "simPdgId[nSim]/I"); - Ntup_->Branch("simCrCode", &crySimCrCode_ ,"simCrCode[nSim]/I"); - Ntup_->Branch("simMom", &crySimMom_ , "simMom[nSim]/F"); - Ntup_->Branch("simStartX", &crySimStartX_ ,"simStartX[nSim]/F"); - Ntup_->Branch("simStartY", &crySimStartY_ ,"simStartY[nSim]/F"); - Ntup_->Branch("simStartZ", &crySimStartZ_ ,"simStartZ[nSim]/F"); - Ntup_->Branch("simStartT", &crySimStartT_ ,"simStartT[nSim]/F"); - Ntup_->Branch("simTime", &crySimTime_ , "simTime[nSim]/F"); - Ntup_->Branch("simEdep", &crySimEdep_ , "simEdep[nSim]/F"); - Ntup_->Branch("simGenIdx", &crySimGenIdx_ ,"simGenIdx[nSim]/I"); + Ntup_->Branch("crySimId", &crySimId_); + Ntup_->Branch("crySimPdgId", &crySimPdgId_); + Ntup_->Branch("crySimCrCode", &crySimCrCode_); + Ntup_->Branch("crySimMom", &crySimMom_); + Ntup_->Branch("crySimTime", &crySimTime_); + Ntup_->Branch("crySimEdep", &crySimEdep_); + Ntup_->Branch("crySimStartX", &crySimStartX_); + Ntup_->Branch("crySimStartY", &crySimStartY_); + Ntup_->Branch("crySimStartZ", &crySimStartZ_); Ntup_->Branch("nCluster", &nCluster_ , "nCluster/I"); Ntup_->Branch("cluEnergy", &cluEnergy_ , "cluEnergy[nCluster]/F"); @@ -169,30 +166,24 @@ namespace mu2e { Ntup_->Branch("cluE2", &cluE2_ , "cluE2[nCluster]/F"); Ntup_->Branch("cluE9", &cluE9_ , "cluE9[nCluster]/F"); Ntup_->Branch("cluE25", &cluE25_ , "cluE25[nCluster]/F"); - Ntup_->Branch("cluDR", &cluDR_ , "cluDR[nCluster]/F"); Ntup_->Branch("cluSecMom", &cluSecMom_ , "cluSecMom[nCluster]/F"); Ntup_->Branch("cluSplit", &cluSplit_ , "cluSplit[nCluster]/I"); Ntup_->Branch("cluConv", &cluConv_ , "cluConv[nCluster]/I"); - Ntup_->Branch("cluSimIdx", &cluSimIdx_ , "cluSimIdx[nCluster]/I"); - Ntup_->Branch("cluSimLen", &cluSimLen_ , "cluSimLen[nCluster]/I"); Ntup_->Branch("cluList", &cluList_); - Ntup_->Branch("nCluSim", &nCluSim_ , "nCluSim/I"); - Ntup_->Branch("cluSimId", &cluSimId_ , "cluSimId[nCluSim]/I"); - Ntup_->Branch("cluSimPdgId", &cluSimPdgId_ , "cluSimPdgId[nCluSim]/I"); - Ntup_->Branch("cluSimGenId", &cluSimGenId_ , "cluSimGenId[nCluSim]/I"); - Ntup_->Branch("cluSimGenPdg", &cluSimGenPdg_, "cluSimGenPdg[nCluSim]/I"); - Ntup_->Branch("cluSimCrCode", &cluSimCrCode_ ,"cluSimCrCode[nCluSim]/I"); - Ntup_->Branch("cluSimMom", &cluSimMom_ , "cluSimMom[nCluSim]/F"); - Ntup_->Branch("cluSimMom2", &cluSimMom2_ , "cluSimMom2[nCluSim]/F"); - Ntup_->Branch("cluSimPosX", &cluSimPosX_ , "cluSimPosX[nCluSim]/F"); - Ntup_->Branch("cluSimPosY", &cluSimPosY_ , "cluSimPosY[nCluSim]/F"); - Ntup_->Branch("cluSimPosZ", &cluSimPosZ_ , "cluSimPosZ[nCluSim]/F"); - Ntup_->Branch("cluSimStartX", &cluSimStartX_ ,"cluSimStartX[nCluSim]/F"); - Ntup_->Branch("cluSimStartY", &cluSimStartY_ ,"cluSimStartY[nCluSim]/F"); - Ntup_->Branch("cluSimStartZ", &cluSimStartZ_ ,"cluSimStartZ[nCluSim]/F"); - Ntup_->Branch("cluSimTime", &cluSimTime_ , "cluSimTime[nCluSim]/F"); - Ntup_->Branch("cluSimEdep", &cluSimEdep_ , "cluSimEdep[nCluSim]/F"); + Ntup_->Branch("cluSimId", &cluSimId_); + Ntup_->Branch("cluSimPdgId", &cluSimPdgId_); + Ntup_->Branch("cluSimCrCode", &cluSimCrCode_); + Ntup_->Branch("cluSimMom", &cluSimMom_); + Ntup_->Branch("cluSimMom2", &cluSimMom2_); + Ntup_->Branch("cluSimTime", &cluSimTime_); + Ntup_->Branch("cluSimEdep", &cluSimEdep_); + Ntup_->Branch("cluSimPosX", &cluSimPosX_); + Ntup_->Branch("cluSimPosY", &cluSimPosY_); + Ntup_->Branch("cluSimPosZ", &cluSimPosZ_); + Ntup_->Branch("cluSimStartX", &cluSimStartX_); + Ntup_->Branch("cluSimStartY", &cluSimStartY_); + Ntup_->Branch("cluSimStartZ", &cluSimStartZ_); Ntup_->Branch("nVd", &nVd_ , "nVd/I"); Ntup_->Branch("vdId", &vdId_ , "vdId[nVd]/I"); @@ -226,11 +217,13 @@ namespace mu2e { //Calorimeter crystal hits (average from readouts) art::Handle CaloHitsHandle; event.getByLabel(caloHitTag_, CaloHitsHandle); + if (!CaloHitsHandle.isValid()) return; const CaloHitCollection& CaloHits(*CaloHitsHandle); //Calorimeter clusters art::Handle caloClustersHandle; event.getByLabel(caloClusterTag_, caloClustersHandle); + if (!caloClustersHandle.isValid()) return; const CaloClusterCollection& caloClusters(*caloClustersHandle); //Virtual detector hits @@ -240,14 +233,23 @@ namespace mu2e { //Calo digi truth assignment art::Handle caloDigiTruthHandle; event.getByLabel(caloHitTruthTag_, caloDigiTruthHandle); + if (!caloDigiTruthHandle.isValid()) return; const CaloHitMCTruthAssn& caloDigiTruth(*caloDigiTruthHandle); //Calo cluster truth assignment art::Handle caloClusterTruthHandle; event.getByLabel(caloClusterTruthTag_, caloClusterTruthHandle); + if (!caloClusterTruthHandle.isValid()) return; const CaloClusterMCTruthAssn& caloClusterTruth(*caloClusterTruthHandle); + //SELECT EVENTS WITH AT LEAST A CLUSTER ABOVE MINCLUENERGY_ + bool select(false); + for (const auto& cluster : caloClusters) { + if (cluster.energyDep()>minCluEnergy_) {select=true; break;} + } + if (!select) return; + std::map, const StepPointMC*> vdMap; @@ -271,10 +273,11 @@ namespace mu2e { //-------------------------- Do calorimeter hits -------------------------------- - nHits_ = nSimHit_ = 0; + nHits_ = 0; + crySimId_.clear();crySimPdgId_.clear();crySimCrCode_.clear();crySimTime_.clear();crySimEdep_.clear(); + crySimMom_.clear();crySimStartX_.clear();crySimStartY_.clear();crySimStartZ_.clear(); - for (unsigned int ic=0; ic sId,sPdg,sCr; + std::vectorst,sed,sm,stx,sty,stz; for (unsigned i=0;i< nCrySims;++i) { const auto& eDepMC = itMC->second->energyDeposit(i); auto parent(eDepMC.sim()); while (parent->hasParent()) parent = parent->parent(); - int genId=-1; - if (parent->genParticle()) genId = parent->genParticle()->generatorId().id(); - - crySimId_[nSimHit_] = eDepMC.sim()->id().asInt(); - crySimPdgId_[nSimHit_] = eDepMC.sim()->pdgId(); - crySimCrCode_[nSimHit_] = eDepMC.sim()->creationCode(); - crySimTime_[nSimHit_] = eDepMC.time(); - crySimEdep_[nSimHit_] = eDepMC.energyDep(); - crySimMom_[nSimHit_] = eDepMC.momentumIn(); - crySimStartX_[nSimHit_] = parent->startPosition().x(); - crySimStartY_[nSimHit_] = parent->startPosition().y(); - crySimStartZ_[nSimHit_] = parent->startPosition().z(); - crySimStartT_[nSimHit_] = parent->startGlobalTime(); - crySimGenIdx_[nSimHit_] = genId; - ++nSimHit_; + + sId.push_back(eDepMC.sim()->id().asInt()); + sPdg.push_back(eDepMC.sim()->pdgId()); + sCr.push_back(eDepMC.sim()->creationCode()); + sm.push_back(eDepMC.momentumIn()); + sed.push_back(eDepMC.energyDep()); + st.push_back(eDepMC.time()); + stx.push_back(parent->startPosition().x()); + sty.push_back(parent->startPosition().y()); + stz.push_back(parent->startPosition().z()); } + crySimId_.push_back(sId); + crySimPdgId_.push_back(sPdg); + crySimCrCode_.push_back(sCr); + crySimMom_.push_back(sm); + crySimEdep_.push_back(sed); + crySimTime_.push_back(st); + crySimStartX_.push_back(stx); + crySimStartY_.push_back(sty); + crySimStartZ_.push_back(stz); + ++nHits_; } //-------------------------- Do clusters -------------------------------- - nCluster_ = nCluSim_ = 0; + nCluster_ = 0; cluList_.clear(); - - - //find the most energetic CE cluster - unsigned icMCIdx(-1); double convEnergy(0); - for (unsigned ic=0; icfirst.get() == &cluster) break; ++itMC;} - - if (itMC != caloClusterTruth.end()) - { - for (auto& edep : itMC->second->energyDeposits()) - { - auto parent(edep.sim()); - while (parent->hasParent()) parent = parent->parent(); - if (parent->genParticle() && parent->genParticle()->generatorId().isConversion() && cluster.energyDep() > convEnergy) - { - convEnergy = cluster.energyDep(); - icMCIdx = ic; - } - } - } - } - - + cluSimId_.clear();cluSimPdgId_.clear();cluSimCrCode_.clear();cluSimTime_.clear();cluSimEdep_.clear(); + cluSimMom_.clear();cluSimMom2_.clear();cluSimPosX_.clear();cluSimPosY_.clear();cluSimPosZ_.clear(); + cluSimStartX_.clear();cluSimStartY_.clear();cluSimStartZ_.clear(); for (unsigned ic=0; ic cryList; for (auto cryPtr : cluster.caloHitsPtrVector()) cryList.push_back(std::distance(&CaloHits.at(0),cryPtr.get())); ClusterUtils cluUtil(cal, cluster); - auto cog = cluUtil.cog3Vector(); auto itMC = caloClusterTruth.begin(); while (itMC != caloClusterTruth.end()) {if (itMC->first.get() == &cluster) break; ++itMC;} const auto eDepMCs = (itMC != caloClusterTruth.end()) ? itMC->second->energyDeposits() : std::vector{}; - const CaloHit& seedHit = CaloHits.at(cryList[0]); - CLHEP::Hep3Vector seedPos = cal.geomUtil().mu2eToDiskFF(cluster.diskID(),cal.crystal(seedHit.crystalID()).position()); - - double dr(0); - if (cryList.size()>1) - { - const CaloHit& hit = CaloHits.at(cryList[1]); - CLHEP::Hep3Vector hitPos = cal.geomUtil().mu2eToDiskFF(cluster.diskID(),cal.crystal(hit.crystalID()).position()); - double dx = hitPos.x()-seedPos.x(); - double dy = hitPos.y()-seedPos.y(); - dr = sqrt(dx*dx+dy*dy); + bool isConv(false); + if (itMC != caloClusterTruth.end()){ + for (auto& edep : itMC->second->energyDeposits()){ + if (edep.sim()->creationCode() == ProcessCode::mu2eCeMinusEndpoint || + edep.sim()->creationCode() == ProcessCode::mu2eCePlusEndpoint || + edep.sim()->creationCode() == ProcessCode::mu2eCeMinusLeadingLog || + edep.sim()->creationCode() == ProcessCode::mu2eCePlusLeadingLog){ + isConv = true; + break; + } + } } + const CaloHit& seedHit = CaloHits.at(cryList[0]); + CLHEP::Hep3Vector seedPos = cal.geomUtil().mu2eToDiskFF(cluster.diskID(),cal.crystal(seedHit.crystalID()).position()); cluDisk_[nCluster_] = cluster.diskID(); cluEnergy_[nCluster_] = cluster.energyDep(); @@ -421,15 +408,14 @@ namespace mu2e { cluE2_[nCluster_] = cluUtil.e2(); cluE9_[nCluster_] = cluUtil.e9(); cluE25_[nCluster_] = cluUtil.e25(); - cluDR_[nCluster_] = dr; cluSecMom_[nCluster_] = cluUtil.secondMoment(); cluSplit_[nCluster_] = cluster.isSplit(); - cluConv_[nCluster_] = ic == icMCIdx; + cluConv_[nCluster_] = isConv;//ic == icMCIdx; cluList_.push_back(cryList); - cluSimIdx_[nCluster_] = nCluSim_; - cluSimLen_[nCluster_] = eDepMCs.size(); + std::vector sId,sPdg,sCr; + std::vectorst,sed,sm,sm2,spx,spy,spz,stx,sty,stz; for (unsigned i=0;i< eDepMCs.size();++i) { const auto& eDepMC = eDepMCs[i]; @@ -437,11 +423,6 @@ namespace mu2e { art::Ptr smother(sim); while (smother->hasParent() && !smother->genParticle() ) smother = smother->parent(); - int genId=-1; - if (smother->genParticle()) genId = smother->genParticle()->generatorId().id(); - int genPdg=-1; - if (smother->genParticle()) genPdg = smother->genParticle()->pdgId(); - double simMom(-1); CLHEP::Hep3Vector simPos(0,0,0); @@ -449,27 +430,37 @@ namespace mu2e { if (vdMapEntry != vdMap.end()) { simMom = vdMapEntry->second->momentum().mag(); - CLHEP::Hep3Vector simPos = cal.geomUtil().mu2eToDiskFF(cluster.diskID(), vdMapEntry->second->position()); + simPos = cal.geomUtil().mu2eToDiskFF(cluster.diskID(), vdMapEntry->second->position()); } - cluSimId_[nCluSim_] = sim->id().asInt(); - cluSimPdgId_[nCluSim_] = sim->pdgId(); - cluSimCrCode_[nCluSim_] = sim->creationCode(); - cluSimGenId_[nCluSim_] = genId; - cluSimGenPdg_[nCluSim_] = genPdg; - cluSimTime_[nCluSim_] = eDepMC.time(); - cluSimEdep_[nCluSim_] = eDepMC.energyDep(); - cluSimMom_[nCluSim_] = eDepMC.momentumIn(); - cluSimMom2_[nCluSim_] = simMom; - cluSimPosX_[nCluSim_] = simPos.x(); // in disk FF frame - cluSimPosY_[nCluSim_] = simPos.y(); - cluSimPosZ_[nCluSim_] = simPos.z(); - cluSimStartX_[nCluSim_] = sim->startPosition().x(); // in disk FF frame - cluSimStartY_[nCluSim_] = sim->startPosition().y(); - cluSimStartZ_[nCluSim_] = sim->startPosition().z(); - - ++nCluSim_; + sId.push_back(sim->id().asInt()); + sPdg.push_back(sim->pdgId()); + sCr.push_back(sim->creationCode()); + st.push_back(eDepMC.time()); + sed.push_back(eDepMC.energyDep()); + sm.push_back(eDepMC.momentumIn()); + sm2.push_back(simMom); + spx.push_back(simPos.x()); // in disk FF frame + spy.push_back(simPos.y()); + spz.push_back(simPos.z()); + stx.push_back(sim->startPosition().x()); // in disk FF frame + sty.push_back(sim->startPosition().y()); + stz.push_back(sim->startPosition().z()); } + cluSimId_.push_back(sId); + cluSimPdgId_.push_back(sPdg); + cluSimCrCode_.push_back(sCr); + cluSimMom_.push_back(sm); + cluSimMom2_.push_back(sm2); + cluSimEdep_.push_back(sed); + cluSimTime_.push_back(st); + cluSimPosX_.push_back(spx); + cluSimPosY_.push_back(spy); + cluSimPosZ_.push_back(spz); + cluSimStartX_.push_back(stx); + cluSimStartY_.push_back(sty); + cluSimStartZ_.push_back(stz); + ++nCluster_; } diff --git a/CaloDiag/src/CaloNNTrain_module.cc b/CaloDiag/src/CaloNNTrain_module.cc new file mode 100644 index 0000000000..85d3449d0f --- /dev/null +++ b/CaloDiag/src/CaloNNTrain_module.cc @@ -0,0 +1,145 @@ +#include "art_root_io/TFileService.h" +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "cetlib_except/exception.h" +#include "fhiclcpp/types/Atom.h" + +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CaloCluster/inc/ClusterUtils.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/MCDataProducts/inc/CaloMCTruthAssns.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +#include "TDirectory.h" +#include "TNtuple.h" +#include "TTree.h" + +#include +#include + + + +namespace mu2e { + + class CaloNNTrain : public art::EDAnalyzer + { + public: + struct Config + { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom caloClusterCollection { Name("caloClusterCollection"), Comment("Calo cluster collection name") }; + fhicl::Atom caloClusterMCCollection { Name("caloClusterMCCollection"), Comment("Calo cluster MC collection name") }; + fhicl::Atom minEtoTest { Name("minEtoTest"), Comment("Minimum Energy to run the MVA") }; + fhicl::Atom MCEdepCut { Name("MCEdepCut"), Comment("Min E cut for MC contribution to cluster") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; + }; + + explicit CaloNNTrain(const art::EDAnalyzer::Table& config) : + EDAnalyzer{config}, + caloClusterToken_ {consumes (config().caloClusterCollection())}, + caloClusterMCToken_{consumes(config().caloClusterMCCollection())}, + minEtoTest_ (config().minEtoTest()), + MCEdepCut_ (config().MCEdepCut()), + diagLevel_ (config().diagLevel()) + {} + + virtual void analyze(const art::Event& event) override; + virtual void beginJob() override; + + + private: + art::ProductToken caloClusterToken_; + art::ProductToken caloClusterMCToken_; + float minEtoTest_; + float MCEdepCut_; + int diagLevel_; + + TTree* Ntup_; + int evt_,cluNcrys_,cluDisk_; + float cluEnergy_,cluTime_,cluCogR_,cluE1_,cluE2_,cluE9_,cluE25_,cluSec_,cluMCEtot_; + std::vector cluMCEPdg_,cluMCCr_; + std::vector cluMCEdep_; + }; + + + void CaloNNTrain::beginJob() + { + art::ServiceHandle tfs; + Ntup_ = tfs->make("Calo", "Calo"); + Ntup_->Branch("evt", &evt_ , "evt/I"); + Ntup_->Branch("cluNcrys", &cluNcrys_ , "cluNcrys/I"); + Ntup_->Branch("cluDisk", &cluDisk_ , "cluDisk/I"); + Ntup_->Branch("cluEnergy", &cluEnergy_ , "cluEnergy/F"); + Ntup_->Branch("cluTime", &cluTime_ , "cluTime/F"); + Ntup_->Branch("cluCogR", &cluCogR_ , "cluCogR/F"); + Ntup_->Branch("cluE1", &cluE1_ , "cluE1/F"); + Ntup_->Branch("cluE2", &cluE2_ , "cluE2/F"); + Ntup_->Branch("cluE9", &cluE9_ , "cluE9/F"); + Ntup_->Branch("cluE25", &cluE25_ , "cluE25/F"); + Ntup_->Branch("cluSec", &cluSec_ , "cluSec/F"); + Ntup_->Branch("cluMCEtot", &cluMCEtot_ , "cluMCEtot/F"); + Ntup_->Branch("cluMCEdep", &cluMCEdep_); + Ntup_->Branch("cluMCEPdg", &cluMCEPdg_); + Ntup_->Branch("cluMCCr", &cluMCCr_); + } + + + void CaloNNTrain::analyze(const art::Event& event) + { + evt_ = event.id().event(); + art::Handle caloClustersHandle = event.getHandle(caloClusterToken_); + art::Handle caloClustersMCHandle = event.getHandle(caloClusterMCToken_); + + if (!caloClustersHandle.isValid()) return; + if (!caloClustersMCHandle.isValid()) return; + + const Calorimeter& cal = *(GeomHandle()); + const CaloClusterCollection& caloClusters(*caloClustersHandle); + const CaloClusterMCTruthAssn& caloClustersMC(*caloClustersMCHandle); + + for (const auto& cluster : caloClusters) + { + if (cluster.energyDep() < minEtoTest_) continue; + + //MC analyis + float MCEdepTot(0); + cluMCEdep_.clear();cluMCEPdg_.clear();cluMCCr_.clear(); + + auto itMC = caloClustersMC.begin(); + while (itMC != caloClustersMC.end()) {if (itMC->first.get() == &cluster) break; ++itMC;} + const auto eDepMCs = (itMC != caloClustersMC.end()) ? itMC->second->energyDeposits() : std::vector{}; + + if (itMC != caloClustersMC.end()){ + MCEdepTot= itMC->second->totalEnergyDep(); + for (auto& edep : itMC->second->energyDeposits()) { + if (edep.energyDep() < MCEdepCut_) continue; + cluMCEdep_.push_back(edep.energyDep()); + cluMCEPdg_.push_back(edep.sim()->pdgId()); + cluMCCr_.push_back(edep.sim()->creationCode()); + } + } + + ClusterUtils cluUtil(cal, cluster); + + cluNcrys_ = cluster.size(); + cluDisk_ = cluster.diskID(); + cluEnergy_ = cluster.energyDep(); + cluTime_ = cluster.time(); + cluCogR_ = cluster.cog3Vector().perp(); + cluE1_ = cluUtil.e1()/cluster.energyDep(); + cluE2_ = cluUtil.e2()/cluster.energyDep(); + cluE9_ = cluUtil.e9()/cluster.energyDep(); + cluE25_ = cluUtil.e25()/cluster.energyDep(); + cluSec_ = cluUtil.secondMoment(); + cluMCEtot_ = MCEdepTot; + + Ntup_->Fill(); + } + } +} + +DEFINE_ART_MODULE(mu2e::CaloNNTrain) diff --git a/CaloFilters/CMakeLists.txt b/CaloFilters/CMakeLists.txt index cc95ca6973..f2b80fbab4 100644 --- a/CaloFilters/CMakeLists.txt +++ b/CaloFilters/CMakeLists.txt @@ -49,8 +49,8 @@ cet_build_plugin(FilterEcalMVATrigger art::module Offline::RecoDataProducts ) -cet_build_plugin(FilterEcalNNTrigger art::module - REG_SOURCE src/FilterEcalNNTrigger_module.cc +cet_build_plugin(CaloNNEval art::module + REG_SOURCE src/CaloNNEval_module.cc LIBRARIES REG Offline::CalorimeterGeom Offline::GeometryService @@ -58,9 +58,17 @@ cet_build_plugin(FilterEcalNNTrigger art::module Offline::RecoDataProducts ) +cet_build_plugin(FilterEcalNNTrigger art::module + REG_SOURCE src/FilterEcalNNTrigger_module.cc + LIBRARIES REG + Offline::RecoDataProducts +) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/ce_bkg_20_BDT.weights.xml ${CURRENT_BINARY_DIR} data/ce_bkg_20_BDT.weights.xml) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/ce_bkg_ecal_20_BDT.weights.xml ${CURRENT_BINARY_DIR} data/ce_bkg_ecal_20_BDT.weights.xml) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/CE_NN_ReLU.weights.xml ${CURRENT_BINARY_DIR} data/CE_NN_ReLU.weights.xml) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/CE_NN_ReLU_noET.weights.xml ${CURRENT_BINARY_DIR} data/CE_NN_ReLU_noET.weights.xml) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/FlatPhoton_ReLU_noET.weights.xml ${CURRENT_BINARY_DIR} data/FlatPhoton_ReLU_noET.weights.xml) install(DIRECTORY data DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/Offline/CaloFilters) diff --git a/CaloFilters/data/CE_NN_ReLU.weights.xml b/CaloFilters/data/CE_NN_ReLU.weights.xml index b7116aadf4..765e38f09c 100644 --- a/CaloFilters/data/CE_NN_ReLU.weights.xml +++ b/CaloFilters/data/CE_NN_ReLU.weights.xml @@ -2,17 +2,17 @@ - + - - - - - + + + + + - + @@ -45,15 +45,16 @@ - - - - - - - - - + + + + + + + + + + @@ -63,61 +64,66 @@ - + + - + - + + - + - - - - - - - - + + + + + + + + + - - - - - - - - + + + + + + + + + - - - - - - - - + + + + + + + + + @@ -125,62 +131,68 @@ - - - 2.8897795944819955e+00 2.5104566021939445e+00 2.1457936397983324e+00 1.7048419191123221e+00 -2.2486157609628603e+00 -1.5647067228141165e+00 -1.7636261356482714e+00 2.0141966735648476e+00 + + + -4.1642418919049700e+00 2.7645176241795246e+00 1.2794422280626394e+00 2.4712413865124412e+00 -7.3270015430612387e+00 1.9360689763290364e+00 2.3811432281522396e+00 2.0936920066048361e+00 -1.3270824675089050e+00 - - 1.0998059727997564e+00 8.7961583177187741e-01 -2.8624732453804640e+00 -4.2983185042366590e-01 -3.6342359953202868e-02 -1.2786862242180601e+00 -2.0244226802661838e+00 7.3494724059290872e-01 + + -3.2300119855515246e+00 -3.1738441595990778e+00 -7.7194705835438332e-02 -6.3111718655345612e+00 -6.3279875230191476e-02 2.3101163597291285e+00 2.8206948001773857e+00 -5.0862380152545739e-01 1.6963640171126451e+00 - - 9.2003226739360011e-01 1.4213913790687025e+00 -8.2021409848535154e-01 1.3932213321106683e+00 1.4602712349186537e+00 -1.3889168774174112e+00 1.4209328711037428e+00 4.5167241262119170e-02 + + -1.2101791029389930e+00 3.0496153459223723e+00 3.5168706140026684e-01 -7.0756892061737107e+00 4.2312715833513548e-01 3.4361327893526714e+00 -6.7964245988020089e-02 -2.9025998035617490e-01 1.4606468835713085e+00 - - -3.8769480756134080e+00 -1.3871498901864074e+00 5.9953583026955537e-01 -1.0413054952287966e+00 5.2964969892507013e-02 9.2904108739729918e-01 1.1649139329580944e+00 1.6557254880210439e+00 + + 3.9703813710811631e-01 -6.9309507102143875e-01 -2.0277505439066429e-01 -8.9640770369030032e-01 3.1578385743649644e+00 1.5770507622333201e+00 -3.1674491827676530e+00 4.3911891943572967e-01 -2.3107431651618954e-01 - - 1.8705324584403085e+00 -4.9547688312002153e-02 2.4153002128106666e+00 -2.6735776384365839e-01 -3.4327047424104489e+00 1.5296122367790568e+00 -1.5079083121908514e+00 1.1376706315566725e+00 + + -9.3827570695613344e-02 -8.6070492639341678e-01 -1.8327139295707826e+00 -9.1290823410556732e-01 -1.1290397149903059e-01 -6.8107040927625517e-01 1.1569161730281838e+00 1.1068070757770845e+00 -6.9115961211822241e-01 - - -1.3338162408965160e+00 8.9836060953298780e-01 9.2274760429901748e-01 1.2152235739597372e+00 -1.2393507823454075e+00 -1.4906180224127805e+00 6.2626340648151868e-01 8.0414418976465341e-02 + + 2.1621285611973096e-01 -2.9290991045512871e-01 -2.0767387227784967e-01 -1.9547612799583874e-01 -1.2971478669533756e+00 -8.4483402500461513e-01 5.5596011029380898e-01 -1.8933762441625630e+00 1.3069590162468725e+00 - - 1.0203207262354812e+00 2.3119723991326104e+00 2.3169808800717995e+00 1.3019763790848873e+00 7.4513711852358633e-01 1.5082065040697887e+00 -2.8514457157970750e+00 -8.9683871008565574e-01 + + -2.4695728104297525e+00 5.3882159027560983e-01 -9.0844166472807375e-01 -1.1288029047826229e-01 1.2384133566007161e+00 3.0390961224429286e-01 1.0728195394251061e+00 2.0828468603912970e-01 -2.7575815311376545e-01 - - 4.9882418442364668e-01 -1.5236848876205535e+00 1.0669649844776627e+00 6.1111406945463942e-01 2.3583520881740413e-01 8.1768194367920299e-01 6.5426499128950277e-01 7.7304152303739671e-01 + + -2.4504142989956695e-01 2.0002365521970622e-01 1.2062453171669854e-01 5.8374790568077617e-01 1.2379520318934836e+00 1.4296347564101666e+00 9.3890873437796640e-01 -7.6415792532915616e-01 -1.4019018675094435e+00 - - 1.3077975563638427e+00 -5.7076923465350948e-01 -1.0354610312543981e+00 -1.2831392761380016e+00 5.7197278781685112e-01 -1.8504495160794860e+00 -1.7199980555095071e+00 -1.3432148831297448e+00 + + 1.2521008238869695e+00 1.8002829431301337e-01 2.0546220057046196e+00 5.1261693213073722e-01 -6.9251130428832630e-01 -5.1451860308551627e+00 6.6540032674205185e-01 5.4340486731223436e-01 2.0279231304087402e+00 + + + -1.7904585823358374e+00 1.4396426950171270e+00 -4.9901324709834366e-01 -7.7038955964190237e+00 -4.3047184010750303e+00 2.7915819408243667e+00 -9.2322776185295929e-01 -1.5998252072304731e+00 -1.8229957973662343e+00 - + + + -8.6158119913364828e-01 + - 1.8511536272246372e+00 + -8.6997743019507146e-01 - 2.5243429176175352e+00 + 2.3981975128116382e+00 - 3.7011547443920016e+00 + -1.2457566034425143e+00 - -2.0691181968701661e-01 + -9.9285650355948896e-01 - -1.5898787200356739e+00 + 1.1908740557342365e+00 - -1.1822992448160146e+00 + -1.0909248141131882e+00 - 1.3233799326226614e+00 + -1.8916380320626034e+00 - 5.5350300555079823e-01 + -1.4498296503183823e+00 - 1.1279142088093761e+00 + 2.3333924914467783e+00 diff --git a/CaloFilters/data/CE_NN_ReLU_noET.weights.xml b/CaloFilters/data/CE_NN_ReLU_noET.weights.xml new file mode 100644 index 0000000000..27b1665601 --- /dev/null +++ b/CaloFilters/data/CE_NN_ReLU_noET.weights.xml @@ -0,0 +1,179 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 6.6457789549253776e-01 3.9745383736480511e+00 5.2579520330341225e+00 8.9864568504144082e-01 -9.9927584888235828e+00 -4.2233060506304049e-01 -1.1391530757812680e+00 + + + -5.9907367073300910e-01 -4.8040947782253196e-01 2.3099633100001182e+00 -1.9030861753689639e+00 -1.2343265450823999e+00 -1.8423760724888383e+00 9.3950975405978521e-01 + + + 2.0706163238781730e+00 -2.0960333236732789e-03 -4.4932389070498457e-01 3.7391208418286048e+00 -1.5689407333447762e+00 1.9108535274361254e+00 1.4280742628463554e+00 + + + 3.0651560858249477e-01 1.8767394906515302e+00 9.0237280860621771e-01 -1.3999823713990713e+00 7.4556183738118853e-01 -3.3719959654437459e-01 -1.8377545622477682e+00 + + + -3.2621373217605023e+00 3.1628371101550896e+00 3.1744538235674380e-01 1.1917900797971521e+00 2.3526168183602750e-01 -6.9160460664308332e-01 -1.7545723718713013e+00 + + + -7.4094287926734981e-01 2.6377018648696499e-01 -1.4901989081501810e-01 -2.0054407861434758e+00 -5.5539086961084683e-01 8.9182669578022566e-01 2.6492297253172503e+00 + + + 5.1373541849242740e-02 2.5312983104460396e+00 -9.1276859073759216e+00 -6.0041477118315811e-01 4.4649449310393097e+00 5.1813568540466381e-01 7.5096056930460742e-01 + + + 2.0932016243092981e+00 -3.2452403852091161e-01 -3.0287002708624851e+00 -2.9435237868795006e+00 -2.4708373741958232e+00 -1.4834889562306977e+00 -5.6644366189909390e-01 + + + + + -6.4659259323887919e-01 + + + -1.4734579381772510e+00 + + + 1.2714410406766994e+00 + + + -1.0037354125666649e+00 + + + -8.3655432403244134e-01 + + + -3.5073967564578412e-01 + + + -7.5346418078087074e-01 + + + 1.9326707139761903e+00 + + + + + + + + diff --git a/CaloFilters/data/FlatPhoton_ReLU_noET.weights.xml b/CaloFilters/data/FlatPhoton_ReLU_noET.weights.xml new file mode 100644 index 0000000000..ce01fb76a4 --- /dev/null +++ b/CaloFilters/data/FlatPhoton_ReLU_noET.weights.xml @@ -0,0 +1,179 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 3.3107537831762990e-01 4.2145188861701097e+00 6.3213385284422241e+00 7.7975212511244896e-01 -9.4225123849402568e+00 -2.3576014339957445e+00 -4.5104821547667218e-01 + + + 2.8540299771303412e+00 -8.7000105543925155e-01 1.9690888222373553e+00 -3.0748114635564998e+00 -2.5000259966076337e+00 -1.5267973356496654e+00 4.1964598631512007e-01 + + + -1.7033379997744644e-01 -1.9588023151334316e-01 6.2972025122567332e-01 2.8359854922619032e+00 -9.9715012768932620e-01 -4.0144221992167323e-01 9.0980589348653618e-01 + + + -4.8205901909475518e-01 4.0289119044052457e-01 1.2231746580498558e-01 -1.7320561500389959e+00 -1.0319648143270015e+00 3.7468112296828471e-01 -3.0992022337672416e+00 + + + -8.1931531037141980e-01 3.3478542902645164e-01 1.7208670974275728e+00 9.7905100679002066e-01 -2.0719375157242998e+00 5.1528157825641341e-01 2.3767376627325398e-01 + + + 1.8696275596837251e+00 -6.3659175845593619e-01 2.0504382162649351e+00 -2.1472523920694324e+00 -6.4310371230961216e-01 2.9905645649044934e-02 8.1428096327316002e-01 + + + -2.5316571145409583e-01 -9.3459404974357729e-01 -3.5656937764344536e+00 -7.4340353271177073e-01 3.5373090443245849e+00 9.7045098758637727e-01 8.0482003539710426e-02 + + + 5.0614359375080853e-01 1.4542084170435916e+00 -7.4495236437054940e-02 -2.7126942451480867e+00 -7.2476166826788091e-01 -2.8426017144950264e+00 -9.3620659546555507e-01 + + + + + -6.2233684573909631e-01 + + + -1.4572878644552416e+00 + + + 1.2809016912853439e+00 + + + -1.5243470766136364e+00 + + + -1.0065018451303442e+00 + + + -6.2418174896974665e-01 + + + -1.6612861249436965e+00 + + + 2.1004564023688328e+00 + + + + + + + + diff --git a/CaloFilters/src/CaloNNEval_module.cc b/CaloFilters/src/CaloNNEval_module.cc new file mode 100644 index 0000000000..e3930e52c1 --- /dev/null +++ b/CaloFilters/src/CaloNNEval_module.cc @@ -0,0 +1,133 @@ +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "cetlib_except/exception.h" +#include "fhiclcpp/types/Atom.h" + +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CaloCluster/inc/ClusterUtils.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/Mu2eUtilities/inc/MVATools.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/TriggerInfo.hh" +#include "Offline/RecoDataProducts/inc/MVAResult.hh" + +#include +#include + + +namespace mu2e { + + class CaloNNEval : public art::EDProducer + { + public: + struct Config + { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom caloClusterCollection { Name("caloClusterCollection"), Comment("Calo cluster collection name") }; + fhicl::Table caloBkgMVA { Name("caloBkgMVA"), Comment("MVA configuration") }; + fhicl::Atom minEtoTest { Name("minEtoTest"), Comment("Minimum cluster energy to run the MVA") }; + fhicl::Atom minRtoTest { Name("minRtoTest"), Comment("Minimum cluster radius to run the MVA") }; + fhicl::Atom minTtoTest { Name("minTtoTest"), Comment("Minimum cluster time to run the MVA") }; + fhicl::Atom maxEtoTest { Name("maxEtoTest"), Comment("Maximum cluster energy to run the MVA") }; + fhicl::Atom maxRtoTest { Name("maxRtoTest"), Comment("Maximum cluster radius to run the MVA") }; + fhicl::Atom maxTtoTest { Name("maxTtoTest"), Comment("Maximum cluster time to run the MVA") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; + }; + + explicit CaloNNEval(const art::EDProducer::Table& config) : + EDProducer{config}, + caloClusterToken_{consumes(config().caloClusterCollection())}, + caloBkgMVA_ (config().caloBkgMVA()), + minEtoTest_ (std::max(config().minEtoTest(),0.001f)), + minRtoTest_ (config().minRtoTest()), + minTtoTest_ (config().minTtoTest()), + maxEtoTest_ (config().maxEtoTest()), + maxRtoTest_ (config().maxRtoTest()), + maxTtoTest_ (config().maxTtoTest()), + diagLevel_ (config().diagLevel()) + { + produces(); + } + + void produce(art::Event& event) override; + void beginJob() override; + + + private: + art::ProductToken caloClusterToken_; + MVATools caloBkgMVA_; + float minEtoTest_; + float minRtoTest_; + float minTtoTest_; + float maxEtoTest_; + float maxRtoTest_; + float maxTtoTest_; + int diagLevel_; + + void evalClusters(const art::Handle& caloClustersHandle, MVAResultCollection& MVAresults); + }; + + + void CaloNNEval::beginJob() + { + caloBkgMVA_.initMVA(); + } + + + + void CaloNNEval::produce(art::Event& event) + { + art::Handle caloClustersHandle = event.getHandle(caloClusterToken_); + auto MVAresults = std::make_unique(); + evalClusters(caloClustersHandle, *MVAresults); + event.put(std::move(MVAresults)); + } + + + //---------------------------------------------------------------------------------------------------------- + void CaloNNEval::evalClusters(const art::Handle& caloClustersHandle, MVAResultCollection& MVAresults) + { + if (!caloClustersHandle.isValid()) return; + + const Calorimeter& cal = *(GeomHandle()); + const CaloClusterCollection& caloClusters(*caloClustersHandle); + + for (auto clusterIt=caloClusters.begin(); clusterIt != caloClusters.end();++clusterIt){ + + float cluE = clusterIt->energyDep(); + float cluR = clusterIt->cog3Vector().perp(); + float cluT = clusterIt->time(); + if (cluE < minEtoTest_ || cluE > maxEtoTest_ || cluR < minRtoTest_ || cluR > maxRtoTest_ || + cluT < minTtoTest_ || cluT > maxTtoTest_) + { + MVAresults.push_back(MVAResult(-1.0)); + if (diagLevel_>0) std::cout<<"CaloNNEval skip cluster with e="<energyDep() + <<" r="<< clusterIt->cog3Vector().perp()<<" t="<time()<<"\n"; + continue; + } + + ClusterUtils cluUtil(cal,*clusterIt); + + std::vector mvavars(7,0.0); + mvavars[0] = clusterIt->cog3Vector().perp(); + mvavars[1] = clusterIt->size(); + mvavars[2] = cluUtil.e1()/clusterIt->energyDep(); + mvavars[3] = cluUtil.e2()/clusterIt->energyDep(); + mvavars[4] = cluUtil.e9()/clusterIt->energyDep(); + mvavars[5] = cluUtil.e25()/clusterIt->energyDep(); + mvavars[6] = cluUtil.secondMoment(); + + float mvaout = caloBkgMVA_.evalMVA(mvavars); + MVAresults.push_back(MVAResult(mvaout)); + if (diagLevel_>0) std::cout<<"CaloNNEval cluster with e="<energyDep() + <<" and r="<< clusterIt->cog3Vector().perp() + <<"has mvaout="< -#include - - namespace mu2e { @@ -25,18 +17,28 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - fhicl::Atom caloClusterCollection { Name("caloClusterCollection"), Comment("Calo cluster collection name") }; - fhicl::Table caloBkgMVA { Name("caloBkgMVA"), Comment("MVA Configuration") }; - fhicl::Atom minEtoTest { Name("minEtoTest"), Comment("Minimum Energy to run the MVA") }; - fhicl::Atom minMVAScore { Name("minMVAScore"), Comment("MVA cut for signal") }; - fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; + fhicl::Atom caloClusterCollection { Name("caloClusterCollection"), Comment("Calo cluster collection name") }; + fhicl::Atom caloMVACollection { Name("caloMVACollection"), Comment("Calo MVA cluster collection name") }; + fhicl::Atom minEtoTest { Name("minEtoTest"), Comment("Minimum cluster energy to run the MVA") }; + fhicl::Atom minRtoTest { Name("minRtoTest"), Comment("Minimum cluster radius to run the MVA") }; + fhicl::Atom minTtoTest { Name("minTtoTest"), Comment("Minimum cluster time to run the MVA") }; + fhicl::Atom maxEtoTest { Name("maxEtoTest"), Comment("Maximum cluster energy to run the MVA") }; + fhicl::Atom maxRtoTest { Name("maxRtoTest"), Comment("Maximum cluster radius to run the MVA") }; + fhicl::Atom maxTtoTest { Name("maxTtoTest"), Comment("Maximum cluster time to run the MVA") }; + fhicl::Atom minMVAScore { Name("minMVAScore"), Comment("MVA cut for signal") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; }; explicit FilterEcalNNTrigger(const art::EDFilter::Table& config) : EDFilter{config}, caloClusterToken_{consumes(config().caloClusterCollection())}, - caloBkgMVA_ (config().caloBkgMVA()), + caloMVAToken_ {consumes(config().caloMVACollection())}, minEtoTest_ (config().minEtoTest()), + minRtoTest_ (config().minRtoTest()), + minTtoTest_ (config().minTtoTest()), + maxEtoTest_ (config().maxEtoTest()), + maxRtoTest_ (config().maxRtoTest()), + maxTtoTest_ (config().maxTtoTest()), minMVAScore_ (config().minMVAScore()), diagLevel_ (config().diagLevel()) { @@ -44,78 +46,69 @@ namespace mu2e { } bool filter(art::Event& event) override; - void beginJob() override; private: art::ProductToken caloClusterToken_; - MVATools caloBkgMVA_; - float minEtoTest_; - float minMVAScore_; - int diagLevel_; - - bool filterClusters(const art::Handle& caloClustersHandle, TriggerInfo& trigInfo); + art::ProductToken caloMVAToken_; + float minEtoTest_; + float minRtoTest_; + float minTtoTest_; + float maxEtoTest_; + float maxRtoTest_; + float maxTtoTest_; + float minMVAScore_; + int diagLevel_; + + bool filterClusters(const art::Handle& caloClustersHandle, + const art::Handle& caloMVAHandle, + TriggerInfo& trigInfo); }; - void FilterEcalNNTrigger::beginJob() - { - caloBkgMVA_.initMVA(); - } - - - bool FilterEcalNNTrigger::filter(art::Event& event) { - art::Handle caloClustersHandle = event.getHandle(caloClusterToken_); + art::Handle caloClustersHandle = event.getHandle(caloClusterToken_); + art::Handle caloMVAHandle = event.getHandle(caloMVAToken_); - auto trigInfo = std::make_unique(); - bool retval = filterClusters(caloClustersHandle, *trigInfo); - event.put(std::move(trigInfo)); + auto trigInfo = std::make_unique(); + bool retval = filterClusters(caloClustersHandle, caloMVAHandle, *trigInfo); + event.put(std::move(trigInfo)); - return retval; + return retval; } //---------------------------------------------------------------------------------------------------------- - bool FilterEcalNNTrigger::filterClusters(const art::Handle& caloClustersHandle, TriggerInfo& trigInfo) + bool FilterEcalNNTrigger::filterClusters(const art::Handle& caloClustersHandle, + const art::Handle& caloMVAHandle, + TriggerInfo& trigInfo) { - const Calorimeter& cal = *(GeomHandle()); - const CaloClusterCollection& caloClusters(*caloClustersHandle); - - bool select(false); - for (auto clusterIt=caloClusters.begin(); clusterIt != caloClusters.end();++clusterIt) - { - if (clusterIt->energyDep() < minEtoTest_) continue; - - const auto& hits = clusterIt->caloHitsPtrVector(); - const auto& neighborsId = cal.crystal(hits[0]->crystalID()).neighbors(); - const auto& nneighborsId = cal.crystal(hits[0]->crystalID()).nextNeighbors(); - - double e9(hits[0]->energyDep()),e25(hits[0]->energyDep()); - for (auto hit : hits) - { - if (std::find(neighborsId.begin(), neighborsId.end(), hit->crystalID()) != neighborsId.end()) {e9 += hit->energyDep();e25 += hit->energyDep();} - if (std::find(nneighborsId.begin(), nneighborsId.end(), hit->crystalID()) != nneighborsId.end()) {e25 += hit->energyDep();} - } - - std::vector mvavars(8,0.0); - mvavars[0] = clusterIt->energyDep(); - mvavars[1] = clusterIt->cog3Vector().perp(); - mvavars[2] = clusterIt->size(); - mvavars[3] = hits[0]->energyDep(); - mvavars[4] = (hits.size()>1) ? hits[0]->energyDep() + hits[1]->energyDep() : hits[0]->energyDep(); - mvavars[5] = e9; - mvavars[6] = e25; - mvavars[7] = clusterIt->diskID(); - - float mvaout = caloBkgMVA_.evalMVA(mvavars); - if (mvaout < minMVAScore_) continue; - - select = true; - size_t index = std::distance(caloClusters.begin(),clusterIt); - trigInfo._caloClusters.push_back(art::Ptr(caloClustersHandle,index)); - } + if (!caloClustersHandle.isValid() || !caloMVAHandle.isValid()) return false; + + const auto& caloClusters(*caloClustersHandle); + const auto& caloMVAs(*caloMVAHandle); + + if (caloClusters.size() != caloMVAs.size()){ + throw cet::exception("FILTER")<< "FilterEcalNNTrigger: Clusters and MVA collection sizes incpmpatible\n"; + } + + bool select(false); + for (size_t index=0;index maxEtoTest_) continue; + if (cluR < minRtoTest_ || cluR > maxRtoTest_) continue; + if (cluT < minTtoTest_ || cluT > maxTtoTest_) continue; + if (mvaout._value < minMVAScore_) continue; + + select = true; + trigInfo._caloClusters.push_back(art::Ptr(caloClustersHandle,index)); + } return select; } diff --git a/CaloReco/fcl/prolog.fcl b/CaloReco/fcl/prolog.fcl index 2e4bf67905..d075acb7db 100644 --- a/CaloReco/fcl/prolog.fcl +++ b/CaloReco/fcl/prolog.fcl @@ -34,7 +34,8 @@ CaloReco : { } } -CaloReco : { @table::CaloReco +CaloReco : { + @table::CaloReco CaloRecoDigiMaker : { @@ -58,15 +59,35 @@ CaloReco : { @table::CaloReco time4Merge : 6.0 #ns diagLevel : 0 } + + CaloHitFastMaker : { + module_type : "CaloHitMakerFast" + ProtonBunchTimeTag : "EWMProducer" + caloDigiCollection : "CaloDigiMaker" + caphriEDepMax : 100 + caphriEDepMin : 1 + deltaTPulses : 11 + digiSampling : 5 + nPEperMeV : 30 + nSigmaNoise : 4 + noiseLevelMeV : 5.5e-1 + ADCToMeV : @local::ADCToMeV + diagLevel : 0 + } } -CaloReco : { @table::CaloReco +CaloReco : { + @table::CaloReco + producers : { CaloRecoDigiMaker : { @table::CaloReco.CaloRecoDigiMaker} CaloHitMaker : { @table::CaloReco.CaloHitMaker} + CaloHitFastMaker : { @table::CaloReco.CaloHitFastMaker} } - Reco : [ CaloRecoDigiMaker, CaloHitMaker ] + Reco : [ CaloRecoDigiMaker, CaloHitMaker, CaloHitFastMaker ] } + + END_PROLOG diff --git a/CaloReco/src/CaloHitMakerFast_module.cc b/CaloReco/src/CaloHitMakerFast_module.cc index 7f4ff2e3f3..1d0b9e044f 100644 --- a/CaloReco/src/CaloHitMakerFast_module.cc +++ b/CaloReco/src/CaloHitMakerFast_module.cc @@ -15,6 +15,8 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/GeometryService.hh" +#include + namespace { @@ -41,11 +43,11 @@ namespace mu2e { using Comment = fhicl::Comment; fhicl::Atom caloDigiCollection { Name("caloDigiCollection"), Comment("CaloDigi collection name") }; fhicl::Atom pbttoken { Name("ProtonBunchTimeTag"), Comment("ProtonBunchTime producer")}; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom deltaTPulses { Name("deltaTPulses"), Comment("Maximum time difference between two signals") }; - fhicl::Atom nPEperMeV { Name("nPEperMeV"), Comment("number of photo-electrons per MeV") }; - fhicl::Atom noiseLevelMeV { Name("noiseLevelMeV"), Comment("Noise level in MeV") }; - fhicl::Atom nSigmaNoise { Name("nSigmaNoise"), Comment("Maxnumber of sigma Noise to combine digi") }; + fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; + fhicl::Atom deltaTPulses { Name("deltaTPulses"), Comment("Maximum time difference between two signals") }; + fhicl::Atom nPEperMeV { Name("nPEperMeV"), Comment("number of photo-electrons per MeV") }; + fhicl::Atom noiseLevelMeV { Name("noiseLevelMeV"), Comment("Noise level in MeV") }; + fhicl::Atom nSigmaNoise { Name("nSigmaNoise"), Comment("Maxnumber of sigma Noise to combine digi") }; fhicl::Atom caphriEDepMax { Name("caphriEDepMax"), Comment("Maximum CAPHRI hit energy in MeV")}; fhicl::Atom caphriEDepMin { Name("caphriEDepMin"), Comment("Minimum CAPHRI hit energy in MeV")}; fhicl::Atom ADCToMeV { Name("ADCToMeV"), Comment("ADC to MeV conversion factor") }; @@ -75,22 +77,22 @@ namespace mu2e { private: - using pulseMapType = std::unordered_map>; + using pulseMapType = std::map>; void extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intCalo, double pbtOffset); void addPulse(pulseMapType& pulseMap, unsigned crystalID, float time, float eDep); - art::ProductToken caloDigisToken_; - const art::ProductToken pbttoken_; - double digiSampling_; - double deltaTPulses_; - double nPEperMeV_; - double noise2_; - double nSigmaNoise_; - float caphriEDepMax_; - float caphriEDepMin_; - double ADCToMeV_; - int diagLevel_; + art::ProductToken caloDigisToken_; + const art::ProductToken pbttoken_; + float digiSampling_; + float deltaTPulses_; + float nPEperMeV_; + float noise2_; + float nSigmaNoise_; + float caphriEDepMax_; + float caphriEDepMin_; + float ADCToMeV_; + int diagLevel_; }; @@ -124,18 +126,30 @@ namespace mu2e { //-------------------------------------------------------------------------------------------------------------- void CaloHitMakerFast::extractHits(const CaloDigiCollection& caloDigis, CaloHitCollection& caloHitsColl, CaloHitCollection& caphriHitsColl, IntensityInfoCalo& intInfo, double pbtOffset) { + std::vector digis(caloDigis.size()); + std::iota(digis.begin(),digis.end(),0); + + auto functorTime = [&caloDigis](auto a, auto b) { + auto dSi = caloDigis[a].SiPMID() - caloDigis[b].SiPMID(); + if (dSi==0) return caloDigis[a].t0() < caloDigis[b].t0(); + return dSi<0; + }; + std::sort(digis.begin(),digis.end(),functorTime); + pulseMapType pulseMap; - for (const auto& caloDigi : caloDigis) + for (const auto index : digis) { + const auto& caloDigi = caloDigis[index]; int crystalID = CaloSiPMId(caloDigi.SiPMID()).crystal().id(); size_t nSamPed = caloDigi.peakpos() > 3 ? 4 : std::max(caloDigi.peakpos()-1, 1); double baseline(0); for (size_t i=0; i 2) std::cout<<"[CaloHitMakerFast] extracted Digi with crystalID="<