#include <ME0SimpleModel.h>

Inheritance diagram for ME0SimpleModel:

Public Member Functions
int	getSimHitBx (const PSimHit , CLHEP::HepRandomEngine )

	ME0SimpleModel (const edm::ParameterSet &)

void	setup () override

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

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

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

	~ME0SimpleModel ()

Public Member Functions inherited from ME0DigiModel
void	fillDigis (int rollDetId, ME0DigiCollection &)

const ME0Geometry *	getGeometry ()

const ME0DigiSimLinks &	me0DigiSimLinks () const

void	setGeometry (const ME0Geometry *geom)

const StripDigiSimLinks &	stripDigiSimLinks () const

virtual	~ME0DigiModel ()

Private Attributes
double	averageEfficiency_

double	averageNoiseRate_

double	averageShapingTime_

int	bxwidth_

bool	digitizeOnlyMuons_

bool	doBkgNoise_

bool	doNoiseCLS_

bool	fixedRollRadius_

double	instLumi_

int	maxBunch_

double	ME0ElecBkgParam0_

double	ME0ElecBkgParam1_

double	ME0ElecBkgParam2_

double	ME0ElecBkgParam3_

double	ME0NeuBkgParam0_

double	ME0NeuBkgParam1_

double	ME0NeuBkgParam2_

double	ME0NeuBkgParam3_

int	minBunch_

double	rateFact_

double	referenceInstLumi_

double	signalPropagationSpeed_

bool	simulateElectronBkg_

bool	simulateIntrinsicNoise_

double	timeJitter_

double	timeResolution_

Additional Inherited Members
Public Types inherited from ME0DigiModel
typedef edm::DetSet< ME0DigiSimLink >	ME0DigiSimLinks

typedef edm::DetSet< StripDigiSimLink >	StripDigiSimLinks

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

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

void	addLinksWithPartId (unsigned int strip, int bx)

	ME0DigiModel (const edm::ParameterSet &)

Protected Attributes inherited from ME0DigiModel
DetectorHitMap	detectorHitMap_

const ME0Geometry *	geometry_

StripDigiSimLinks	stripDigiSimLinks_

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

ME0DigiSimLinks	theME0DigiSimLinks_

Detailed Description

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

Author: Sven Dildick by Roumyana Hadjiiska

Definition at line 24 of file ME0SimpleModel.h.

Constructor & Destructor Documentation

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

Definition at line 17 of file ME0SimpleModel.cc.

                                                             :
 ME0DigiModel(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"))
 , ME0ElecBkgParam0_(config.getParameter<double> ("ME0ElecBkgParam0"))
 , ME0ElecBkgParam1_(config.getParameter<double> ("ME0ElecBkgParam1"))
 , ME0ElecBkgParam2_(config.getParameter<double> ("ME0ElecBkgParam2"))
 , ME0ElecBkgParam3_(config.getParameter<double> ("ME0ElecBkgParam3"))
 , ME0NeuBkgParam0_(config.getParameter<double> ("ME0NeuBkgParam0"))
 , ME0NeuBkgParam1_(config.getParameter<double> ("ME0NeuBkgParam1"))
 , ME0NeuBkgParam2_(config.getParameter<double> ("ME0NeuBkgParam2"))
 , ME0NeuBkgParam3_(config.getParameter<double> ("ME0NeuBkgParam3"))
 {
 }

ME0SimpleModel::~ME0SimpleModel ( )

Definition at line 47 of file ME0SimpleModel.cc.

48 {

49 }

Member Function Documentation

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

Definition at line 100 of file ME0SimpleModel.cc.

References averageShapingTime_, bxwidth_, debug, PSimHit::detUnitId(), ME0Geometry::etaPartition(), Exception, ME0DigiModel::geometry_, hcalTTPDigis_cfi::id, PSimHit::localPosition(), LogDebug, 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 ME0DetId& id(simhit->detUnitId());
   const ME0EtaPartition* roll(geometry_->etaPartition(id));
   if (!roll)
   {
     throw cms::Exception("Geometry")<< "ME0SimpleModel::getSimHitBx() - ME0 simhit id does not match any ME0 roll id: " << id << "\n";
     return 999;
   }
   if (roll->id().region() == 0)
   {
     throw cms::Exception("Geometry") << "ME0SimpleModel::getSimHitBx() - this ME0 id is from barrel, which cannot happen: " << roll->id() << "\n";
   }
   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*CLHEP::ns*CLHEP::cm)/(CLHEP::c_light); //[cm/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_ * CLHEP::c_light/(CLHEP::ns*CLHEP::cm);
   // 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)
   {
     LogDebug("ME0SimpleModel")<<  "checktime " << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT =  " << simhitTime 
         << "\trefT =  " << referenceTime << "\ttof = " << tof << "\tavePropT =  " << averagePropagationTime
         << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
   }
   return bx;
 }

