#include <GEMSimpleModel.h>

Inheritance diagram for GEMSimpleModel:

Public Member Functions
	GEMSimpleModel (const edm::ParameterSet &)

int	getSimHitBx (const PSimHit , CLHEP::HepRandomEngine )

void	setup () override

std::vector< std::pair< int, int > >	simulateClustering (const GEMEtaPartition , const PSimHit , const int, CLHEP::HepRandomEngine *) override

void	simulateNoise (const GEMEtaPartition , CLHEP::HepRandomEngine ) override

void	simulateSignal (const GEMEtaPartition , const edm::PSimHitContainer &, CLHEP::HepRandomEngine ) override

	~GEMSimpleModel ()

Public Member Functions inherited from GEMDigiModel
void	fillDigis (int rollDetId, GEMDigiCollection &)

const GEMDigiSimLinks &	gemDigiSimLinks () const

const GEMGeometry *	getGeometry ()

void	setGeometry (const GEMGeometry *geom)

const StripDigiSimLinks &	stripDigiSimLinks () const

virtual	~GEMDigiModel ()

Private Attributes
double	averageEfficiency_

double	averageNoiseRate_

double	averageShapingTime_

int	bxwidth_

bool	digitizeOnlyMuons_

bool	doBkgNoise_

bool	doNoiseCLS_

bool	fixedRollRadius_

double	GE11ElecBkgParam0_

double	GE11ElecBkgParam1_

double	GE11ElecBkgParam2_

double	GE11ModNeuBkgParam0_

double	GE11ModNeuBkgParam1_

double	GE11ModNeuBkgParam2_

double	GE21ElecBkgParam0_

double	GE21ElecBkgParam1_

double	GE21ElecBkgParam2_

double	GE21ModNeuBkgParam0_

double	GE21ModNeuBkgParam1_

double	GE21ModNeuBkgParam2_

double	instLumi_

int	maxBunch_

int	minBunch_

double	rateFact_

const double	referenceInstLumi_

double	signalPropagationSpeed_

bool	simulateElectronBkg_

bool	simulateIntrinsicNoise_

double	timeJitter_

double	timeResolution_

Additional Inherited Members
Public Types inherited from GEMDigiModel
typedef edm::DetSet< GEMDigiSimLink >	GEMDigiSimLinks

typedef edm::DetSet< StripDigiSimLink >	StripDigiSimLinks

Protected Types inherited from GEMDigiModel
typedef std::multimap< std::pair< unsigned int, int >, const PSimHit *, std::less< std::pair< unsigned int, int > > >	DetectorHitMap

Protected Member Functions inherited from GEMDigiModel
void	addLinks (unsigned int strip, int bx)
	creates links from Digi to SimTrack More...

void	addLinksWithPartId (unsigned int strip, int bx)

	GEMDigiModel (const edm::ParameterSet &)

Protected Attributes inherited from GEMDigiModel
DetectorHitMap	detectorHitMap_

const GEMGeometry *	geometry_

StripDigiSimLinks	stripDigiSimLinks_

std::set< std::pair< int, int > >	strips_

GEMDigiSimLinks	theGemDigiSimLinks_

Detailed Description

Class for the GEM strip response simulation based on a very simple model

Author: Sven Dildick by Roumyana Hadjiiska

Definition at line 22 of file GEMSimpleModel.h.

Constructor & Destructor Documentation

GEMSimpleModel::GEMSimpleModel ( const edm::ParameterSet & config )

