800 me0 segment time pr

RecoLocalMuon/GEMRecHit/interface/ME0RecHitBaseAlgo.h

@@ -48,6 +48,7 @@ class ME0RecHitBaseAlgo {
                              LocalPoint& Point,
                              LocalError& error) const = 0;
 
+  bool recOnlyMuons;
 };
 #endif
 

RecoLocalMuon/GEMRecHit/python/me0RecHits_cfi.py

@@ -1,7 +1,9 @@
 import FWCore.ParameterSet.Config as cms
 
 me0RecHits = cms.EDProducer("ME0RecHitProducer",
-    recAlgoConfig = cms.PSet(),
+    recAlgoConfig = cms.PSet(
+        recOnlyMuons = cms.bool(False),
+    ),
     recAlgo = cms.string('ME0RecHitStandardAlgo'),
     me0DigiLabel = cms.InputTag("simMuonME0Digis"),
 )

RecoLocalMuon/GEMRecHit/src/ME0RecHitBaseAlgo.cc

@@ -13,8 +13,8 @@
 #include "DataFormats/GEMDigi/interface/ME0DigiPreRecoCollection.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
-
 ME0RecHitBaseAlgo::ME0RecHitBaseAlgo(const edm::ParameterSet& config) {
+  recOnlyMuons = config.getParameter<bool>("recOnlyMuons");
 }
 
 ME0RecHitBaseAlgo::~ME0RecHitBaseAlgo(){}
@@ -33,12 +33,9 @@ const ME0DigiPreRecoCollection::Range& digiRange){
     // Call the compute method
     bool OK = this->compute(*digi, point, tmpErr);
     if (!OK) continue;
-
-    if (std::abs(digi->pdgid()) == 13) {
-       ME0RecHit* recHit = new ME0RecHit(me0Id,digi->tof(),point,tmpErr);
-       result.push_back(recHit);
-    }
-
+    if (recOnlyMuons && std::abs(digi->pdgid()) != 13)  continue;
+    ME0RecHit* recHit = new ME0RecHit(me0Id,digi->tof(),point,tmpErr);
+    result.push_back(recHit);
   }
   return result;
 }

RecoLocalMuon/GEMSegment/plugins/ME0SegAlgoMM.cc

