#include <ME0PreRecoGaussianModel.h>

Inheritance diagram for ME0PreRecoGaussianModel:

Public Member Functions
double	correctSigmaU (const ME0EtaPartition *, double)

	ME0PreRecoGaussianModel (const edm::ParameterSet &)

void	setup () override

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

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

	~ME0PreRecoGaussianModel () override

Public Member Functions inherited from ME0DigiPreRecoModel
void	fillDigis (int rollDetId, ME0DigiPreRecoCollection &)

const ME0Geometry *	getGeometry ()

void	setGeometry (const ME0Geometry *geom)

virtual	~ME0DigiPreRecoModel ()

Private Attributes
double	averageEfficiency_

double	constPhiSmearing_

bool	corr

bool	digitizeOnlyMuons_

std::vector< double >	eleBkg

double	error_u

double	error_v

bool	etaproj

bool	gaussianSmearing_

double	instLumi_

int	maxBunch_

int	minBunch_

std::vector< double >	neuBkg

double	rateFact_

double	referenceInstLumi_

double	sigma_t

double	sigma_u

double	sigma_v

bool	simulateElectronBkg_

bool	simulateNeutralBkg_

Additional Inherited Members
Protected Member Functions inherited from ME0DigiPreRecoModel
	ME0DigiPreRecoModel (const edm::ParameterSet &)

Protected Attributes inherited from ME0DigiPreRecoModel
std::set< ME0DigiPreReco >	digi_

const ME0Geometry *	geometry_

Detailed Description

Class for the ME0 Gaussian response simulation as pre-reco step

Definition at line 17 of file ME0PreRecoGaussianModel.h.

Constructor & Destructor Documentation

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

Definition at line 19 of file ME0PreRecoGaussianModel.cc.

References eleBkg, and neuBkg.

     : ME0DigiPreRecoModel(config),
       sigma_t(config.getParameter<double>("timeResolution")),
       sigma_u(config.getParameter<double>("phiResolution")),
       sigma_v(config.getParameter<double>("etaResolution")),
       error_u(config.getParameter<double>("phiError")),
       error_v(config.getParameter<double>("etaError")),
       gaussianSmearing_(config.getParameter<bool>("gaussianSmearing")),
       constPhiSmearing_(config.getParameter<bool>("constantPhiSpatialResolution")),
       corr(config.getParameter<bool>("useCorrelation")),
       etaproj(config.getParameter<bool>("useEtaProjectiveGEO")),
       digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
       averageEfficiency_(config.getParameter<double>("averageEfficiency")),
       // simulateIntrinsicNoise_(config.getParameter<bool>("simulateIntrinsicNoise")),
       // averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
       simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
       simulateNeutralBkg_(config.getParameter<bool>("simulateNeutralBkg")),
       minBunch_(config.getParameter<int>("minBunch")),
       maxBunch_(config.getParameter<int>("maxBunch")),
       instLumi_(config.getParameter<double>("instLumi")),
       rateFact_(config.getParameter<double>("rateFact")),
       referenceInstLumi_(config.getParameter<double>("referenceInstLumi")) {
   // polynomial parametrisation of neutral (n+g) and electron background
   // This is the background for an Instantaneous Luminosity of L = 5E34 cm^-2 s^-1
   neuBkg.push_back(0.00386257);
   neuBkg.push_back(6344.65);
   neuBkg.push_back(16627700);
   neuBkg.push_back(-102098);
 
   eleBkg.push_back(0.00171409);
   eleBkg.push_back(4900.56);
   eleBkg.push_back(710909);
   eleBkg.push_back(-4327.25);
 }

ME0PreRecoGaussianModel::~ME0PreRecoGaussianModel ( )

override

Definition at line 54 of file ME0PreRecoGaussianModel.cc.

54 {}

Member Function Documentation

double ME0PreRecoGaussianModel::correctSigmaU	(	const ME0EtaPartition *	roll,
		double	y
	)

Definition at line 336 of file ME0PreRecoGaussianModel.cc.

References ME0EtaPartitionSpecs::parameters(), Gflash::Rmax, sigma_u, ME0EtaPartition::specs(), ME0EtaPartition::topology(), and y.

Referenced by simulateNoise(), and simulateSignal().

                                                                                    {
   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
   auto& parameters(roll->specs()->parameters());
   double height(parameters[2]);        // height     = height from Center of Roll
   double rollRadius = top_->radius();  // rollRadius = Radius at Center of Roll
   double Rmax = rollRadius + height;   // MaxRadius  = Radius at top of Roll
   double Rx = rollRadius + y;          // y in [-height,+height]
   double sigma_u_new = Rx / Rmax * sigma_u;
   return sigma_u_new;
 }