Definition at line 13 of file GEMSimpleModel.cc.

                                                             :
 GEMDigiModel(config)
 , averageEfficiency_(config.getParameter<double> ("averageEfficiency"))
 , averageShapingTime_(config.getParameter<double> ("averageShapingTime"))
 , timeResolution_(config.getParameter<double> ("timeResolution"))
 , timeJitter_(config.getParameter<double> ("timeJitter"))
 , averageNoiseRate_(config.getParameter<double> ("averageNoiseRate"))
 , signalPropagationSpeed_(config.getParameter<double> ("signalPropagationSpeed"))
 , bxwidth_(config.getParameter<int> ("bxwidth"))
 , minBunch_(config.getParameter<int> ("minBunch"))
 , maxBunch_(config.getParameter<int> ("maxBunch"))
 , digitizeOnlyMuons_(config.getParameter<bool> ("digitizeOnlyMuons"))
 , doBkgNoise_(config.getParameter<bool> ("doBkgNoise"))
 , doNoiseCLS_(config.getParameter<bool> ("doNoiseCLS"))
 , fixedRollRadius_(config.getParameter<bool> ("fixedRollRadius"))
 , simulateElectronBkg_(config.getParameter<bool> ("simulateElectronBkg"))
 , instLumi_(config.getParameter<double>("instLumi"))
 , rateFact_(config.getParameter<double>("rateFact"))
 , referenceInstLumi_(config.getParameter<double>("referenceInstLumi"))
 , GE11ElecBkgParam0_(config.getParameter<double>("GE11ElecBkgParam0"))
 , GE11ElecBkgParam1_(config.getParameter<double>("GE11ElecBkgParam1"))
 , GE11ElecBkgParam2_(config.getParameter<double>("GE11ElecBkgParam2"))
 , GE21ElecBkgParam0_(config.getParameter<double>("GE21ElecBkgParam0"))
 , GE21ElecBkgParam1_(config.getParameter<double>("GE21ElecBkgParam1"))
 , GE21ElecBkgParam2_(config.getParameter<double>("GE21ElecBkgParam2"))
 , GE11ModNeuBkgParam0_(config.getParameter<double>("GE11ModNeuBkgParam0"))
 , GE11ModNeuBkgParam1_(config.getParameter<double>("GE11ModNeuBkgParam1"))
 , GE11ModNeuBkgParam2_(config.getParameter<double>("GE11ModNeuBkgParam2"))
 , GE21ModNeuBkgParam0_(config.getParameter<double>("GE11ModNeuBkgParam0"))
 , GE21ModNeuBkgParam1_(config.getParameter<double>("GE11ModNeuBkgParam1"))
 , GE21ModNeuBkgParam2_(config.getParameter<double>("GE11ModNeuBkgParam2"))
 {
 }

GEMSimpleModel::~GEMSimpleModel ( )

Definition at line 47 of file GEMSimpleModel.cc.

48 {

49 }

Member Function Documentation

int GEMSimpleModel::getSimHitBx	(	const PSimHit *	simhit,
		CLHEP::HepRandomEngine *	engine
	)

Definition at line 99 of file GEMSimpleModel.cc.

References averageShapingTime_, bxwidth_, gather_cfg::cout, debug, PSimHit::detUnitId(), MillePedeFileConverter_cfg::e, GEMGeometry::etaPartition(), Exception, GEMDigiModel::geometry_, hcalTTPDigis_cfi::id, PSimHit::localPosition(), signalPropagationSpeed_, mathSSE::sqrt(), timeJitter_, PSimHit::timeOfFlight(), timeResolution_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by simulateSignal().

 {
   int bx = -999;
   const LocalPoint simHitPos(simhit->localPosition());
   const float tof(simhit->timeOfFlight());
   // random Gaussian time correction due to electronics jitter
   float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
   const GEMDetId id(simhit->detUnitId());
   const GEMEtaPartition* roll(geometry_->etaPartition(id));
   if (!roll)
   {
     throw cms::Exception("Geometry")<< "GEMSimpleModel::getSimHitBx() - GEM simhit id does not match any GEM roll id: " << id << "\n";
     return 999;
   }
   if (roll->id().region() == 0)
   {
     throw cms::Exception("Geometry") << "GEMSimpleModel::getSimHitBx() - this GEM id is from barrel, which cannot happen: " << roll->id() << "\n";
   }
   const double cspeed = 299792458;   // signal propagation speed in vacuum in [m/s]
   const int nstrips = roll->nstrips();
   float middleStrip = nstrips/2.;
   const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
   const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
   double muRadius = sqrt(globMiddleRol.x()*globMiddleRol.x() + globMiddleRol.y()*globMiddleRol.y() +globMiddleRol.z()*globMiddleRol.z());
   double timeCalibrationOffset_ = (muRadius *1e+9)/(cspeed*1e+2); //[ns]
 
   const TrapezoidalStripTopology* top(dynamic_cast<const TrapezoidalStripTopology*> (&(roll->topology())));
   const float halfStripLength(0.5 * top->stripLength());
   const float distanceFromEdge(halfStripLength - simHitPos.y());
 
   // signal propagation speed in material in [cm/ns]
   double signalPropagationSpeedTrue = signalPropagationSpeed_ * cspeed * 1e-7;  // 1e+2 * 1e-9;
 
   // average time for the signal to propagate from the SimHit to the top of a strip
   const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
   // random Gaussian time correction due to the finite timing resolution of the detector
   float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
   const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
   float referenceTime = 0.;
   referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
   const float timeDifference(simhitTime - referenceTime);
   // assign the bunch crossing
   bx = static_cast<int> (std::round((timeDifference) / bxwidth_));
 
   // check time
   const bool debug(false);
   if (debug)
   {
     std::cout << "checktime " << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT =  " << simhitTime
         << "\trefT =  " << referenceTime << "\ttof = " << tof << "\tavePropT =  " << averagePropagationTime
         << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
   }
   return bx;
 }

void GEMSimpleModel::setup ( )

overridevirtual

Implements GEMDigiModel.

Definition at line 51 of file GEMSimpleModel.cc.

 {
   return;
 }