void ME0SimpleModel::setup ( )

overridevirtual

Implements ME0DigiModel.

Definition at line 51 of file ME0SimpleModel.cc.

 {
   return;
 }

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

overridevirtual

Implements ME0DigiModel.

Definition at line 249 of file ME0SimpleModel.cc.

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

Referenced by simulateSignal().

 {
   const StripTopology& topology = roll->specificTopology();
   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 ME0SimpleModel::simulateNoise	(	const ME0EtaPartition *	roll,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements ME0DigiModel.

Definition at line 152 of file ME0SimpleModel.cc.

References averageEfficiency_, averageNoiseRate_, bxwidth_, doBkgNoise_, doNoiseCLS_, MillePedeFileConverter_cfg::e, Exception, JetChargeProducer_cfi::exp, fixedRollRadius_, mps_fire::i, ME0EtaPartition::id(), instLumi_, gen::k, maxBunch_, ME0ElecBkgParam0_, ME0ElecBkgParam1_, ME0ElecBkgParam2_, ME0ElecBkgParam3_, ME0NeuBkgParam0_, ME0NeuBkgParam1_, ME0NeuBkgParam2_, ME0NeuBkgParam3_, minBunch_, ME0EtaPartition::nstrips(), ME0EtaPartition::pitch(), rateFact_, referenceInstLumi_, simulateElectronBkg_, simulateIntrinsicNoise_, mathSSE::sqrt(), ME0DigiModel::strips_, and ME0EtaPartition::topology().

 {
   if (!doBkgNoise_)
     return;
   const ME0DetId me0Id(roll->id());
   const int nstrips(roll->nstrips());
   double trArea(0.0);
   double trStripArea(0.0);
 
   if (me0Id.region() == 0)
   {
     throw cms::Exception("Geometry")
         << "ME0Synchronizer::simulateNoise() - this ME0 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.));
 
   const float rSqrtR = rollRadius * sqrt(rollRadius);
 
   //calculate noise from model
   double averageNeutralNoiseRatePerRoll = 0.;
   double averageNoiseElectronRatePerRoll = 0.;
   double averageNoiseRatePerRoll = 0.;
   if (me0Id.station() != 1)
   {
     throw cms::Exception("Geometry") << "ME0SimpleModel::simulateNoise() - this ME0 id is from station 1, which cannot happen: " <<roll->id() << "\n";
   }
   else
   {
     averageNeutralNoiseRatePerRoll = ME0NeuBkgParam0_ * rollRadius* std::exp(ME0NeuBkgParam1_/rSqrtR) + ME0NeuBkgParam2_/rSqrtR + ME0NeuBkgParam3_/(sqrt(rollRadius));
     //simulate electron background for ME0
     if (simulateElectronBkg_)
       averageNoiseElectronRatePerRoll = ME0ElecBkgParam0_ * rSqrtR* std::exp(ME0ElecBkgParam1_/rSqrtR) + ME0ElecBkgParam2_/rSqrtR + ME0ElecBkgParam3_/(sqrt(rollRadius));
     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)
     {
       int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, aveIntrinsicNoisePerStrip); 
     
       for (int k = 0; k < randPoissonQ; k++ )
       {
         const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
     strips_.emplace(k+1,time_hit);
       }
     }
   }//end simulate intrinsic noise
 
   //simulate bkg contribution
   const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9);
   int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, averageNoise);
 
   for (int i = 0; i < randPoissonQ; ++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.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 ME0SimpleModel::simulateSignal	(	const ME0EtaPartition *	roll,
		const edm::PSimHitContainer &	simHits,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements ME0DigiModel.

Definition at line 56 of file ME0SimpleModel.cc.

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

 {
   stripDigiSimLinks_.clear();
   detectorHitMap_.clear();
   stripDigiSimLinks_ = StripDigiSimLinks(roll->id().rawId());
 
   theME0DigiSimLinks_.clear();
   theME0DigiSimLinks_ = ME0DigiSimLinks(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 me0 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);
     }
   }
 }