void ME0PreRecoGaussianModel::setup ( )

inlineoverridevirtual

Implements ME0DigiPreRecoModel.

Definition at line 26 of file ME0PreRecoGaussianModel.h.

26 {}

void ME0PreRecoGaussianModel::simulateNoise	(	const ME0EtaPartition *	roll,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements ME0DigiPreRecoModel.

Definition at line 122 of file ME0PreRecoGaussianModel.cc.

References makePileupJSON::bx, bxwidth, constPhiSmearing_, corr, correctSigmaU(), ME0DigiPreRecoModel::digi_, alignCSCRings::e, eleBkg, error_u, error_v, Exception, funct::exp(), ME0EtaPartition::id(), instLumi_, maxBunch_, minBunch_, neuBkg, ME0EtaPartitionSpecs::parameters(), rateFact_, DetId::rawId(), referenceInstLumi_, sigma_u, sigma_v, simulateElectronBkg_, simulateNeutralBkg_, funct::sin(), ME0EtaPartition::specs(), mathSSE::sqrt(), funct::tan(), and ME0EtaPartition::topology().

                                                                                                      {
   double trArea(0.0);
   const ME0DetId me0Id(roll->id());
 
   // Extract detailed information from the Strip Topology:
   // base_bottom, base_top, height, strips, pads
   // note that (0,0) is in the middle of the roll ==> all param are at all half length
   if (me0Id.region() == 0) {
     throw cms::Exception("Geometry") << "Asking TrapezoidalStripTopology from a ME0 will fail";
   }  // not sure we really need this
   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
 
   auto& parameters(roll->specs()->parameters());
   double bottomLength(parameters[0]);
   bottomLength = 2 * bottomLength;  // bottom is largest length, so furtest away from beamline
   double topLength(parameters[1]);
   topLength = 2 * topLength;  // top is shortest length, so closest to beamline
   double height(parameters[2]);
   height = 2 * height;
   double myTanPhi = (topLength - bottomLength) / (height * 2);
   double rollRadius = top_->radius();
   trArea = height * (topLength + bottomLength) / 2.0;
 
   // Divide the detector area in different strips
   // take smearing in y-coord as height for each strip
   double initialHeight = sigma_v;
   if (sigma_v < 1.0)
     initialHeight = 1.0;
   double heightIt = initialHeight;
   int heightbins = height / heightIt;  // round down
 
   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: roll with id = " << roll->id();
   edm::LogVerbatim("ME0PreRecoGaussianModelNoise") << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId()
                                                    << "] :: extracting parameters from the TrapezoidalStripTopology";
   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: bottom = " << bottomLength
       << " [cm] top  = " << topLength << " [cm] height = " << height << " [cm]"
       << " area  = " << trArea << " [cm^2] Rmid = " << rollRadius
       << " [cm] => Rmin = " << rollRadius - height * 1.0 / 2.0 << " [cm] Rmax = " << rollRadius + height * 1.0 / 2.0
       << " [cm]";
   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: heightbins = " << heightbins;
 
   for (int hx = 0; hx < heightbins; ++hx) {
     double bottomIt = bottomLength + hx * 2 * tan(10. / 180 * 3.14) * heightIt;
     double topIt = bottomLength + (hx + 1) * 2 * tan(10. / 180 * 3.14) * heightIt;
     if (hx == heightbins - 1) {
       topIt = topLength;  // last bin ... make strip a bit larger to cover entire roll
       heightIt = height - hx * heightIt;
     }
     double areaIt = heightIt * (bottomIt + topIt) * 1.0 / 2;
 
     edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
         << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: height = " << std::setw(12) << heightIt
         << " [cm] bottom = " << std::setw(12) << bottomIt << " [cm]"
         << " top = " << std::setw(12) << topIt << " [cm] area = " << std::setw(12) << areaIt
         << " [cm^2] || sin(10) = " << sin(10. / 180 * 3.14);
 
     double myRandY = CLHEP::RandFlat::shoot(engine);
     double y0_rand = (hx + myRandY) * heightIt;  // Y coord, measured from the bottom of the roll
     if (hx == heightbins - 1)
       y0_rand = hx * initialHeight + myRandY * heightIt;
     double yy_rand =
         (y0_rand -
          height * 1.0 / 2);  // Y coord, measured from the middle of the roll, which is the Y coord in Local Coords
     double yy_glob = rollRadius + yy_rand;  // R coord in Global Coords
     // max length in x for given y coordinate (cfr trapezoidal eta partition)
     const float rSqrtR = yy_glob * sqrt(yy_glob);
     double xMax = topLength / 2.0 - (height / 2.0 - yy_rand) * myTanPhi;
 
     double sigma_u_new = sigma_u;
     if (constPhiSmearing_)
       sigma_u_new = correctSigmaU(roll, yy_rand);
 
     // 1) Intrinsic Noise ... Not implemented right now
     // ------------------------------------------------
     // if (simulateIntrinsicNoise_)
     // {
     // }
 
     // 2) Background Noise
     // ----------------------------
 
     // 2a) electron background
     // -----------------------
     if (simulateElectronBkg_) {
       // Extract / Calculate the Average Electron Rate
       // for the given global Y coord from Parametrization
       double averageElectronRatePerRoll =
           eleBkg[0] * rSqrtR * std::exp(eleBkg[1] / rSqrtR) + eleBkg[2] / rSqrtR + eleBkg[3] / (sqrt(yy_glob));
       // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
       averageElectronRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
 
       // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll in this bx
       const double averageElecRate(averageElectronRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
                                    areaIt);
 
       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
           << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
           << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageElectronRatePerRoll
           << " [Hz/cm^2]"
           << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
           << averageElecRate << " [hits]";
 
       bool ele_eff = (CLHEP::RandFlat::shoot(engine) < averageElecRate) ? true : false;
 
       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
           << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
           << " => local y = " << std::setw(12) << yy_rand << " [cm]"
           << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
           << averageElecRate << " => efficient? " << ele_eff << std::endl;
 
       if (ele_eff) {
         //calculate xx_rand at a given yy_rand
         double myRandX = CLHEP::RandFlat::shoot(engine);
         double xx_rand = 2 * xMax * (myRandX - 0.5);
         double ex = sigma_u_new;
         double ey = sigma_v;
         double corr = 0.;
         // extract random BX
         double myrandBX = CLHEP::RandFlat::shoot(engine);
         int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
         // extract random time in this BX
         double myrandT = CLHEP::RandFlat::shoot(engine);
         double minBXtime = (bx - 0.5) * bxwidth;  // double maxBXtime = (bx+0.5)*bxwidth;
         double time = myrandT * bxwidth + minBXtime;
         double myrandP = CLHEP::RandFlat::shoot(engine);
         int pdgid = 0;
         if (myrandP <= 0.5)
           pdgid = -11;  // electron
         else
           pdgid = 11;  // positron
         if (ex == 0)
           ex = error_u;  //errors cannot be zero
         if (ey == 0)
           ey = error_v;
 
         digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
 
         edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
             << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] =====> electron hit in " << roll->id()
             << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
             << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
       }
     }  // end if electron bkg
 
     // 2b) neutral (n+g) background
     // ----------------------------
     if (simulateNeutralBkg_) {
       // Extract / Calculate the Average Neutral Rate
       // for the given global Y coord from Parametrization
       double averageNeutralRatePerRoll =
           neuBkg[0] * yy_glob * std::exp(neuBkg[1] / rSqrtR) + neuBkg[2] / rSqrtR + neuBkg[3] / (sqrt(yy_glob));
       // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
       averageNeutralRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
 
       // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll
       const double averageNeutrRate(averageNeutralRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
                                     areaIt);
 
       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
           << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
           << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageNeutralRatePerRoll
           << " [Hz/cm^2]"
           << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
           << averageNeutrRate << " [hits]";
 
       bool neu_eff = (CLHEP::RandFlat::shoot(engine) < averageNeutrRate) ? true : false;
 
       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
           << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
           << " => local y = " << std::setw(12) << yy_rand << " [cm]"
           << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
           << averageNeutrRate << " => efficient? " << neu_eff << std::endl;
 
       if (neu_eff) {
         //calculate xx_rand at a given yy_rand
         double myRandX = CLHEP::RandFlat::shoot(engine);
         double xx_rand = 2 * xMax * (myRandX - 0.5);
         double ex = sigma_u_new;
         double ey = sigma_v;
         double corr = 0.;
         // extract random BX
         double myrandBX = CLHEP::RandFlat::shoot(engine);
         int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
         // extract random time in this BX
         double myrandT = CLHEP::RandFlat::shoot(engine);
         double minBXtime = (bx - 0.5) * bxwidth;
         double time = myrandT * bxwidth + minBXtime;
         int pdgid = 0;
         double myrandP = CLHEP::RandFlat::shoot(engine);
         if (myrandP <= 0.08)
           pdgid = 2112;  // neutrons: GEM sensitivity for neutrons: 0.08%
         else
           pdgid =
               22;  // photons:  GEM sensitivity for photons:  1.04% ==> neutron fraction = (0.08 / 1.04) = 0.077 = 0.08
         if (ex == 0)
           ex = error_u;  //errors cannot be zero
         if (ey == 0)
           ey = error_v;
 
         digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
 
         edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
             << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] ======> neutral hit in " << roll->id()
             << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
             << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
       }
 
     }  // end if neutral bkg
   }    // end loop over strips (= pseudo rolls)
 }

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