std::vector< std::pair< int, int > > GEMSimpleModel::simulateClustering	(	const GEMEtaPartition *	roll,
		const PSimHit *	simHit,
		const int	bx,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements GEMDigiModel.

Definition at line 269 of file GEMSimpleModel.cc.

References averageEfficiency_, GEMEtaPartition::centreOfStrip(), Topology::channel(), PSimHit::localPosition(), GEMEtaPartition::nstrips(), GEMEtaPartition::specificTopology(), GeomDet::toGlobal(), ecaldqm::topology(), and PV3DBase< T, PVType, FrameType >::x().

Referenced by simulateSignal().

 {
   const StripTopology& topology = roll->specificTopology(); // const LocalPoint& entry(simHit->entryPoint());
   const LocalPoint& hit_position(simHit->localPosition());
   const int nstrips(roll->nstrips());
   int centralStrip = 0;
   if (!(topology.channel(hit_position) + 1 > nstrips))
     centralStrip = topology.channel(hit_position) + 1;
   else
     centralStrip = topology.channel(hit_position);
   const GlobalPoint& pointSimHit = roll->toGlobal(hit_position);
   const GlobalPoint& pointDigiHit = roll->toGlobal(roll->centreOfStrip(centralStrip));
   double deltaX = pointSimHit.x() - pointDigiHit.x();
 
   // Add central digi to cluster vector
   std::vector < std::pair<int, int> > cluster_;
   cluster_.clear();
   cluster_.emplace_back(centralStrip, bx);
 
   //simulate cross talk
   int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
   if (clusterSize == 2)
   {
     if (deltaX <= 0)
     {
       if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
         cluster_.emplace_back(centralStrip - 1, bx);
     }
     else
     {
       if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
         cluster_.emplace_back(centralStrip + 1, bx);
     }
   }
   return cluster_;
 }

void GEMSimpleModel::simulateNoise	(	const GEMEtaPartition *	roll,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements GEMDigiModel.

Definition at line 154 of file GEMSimpleModel.cc.

References averageEfficiency_, averageNoiseRate_, bxwidth_, doBkgNoise_, doNoiseCLS_, MillePedeFileConverter_cfg::e, Exception, fixedRollRadius_, GE11ElecBkgParam0_, GE11ElecBkgParam1_, GE11ElecBkgParam2_, GE11ModNeuBkgParam0_, GE11ModNeuBkgParam1_, GE11ModNeuBkgParam2_, GE21ElecBkgParam0_, GE21ElecBkgParam1_, GE21ElecBkgParam2_, GE21ModNeuBkgParam0_, GE21ModNeuBkgParam1_, GE21ModNeuBkgParam2_, mps_fire::i, GEMEtaPartition::id(), instLumi_, gen::k, maxBunch_, minBunch_, GEMEtaPartition::nstrips(), GEMEtaPartition::pitch(), rateFact_, referenceInstLumi_, simulateElectronBkg_, simulateIntrinsicNoise_, GEMDigiModel::strips_, and GEMEtaPartition::topology().

 {
   if (!doBkgNoise_)
     return;
   const GEMDetId& gemId(roll->id());
   const int nstrips(roll->nstrips());
   double trArea(0.0);
   double trStripArea(0.0);
   if (gemId.region() == 0)
   {
     throw cms::Exception("Geometry")
         << "GEMSynchronizer::simulateNoise() - this GEM id is from barrel, which cannot happen.";
   }
   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
   const float striplength(top_->stripLength());
   trStripArea = (roll->pitch()) * striplength;
   trArea = trStripArea * nstrips;
   const int nBxing(maxBunch_ - minBunch_ + 1);
   const float rollRadius(fixedRollRadius_ ? top_->radius() : 
              top_->radius() + CLHEP::RandFlat::shoot(engine, -1.*top_->stripLength()/2., top_->stripLength()/2.));
 
   //calculate noise from model
   double averageNeutralNoiseRatePerRoll = 0.;
   double averageNoiseElectronRatePerRoll = 0.;
   double averageNoiseRatePerRoll = 0.;
   if (gemId.station() == 1)
   {
     //simulate neutral background for GE1/1
     averageNeutralNoiseRatePerRoll = (GE11ModNeuBkgParam0_
                       + GE11ModNeuBkgParam1_ * rollRadius
                       + GE11ModNeuBkgParam2_ * rollRadius * rollRadius);    //simulate electron background for GE1/1
     if (simulateElectronBkg_)
       averageNoiseElectronRatePerRoll = (GE11ElecBkgParam0_
                      + GE11ElecBkgParam1_ * rollRadius
                      + GE11ElecBkgParam2_ * rollRadius * rollRadius);
 
     // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
     averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
     averageNoiseRatePerRoll *= instLumi_*rateFact_*1.0/referenceInstLumi_;
   }
   if (gemId.station() == 2)
   {
     //simulate neutral background for GE2/1
     averageNeutralNoiseRatePerRoll = (GE21ModNeuBkgParam0_
                       + GE21ModNeuBkgParam1_ * rollRadius
                       + GE21ModNeuBkgParam2_ * rollRadius * rollRadius);
     //simulate electron background for GE2/1
     if (simulateElectronBkg_)
       averageNoiseElectronRatePerRoll = (GE21ElecBkgParam0_
                      + GE21ElecBkgParam1_ * rollRadius
                      + GE21ElecBkgParam2_ * rollRadius * rollRadius);
 
     // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
     averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
     averageNoiseRatePerRoll *= instLumi_*rateFact_*1.0/referenceInstLumi_;
   }
 
   //simulate intrinsic noise
   if(simulateIntrinsicNoise_)
   {
     const double aveIntrinsicNoisePerStrip(averageNoiseRate_ * nBxing * bxwidth_ * trStripArea * 1.0e-9);
     for(int j = 0; j < nstrips; ++j)
     {
       CLHEP::RandPoissonQ randPoissonQ(*engine, aveIntrinsicNoisePerStrip);
       const int n_intrHits(randPoissonQ.fire());
     
       for (int k = 0; k < n_intrHits; k++ )
     {
         const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
         std::pair<int, int> digi(k+1,time_hit);
         strips_.emplace(digi);
       }
     }
   }//end simulate intrinsic noise
 
   //simulate bkg contribution
   const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9);
   CLHEP::RandPoissonQ randPoissonQ(*engine, averageNoise);
   const int n_hits(randPoissonQ.fire());
   for (int i = 0; i < n_hits; ++i)
   {
     const int centralStrip(static_cast<int> (CLHEP::RandFlat::shoot(engine, 1, nstrips)));
     const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
     if (doNoiseCLS_)
     {
       std::vector < std::pair<int, int> > cluster_;
       cluster_.clear();
       cluster_.emplace_back(centralStrip, time_hit);
       int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
       if (clusterSize == 2)
       {
         if(CLHEP::RandFlat::shoot(engine) < 0.5)
         {
           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
             cluster_.emplace_back(centralStrip - 1, time_hit);
         }
         else
         {
           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
             cluster_.emplace_back(centralStrip + 1, time_hit);
         }
       }
       for (const auto& digi : cluster_)
       {
         strips_.emplace(digi);
       }
     } //end doNoiseCLS_
     else
     {
       strips_.emplace(centralStrip, time_hit);
     }
   }
   return;
 }

void GEMSimpleModel::simulateSignal	(	const GEMEtaPartition *	roll,
		const edm::PSimHitContainer &	simHits,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements GEMDigiModel.

Definition at line 56 of file GEMSimpleModel.cc.

References funct::abs(), averageEfficiency_, edm::DetSet< T >::clear(), GEMDigiModel::detectorHitMap_, digitizeOnlyMuons_, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, getSimHitBx(), GEMEtaPartition::id(), cmsBatch::log, DetId::rawId(), simulateClustering(), GEMDigiModel::stripDigiSimLinks_, GEMDigiModel::strips_, and GEMDigiModel::theGemDigiSimLinks_.

 {
   stripDigiSimLinks_.clear();
   detectorHitMap_.clear();
   stripDigiSimLinks_ = StripDigiSimLinks(roll->id().rawId());
   theGemDigiSimLinks_.clear();
   theGemDigiSimLinks_ = GEMDigiSimLinks(roll->id().rawId());
   bool digiMuon = false;
   bool digiElec = false;
   for (const auto& hit : simHits)
   {
     if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
       continue;
     double elecEff = 0.;
     double partMom = hit.pabs();
     double checkMuonEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
     double checkElecEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
     if (std::abs(hit.particleType()) == 13 && checkMuonEff < averageEfficiency_)
       digiMuon = true;
     if (std::abs(hit.particleType()) != 13) //consider all non muon particles with gem efficiency to electrons
     {
       if (partMom <= 1.95e-03)
         elecEff = 1.7e-05 * std::exp(2.1 * partMom * 1000.);
       if (partMom > 1.95e-03 && partMom < 10.e-03)
         elecEff = 1.34 * log(7.96e-01 * partMom * 1000. - 5.75e-01)
             / (1.34 + log(7.96e-01 * partMom * 1000. - 5.75e-01));
       if (partMom > 10.e-03)
         elecEff = 1.;
       if (checkElecEff < elecEff)
         digiElec = true;
     }
    if (!(digiMuon || digiElec))
       continue;
     const int bx(getSimHitBx(&hit, engine));
     const std::vector<std::pair<int, int> >& cluster(simulateClustering(roll, &hit, bx, engine));
     for (const auto & digi : cluster)
     {
       detectorHitMap_.emplace(digi,&hit);
       strips_.emplace(digi);
     }
   }
 }