#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 18 of file ME0PreRecoGaussianModel.h.

Constructor & Destructor Documentation

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

Definition at line 20 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 56 of file ME0PreRecoGaussianModel.cc.

57 {

58 }

Member Function Documentation

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

Definition at line 290 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 29 of file ME0PreRecoGaussianModel.h.

29 {}

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

overridevirtual

Implements ME0DigiPreRecoModel.

Definition at line 112 of file ME0PreRecoGaussianModel.cc.

References bxwidth, constPhiSmearing_, corr, correctSigmaU(), ME0DigiPreRecoModel::digi_, MillePedeFileConverter_cfg::e, eleBkg, error_u, error_v, Exception, JetChargeProducer_cfi::exp, ME0EtaPartition::id(), instLumi_, createfilelist::int, maxBunch_, minBunch_, neuBkg, ME0EtaPartitionSpecs::parameters(), BPhysicsValidation_cfi::pdgid, rateFact_, DetId::rawId(), referenceInstLumi_, sigma_u, sigma_v, simulateElectronBkg_, simulateNeutralBkg_, funct::sin(), ME0EtaPartition::specs(), mathSSE::sqrt(), funct::tan(), ntuplemaker::time, ME0EtaPartition::topology(), and anotherprimaryvertexanalyzer_cfi::xMax.

 {
   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 60 of file ME0PreRecoGaussianModel.cc.

References funct::abs(), averageEfficiency_, bxwidth, constPhiSmearing_, corr, correctSigmaU(), ME0DigiPreRecoModel::digi_, digitizeOnlyMuons_, mps_splice::entry, error_u, error_v, gaussianSmearing_, ME0EtaPartition::id(), maxBunch_, minBunch_, BPhysicsValidation_cfi::pdgid, 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]";
 
 }
 }