overridevirtual

Implements ME0DigiPreRecoModel.

Definition at line 56 of file ME0PreRecoGaussianModel.cc.

References funct::abs(), averageEfficiency_, makePileupJSON::bx, bxwidth, constPhiSmearing_, corr, correctSigmaU(), ME0DigiPreRecoModel::digi_, digitizeOnlyMuons_, mps_splice::entry, error_u, error_v, gaussianSmearing_, ME0EtaPartition::id(), maxBunch_, minBunch_, sigma_t, sigma_u, sigma_v, x, and y.

                                                                            {
   for (const auto& hit : simHits) {
     // Digitize only Muons?
     if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
       continue;
     // Digitize only in [minBunch,maxBunch] window
     // window is: [(2n+1)*bxw/2, (2n+3)*bxw/2], n = [minBunch, maxBunch]
     if (hit.timeOfFlight() < (2 * minBunch_ + 1) * bxwidth * 1.0 / 2 ||
         hit.timeOfFlight() > (2 * maxBunch_ + 3) * bxwidth * 1.0 / 2)
       continue;
     // is GEM efficient?
     if (CLHEP::RandFlat::shoot(engine) > averageEfficiency_)
       continue;
     // create digi
     auto entry = hit.entryPoint();
     double x = 0.0, y = 0.0;
 
     double sigma_u_new = sigma_u;
     if (constPhiSmearing_)
       sigma_u_new = correctSigmaU(roll, entry.y());
 
     if (gaussianSmearing_) {  // Gaussian Smearing
       x = CLHEP::RandGaussQ::shoot(engine, entry.x(), sigma_u_new);
       y = CLHEP::RandGaussQ::shoot(engine, entry.y(), sigma_v);
     } else {  // Uniform Smearing ... use the sigmas as boundaries
       x = entry.x() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_u_new;
       y = entry.y() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_v;
     }
     double ex = sigma_u_new;
     double ey = sigma_v;
     double corr = 0.;
     double tof = CLHEP::RandGaussQ::shoot(engine, hit.timeOfFlight(), sigma_t);
     int pdgid = hit.particleType();
     if (ex == 0)
       ex = error_u;  //errors cannot be zero
     if (ey == 0)
       ey = error_v;
 
     int evtId = hit.eventId().event();
     int bx = hit.eventId().bunchCrossing();
     int procType = hit.processType();
     int res = 1;
     if (!(evtId == 0 && bx == 0 && procType == 0))
       res = 2;
 
     digi_.emplace(x, y, ex, ey, corr, tof, pdgid, res);
 
     edm::LogVerbatim("ME0PreRecoGaussianModel")
         << "[ME0PreRecoDigi :: simulateSignal] :: simhit in " << roll->id() << " at loc x = " << std::setw(8)
         << entry.x() << " [cm]"
         << " loc y = " << std::setw(8) << entry.y() << " [cm] time = " << std::setw(8) << hit.timeOfFlight()
         << " [ns] pdgid = " << std::showpos << std::setw(4) << pdgid;
     edm::LogVerbatim("ME0PreRecoGaussianModel")
         << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " at loc x = " << std::setw(8) << x
         << " [cm] loc y = " << std::setw(8) << y << " [cm]"
         << " time = " << std::setw(8) << tof << " [ns]";
     edm::LogVerbatim("ME0PreRecoGaussianModel")
         << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " with DX = " << std::setw(8)
         << (entry.x() - x) << " [cm]"
         << " DY = " << std::setw(8) << (entry.y() - y) << " [cm] DT = " << std::setw(8) << (hit.timeOfFlight() - tof)
         << " [ns]";
   }
 }