@@ -32,7 +32,13 @@ ME0SegAlgoMM::ME0SegAlgoMM(const edm::ParameterSet& ps) : ME0SegmentAlgorithm(ps
   preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
   dPhiChainBoxMax           = ps.getParameter<double>("dPhiChainBoxMax");
   dEtaChainBoxMax           = ps.getParameter<double>("dEtaChainBoxMax");
+  dTimeChainBoxMax          = ps.getParameter<double>("dTimeChainBoxMax");
   maxRecHitsInCluster       = ps.getParameter<int>("maxRecHitsInCluster");
+
+  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::ctor] Parameters to build segments :: "
+				   << "preClustering = "<<preClustering<<" preClustering_useChaining = "<<preClustering_useChaining
+				   <<" dPhiChainBoxMax = "<<dPhiChainBoxMax<<" dEtaChainBoxMax = "<<dEtaChainBoxMax<<" dTimeChainBoxMax = "<<dTimeChainBoxMax
+				   <<" minHitsPerSegment = "<<minHitsPerSegment<<" maxRecHitsInCluster = "<<maxRecHitsInCluster;
 }
 
 /* Destructor
@@ -45,10 +51,18 @@ ME0SegAlgoMM::~ME0SegAlgoMM() {
 std::vector<ME0Segment> ME0SegAlgoMM::run(const ME0Ensemble& ensemble, const EnsembleHitContainer& rechits) {
 
   theEnsemble = ensemble;
-
-  ME0DetId enId((theEnsemble.first)->id());
-  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::run] build segments in chamber " << enId;
-
+
+  #ifdef EDM_ML_DEBUG // have lines below only compiled when in debug mode
+  ME0DetId chId((theEnsamble.first)->id());  
+  edm::LogVerbatim("GEMSegAlgoMM") << "[ME0SegAlgoMM::run] build segments in chamber " << chId << " which contains "<<rechits.size()<<" rechits";
+  for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
+    auto me0id = (*rh)->me0Id();
+    auto rhLP = (*rh)->localPosition();
+    edm::LogVerbatim("ME0SegAlgoMM") << "[RecHit :: Loc x = "<<std::showpos<<std::setw(9)<<rhLP.x()<<" Loc y = "<<std::showpos<<std::setw(9)<<rhLP.y()
+                                     <<" Time = "<<std::showpos<<(*rh)->tof()<<" -- "<<me0id.rawId()<<" = "<<me0id<<" ]";
+  }
+  #endif
+
   // pre-cluster rechits and loop over all sub clusters separately
   std::vector<ME0Segment>          segments_temp;
   std::vector<ME0Segment>          segments;
@@ -256,20 +270,22 @@ ME0SegAlgoMM::chainHits(const EnsembleHitContainer& rechits) {
 
 bool ME0SegAlgoMM::isGoodToMerge(const EnsembleHitContainer& newChain, const EnsembleHitContainer& oldChain) {
 
-  std::vector<float> phi_new, eta_new, phi_old, eta_old;
+  std::vector<float> phi_new, eta_new, time_new, phi_old, eta_old, time_old;
   std::vector<int> layer_new, layer_old;
 
   for(size_t iRH_new = 0;iRH_new<newChain.size();++iRH_new){
     GlobalPoint pos_new = theEnsemble.first->toGlobal(newChain[iRH_new]->localPosition());
     layer_new.push_back(newChain[iRH_new]->me0Id().layer());
     phi_new.push_back(pos_new.phi());
     eta_new.push_back(pos_new.eta());
+    time_new.push_back(newChain[iRH_new]->tof());
   }  
   for(size_t iRH_old = 0;iRH_old<oldChain.size();++iRH_old){
     GlobalPoint pos_old = theEnsemble.first->toGlobal(oldChain[iRH_old]->localPosition());
     layer_old.push_back(oldChain[iRH_old]->me0Id().layer());
     phi_old.push_back(pos_old.phi());
     eta_old.push_back(pos_old.eta());
+    time_old.push_back(oldChain[iRH_old]->tof());
   }
 
   for(size_t jRH_new = 0; jRH_new<phi_new.size(); ++jRH_new){
@@ -288,8 +304,11 @@ bool ME0SegAlgoMM::isGoodToMerge(const EnsembleHitContainer& newChain, const Ens
       bool etaRequirementOK = fabs(eta_new[jRH_new]-eta_old[jRH_old]) < dEtaChainBoxMax;
       // and the difference in layer index should be < (nlayers-1)
       bool layerRequirementOK = abs(layer_new[jRH_new]-layer_old[jRH_old]) < (theEnsemble.first->id().nlayers()-1);
+      // and they should have a time difference compatible with the hypothesis 
+      // that the rechits originate from the same particle, but were detected in different layers
+      bool timeRequirementOK = fabs(time_new[jRH_new] - time_old[jRH_old]) < dTimeChainBoxMax;
 
-      if(layerRequirementOK && phiRequirementOK && etaRequirementOK){
+      if(layerRequirementOK && phiRequirementOK && etaRequirementOK && timeRequirementOK){
         return true;
       }
     }
@@ -304,6 +323,16 @@ bool ME0SegAlgoMM::isGoodToMerge(const EnsembleHitContainer& newChain, const Ens
 std::vector<ME0Segment> ME0SegAlgoMM::buildSegments(const EnsembleHitContainer& rechits) {
   std::vector<ME0Segment> me0segs;
 
+  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] will now try to fit a ME0Segment from collection of "<<rechits.size()<<" ME0 RecHits";
+  #ifdef EDM_ML_DEBUG // have lines below only compiled when in debug mode 
+    for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
+      auto me0id = (*rh)->me0Id();
+      auto rhLP = (*rh)->localPosition();
+      edm::LogVerbatim("ME0SegAlgoMM") << "[RecHit :: Loc x = "<<std::showpos<<std::setw(9)<<rhLP.x()<<" Loc y = "<<std::showpos<<std::setw(9)<<rhLP.y()
+				       <<" Time = "<<std::showpos<<(*rh)->tof()<<" -- "<<me0id.rawId()<<" = "<<me0id<<" ]";
+      }
+  #endif
+
   proto_segment.clear();
   // select hits from the ensemble and sort it 
   for (auto rh=rechits.begin(); rh!=rechits.end();rh++){

RecoLocalMuon/GEMSegment/plugins/ME0SegAlgoMM.h

@@ -68,6 +68,7 @@ class ME0SegAlgoMM : public ME0SegmentAlgorithm {
   bool    preClustering_useChaining;
   double  dPhiChainBoxMax;
   double  dEtaChainBoxMax;
+  double  dTimeChainBoxMax;
   int     maxRecHitsInCluster;
 
  private:

RecoLocalMuon/GEMSegment/plugins/ME0SegmentBuilder.cc

@@ -63,7 +63,7 @@ void ME0SegmentBuilder::build(const ME0RecHitCollection* recHits, ME0SegmentColl
     }    
     ME0SegmentAlgorithm::ME0Ensemble ensemble(std::pair<const ME0EtaPartition*, std::map<uint32_t,const ME0EtaPartition*> >(firstlayer,ens));
 
-    // LogDebug("ME0SegmentBuilder") << "found " << me0RecHits.size() << " rechits in chamber " << *enIt;
+    LogDebug("ME0SegmentBuilder") << "found " << me0RecHits.size() << " rechits in chamber " /*<< *enIt */;
 
     // given the chamber select the appropriate algo... and run it
     std::vector<ME0Segment> segv = algo->run(ensemble, me0RecHits);

RecoLocalMuon/GEMSegment/python/me0Segments_cfi.py

@@ -12,6 +12,7 @@
         preClusteringUseChaining = cms.bool(True),
         dPhiChainBoxMax = cms.double(.02),
         dEtaChainBoxMax = cms.double(.05),
+        dTimeChainBoxMax = cms.double(1.50), # 1ns, +/- time to fly through 30cm thick ME0
         maxRecHitsInCluster = cms.int32(6)
     )
 )

SimMuon/GEMDigitizer/interface/ME0PreRecoGaussianModel.h

@@ -37,32 +37,16 @@ class ME0PreRecoGaussianModel: public ME0DigiPreRecoModel
   bool etaproj;
   bool digitizeOnlyMuons_;
   double averageEfficiency_;
-  bool doBkgNoise_;
-  bool simulateIntrinsicNoise_;
+  // bool simulateIntrinsicNoise_; // not implemented
+  // double averageNoiseRate_;     // not implemented
   bool simulateElectronBkg_;
+  bool simulateNeutralBkg_;
 
-  double averageNoiseRate_;
-  int bxwidth_;
   int minBunch_;
   int maxBunch_;
 
-  //params for the simple pol6 model of neutral bkg for ME0:
-  double ME0ModNeuBkgParam0;
-  double ME0ModNeuBkgParam1;
-  double ME0ModNeuBkgParam2;
-  double ME0ModNeuBkgParam3;
-  double ME0ModNeuBkgParam4;
-  double ME0ModNeuBkgParam5;
-  double ME0ModNeuBkgParam6;
-
-  double ME0ModElecBkgParam0;
-  double ME0ModElecBkgParam1;
-  double ME0ModElecBkgParam2;
-  double ME0ModElecBkgParam3;
-  double ME0ModElecBkgParam4;
-  double ME0ModElecBkgParam5;
-  double ME0ModElecBkgParam6;
-  double ME0ModElecBkgParam7;
+  // params for the simple pol6 model of neutral bkg for ME0:
+  std::vector<double> neuBkg, eleBkg;
 
 };
 #endif

SimMuon/GEMDigitizer/python/muonME0DigisPreReco_cfi.py

@@ -2,21 +2,19 @@
 
 # Module to create simulated ME0 Pre Reco digis.
 simMuonME0Digis = cms.EDProducer("ME0DigiPreRecoProducer",
-    mixLabel = cms.string('mix'),
     inputCollection = cms.string('g4SimHitsMuonME0Hits'),
     digiPreRecoModelString = cms.string('PreRecoGaussian'),
-    timeResolution = cms.double(0.010),      # [in ns] => for now at 10ps 
-    phiResolution = cms.double(0.05),        # [in cm] average resolution along local x in case of no correlation
-    etaResolution = cms.double(1.),          # [in cm] average resolution along local y in case of no correlation
-    useCorrelation  = cms.bool(False),
-    useEtaProjectiveGEO  = cms.bool(False),
+    timeResolution = cms.double(0.001), # in ns
+    phiResolution = cms.double(0.05),   # in cm average resolution along local x in case of no correlation
+    etaResolution = cms.double(1.),     # in cm average resolution along local y in case of no correlation
+    useCorrelation = cms.bool(False),
+    useEtaProjectiveGEO = cms.bool(False),
     averageEfficiency = cms.double(0.98),
-    doBkgNoise = cms.bool(False),            # False => No background noise simulation
     digitizeOnlyMuons = cms.bool(False),
-    simulateIntrinsicNoise = cms.bool(False),
-    simulateElectronBkg = cms.bool(False),   # True => will simulate electron background
-    averageNoiseRate = cms.double(0.001),    # simulation of intrinsic noise
-    bxwidth = cms.int32(25),
-    minBunch = cms.int32(-5),                # in bx-units (x 25 ns)
-    maxBunch = cms.int32(3)
+    # simulateIntrinsicNoise = cms.bool(False), # intrinsic noise --> not implemented
+    # averageNoiseRate = cms.double(0.001),     # intrinsic noise --> not implemented
+    simulateElectronBkg = cms.bool(True),       # True - will simulate electron background
+    simulateNeutralBkg  = cms.bool(True),       # True - will simulate neutral (n+g)  background
+    minBunch = cms.int32(-5),                   # [x 25 ns], forms the readout window together with maxBunch,
+    maxBunch = cms.int32(3),                    # we should think of shrinking this window ...
 )