CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
List of all members | Public Member Functions | Private Types | Private Member Functions | Private Attributes
CastorFastTowerProducer Class Reference

#include <FastSimulation/ForwardDetectors/plugins/CastorFastTowerProducer.cc>

Inheritance diagram for CastorFastTowerProducer:
edm::stream::EDProducer<>

Public Member Functions

 CastorFastTowerProducer (const edm::ParameterSet &)
 
 ~CastorFastTowerProducer () override
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
 EDProducer (const EDProducer &)=delete
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 

Private Types

typedef std::vector< reco::CastorTowerCastorTowerCollection
 
typedef math::XYZPointD Point
 
typedef ROOT::Math::RhoEtaPhiPoint TowerPoint
 

Private Member Functions

double make_noise ()
 
void produce (edm::Event &, const edm::EventSetup &) override
 

Private Attributes

const edm::EDGetTokenT< reco::GenParticleCollectiontokGenPart_
 

Additional Inherited Members

- Public Types inherited from edm::stream::EDProducer<>
using CacheTypes = CacheContexts< T... >
 
using GlobalCache = typename CacheTypes::GlobalCache
 
using HasAbility = AbilityChecker< T... >
 
using InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
 
using LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
 
using LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
 
using LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
 
using RunCache = typename CacheTypes::RunCache
 
using RunContext = RunContextT< RunCache, GlobalCache >
 
using RunSummaryCache = typename CacheTypes::RunSummaryCache
 

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 17 of file CastorFastTowerProducer.h.

Member Typedef Documentation

◆ CastorTowerCollection

Definition at line 29 of file CastorFastTowerProducer.h.

◆ Point

Definition at line 27 of file CastorFastTowerProducer.h.

◆ TowerPoint

typedef ROOT::Math::RhoEtaPhiPoint CastorFastTowerProducer::TowerPoint
private

Definition at line 28 of file CastorFastTowerProducer.h.

Constructor & Destructor Documentation

◆ CastorFastTowerProducer()

CastorFastTowerProducer::CastorFastTowerProducer ( const edm::ParameterSet iConfig)
explicit

Definition at line 52 of file CastorFastTowerProducer.cc.

53  : tokGenPart_(consumes<reco::GenParticleCollection>(edm::InputTag{"genParticles"})) {
54  //register your products
55  produces<CastorTowerCollection>();
56 
57  //now do what ever other initialization is needed
58 }
const edm::EDGetTokenT< reco::GenParticleCollection > tokGenPart_

◆ ~CastorFastTowerProducer()

CastorFastTowerProducer::~CastorFastTowerProducer ( )
override

Definition at line 60 of file CastorFastTowerProducer.cc.

60  {
61  // do anything here that needs to be done at desctruction time
62  // (e.g. close files, deallocate resources etc.)
63 }

Member Function Documentation

◆ make_noise()

double CastorFastTowerProducer::make_noise ( )
private

Definition at line 354 of file CastorFastTowerProducer.cc.

References diffTwoXMLs::r2, and mps_fire::result.

Referenced by produce().

354  {
355  double result = -1.;
356  TRandom3 r2(0);
357  double mu_noise = 0.053; // GeV (from 1.214 ADC) per channel
358  double sigma_noise = 0.027; // GeV (from 0.6168 ADC) per channel
359 
360  while (result < 0.) {
361  result = r2.Gaus(mu_noise, sigma_noise);
362  }
363 
364  return result;
365 }

◆ produce()

void CastorFastTowerProducer::produce ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
overrideprivate

Definition at line 70 of file CastorFastTowerProducer.cc.

References a, funct::abs(), b, d0, d1, hcalRecHitTable_cff::depth, hcalRecHitTable_cff::energy, hcaldqm::quantity::fdepth, AJJGenJetFilter_cfi::genParticles, mps_fire::i, iEvent, dqmiolumiharvest::j, dqm-mbProfile::log, M_PI, make_noise(), SiStripPI::mean, eostools::move(), AlCaHLTBitMon_ParallelJobs::p, funct::pow(), HLT_2024v10_cff::pt1, HLT_2024v10_cff::pt2, edm::RefVector< C, T, F >::push_back(), alignCSCRings::r, nano_mu_digi_cff::sector, tmax, tokGenPart_, and x.

70  {
71  using namespace edm;
72  using namespace reco;
73  using namespace std;
74  using namespace TMath;
75 
76  //
77  // Make CastorTower objects
78  //
79 
80  //cout << "entering event" << endl;
81 
83 
84  // make pointer to towers that will be made
85  unique_ptr<CastorTowerCollection> CastorTowers(new CastorTowerCollection);
86 
87  // declare castor array
88  double castorplus[4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
89  double castormin[4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
90  // set phi values of array sectors and everything else to zero
91  for (int j = 0; j < 16; j++) {
92  castorplus[3][j] = -2.94524 + j * 0.3927;
93  castormin[3][j] = -2.94524 + j * 0.3927;
94  castorplus[0][j] = 0.;
95  castormin[0][j] = 0.;
96  castorplus[1][j] = 0.;
97  castormin[1][j] = 0.;
98  castorplus[2][j] = 0.;
99  castormin[2][j] = 0.;
100  //castorplus[4][j] = 0.;
101  //castormin[4][j] = 0.;
102  }
103 
104  // declare properties vectors
105  vector<double> depthplus[16];
106  vector<double> depthmin[16];
107  vector<double> fhotplus[16];
108  vector<double> fhotmin[16];
109  vector<double> energyplus[16];
110  vector<double> energymin[16];
111  //vector<double> femplus [16];
112  //vector<double> femmin [16];
113 
114  for (int i = 0; i < 16; i++) {
115  depthplus[i].clear();
116  depthmin[i].clear();
117  fhotplus[i].clear();
118  fhotmin[i].clear();
119  energyplus[i].clear();
120  energymin[i].clear();
121  //femplus[i].clear();
122  //femmin[i].clear();
123  }
124 
125  //cout << "declared everything" << endl;
126 
127  // start particle loop
128  for (size_t i = 0; i < genParticles->size(); ++i) {
129  const Candidate& p = (*genParticles)[i];
130 
131  // select particles in castor
132  if (fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) {
133  //cout << "found particle in castor, start calculating" << endl;
134 
135  // declare energies
136  //double gaus_E = -1.;
137  double energy = -1.;
138  double emEnergy = 0.;
139  double hadEnergy = 0.;
140  double fhot = 0.;
141  double depth = 0.;
142  //double fem = 0.;
143 
144  // add energies - em: if particle is e- or gamma
145  if (p.pdgId() == 11 || p.pdgId() == 22) {
146  // calculate primary tower energy for electrons
147  while (energy < 0.) {
148  // apply energy smearing with gaussian random generator
149  TRandom3 r(0);
150  // use sigma/E parametrization from the full simulation
151  double mean = 1.0024 * p.energy() - 0.3859;
152  double sigma = 0.0228 * p.energy() + 2.1061;
153  energy = r.Gaus(mean, sigma);
154  }
155 
156  // calculate electromagnetic electron/photon energy leakage
157  double tmax;
158  double a;
159  double cte;
160  if (p.pdgId() == 11) {
161  cte = -0.5;
162  } else {
163  cte = 0.5;
164  }
165  tmax = 1.0 * (log(energy / 0.0015) + cte);
166  a = tmax * 0.5 + 1;
167  double leakage;
168  double x = 0.5 * 19.38;
169  leakage = energy - energy * Gamma(a, x);
170 
171  // add emEnergy
172  emEnergy = energy - leakage;
173  // add hadEnergy leakage
174  hadEnergy = leakage;
175 
176  // calculate EM depth from parametrization
177  double d0 = 0.2338 * pow(p.energy(), -0.1634);
178  double d1 = 5.4336 * pow(p.energy(), 0.2410) + 14408.1025;
179  double d2 = 1.4692 * pow(p.energy(), 0.1307) - 0.5216;
180  if (d0 < 0.)
181  d0 = abs(d0);
182 
183  TF1* fdepth =
184  new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2])))", 14400., 14460.);
185  fdepth->SetParameters(d0, d1, d2);
186  depth = fdepth->GetRandom();
187  fdepth->Delete();
188  if (p.eta() < 0.)
189  depth = -1 * depth;
190 
191  } else {
192  // calculate primary tower energy for hadrons
193  while (energy < 0.) {
194  // apply energy smearing with gaussian random generator
195  TRandom3 r(0);
196  // use sigma/E parametrization from the full simulation
197  double mean = 0.8340 * p.energy() - 8.5054;
198  double sigma = 0.1595 * p.energy() + 3.1183;
199  energy = r.Gaus(mean, sigma);
200  }
201 
202  // add hadEnergy
203  hadEnergy = energy;
204 
205  // in the near future add fem parametrization
206 
207  // calculate depth for HAD particle from parametrization
208  double d0 = -0.000012 * p.energy() + 0.0661;
209  double d1 = 785.7524 * pow(p.energy(), 0.0262) + 13663.4262;
210  double d2 = 9.8748 * pow(p.energy(), 0.1720) + 37.0187;
211  if (d0 < 0.)
212  d0 = abs(d0);
213 
214  TF1* fdepth =
215  new TF1("fdepth", "[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2]) ))", 14400., 15500.);
216  fdepth->SetParameters(d0, d1, d2);
217  depth = fdepth->GetRandom();
218  fdepth->Delete();
219  if (p.eta() < 0.)
220  depth = -1 * depth;
221  }
222 
223  // make tower
224 
225  // set sector
226  int sector = -1;
227  for (int j = 0; j < 16; j++) {
228  double a = -M_PI + j * 0.3927;
229  double b = -M_PI + (j + 1) * 0.3927;
230  if ((p.phi() > a) && (p.phi() < b)) {
231  sector = j;
232  }
233  }
234 
235  // set eta
236  if (p.eta() > 0) {
237  castorplus[0][sector] = castorplus[0][sector] + energy;
238  castorplus[1][sector] = castorplus[1][sector] + emEnergy;
239  castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
240 
241  depthplus[sector].push_back(depth);
242  fhotplus[sector].push_back(fhot);
243  energyplus[sector].push_back(energy);
244  //cout << "filled vectors" << endl;
245  //cout << "energyplus size = " << energyplus[sector].size() << endl;
246  //cout << "depthplus size = " << depthplus[sector].size() << endl;
247  //cout << "fhotplus size = " << fhotplus[sector].size() << endl;
248 
249  } else {
250  castormin[0][sector] = castormin[0][sector] + energy;
251  castormin[1][sector] = castormin[1][sector] + emEnergy;
252  castormin[2][sector] = castormin[2][sector] + hadEnergy;
253 
254  depthmin[sector].push_back(depth);
255  fhotmin[sector].push_back(fhot);
256  energymin[sector].push_back(energy);
257  //cout << "filled vectors" << endl;
258  }
259  }
260  }
261 
262  // add and substract pedestals/noise
263  for (int j = 0; j < 16; j++) {
264  double hadnoise = 0.;
265  double emnoise = 0.;
266  for (int i = 0; i < 12; i++) {
267  hadnoise = hadnoise + make_noise();
268  if (i < 2)
269  emnoise = emnoise + make_noise();
270  }
271 
272  hadnoise = hadnoise - 12 * 0.053;
273  emnoise = emnoise - 2 * 0.053;
274  if (hadnoise < 0.)
275  hadnoise = 0.;
276  if (emnoise < 0.)
277  emnoise = 0.;
278  double totnoise = hadnoise + emnoise;
279 
280  // add random noise
281  castorplus[0][j] = castorplus[0][j] + totnoise;
282  castormin[0][j] = castormin[0][j] + totnoise;
283  castorplus[1][j] = castorplus[1][j] + emnoise;
284  castormin[1][j] = castormin[1][j] + emnoise;
285  castorplus[2][j] = castorplus[2][j] + hadnoise;
286  castormin[2][j] = castormin[2][j] + hadnoise;
287 
288  //cout << "after constant substraction" << endl;
289  //cout << "total noise = " << castorplus[0][j] << " em noise = " << castorplus[1][j] << " had noise = " << castorplus[2][j] << endl;
290  //cout << "fem should be = " << castorplus[1][j]/castorplus[0][j] << endl;
291  }
292 
293  // store towers from castor arrays
294  // eta = 5.9
295  for (int j = 0; j < 16; j++) {
296  if (castorplus[0][j] != 0.) {
297  double fem = 0.;
298  fem = castorplus[1][j] / castorplus[0][j];
299  TowerPoint pt1(88.5, 5.9, castorplus[3][j]);
300  Point pt2(pt1);
301 
302  //cout << "fem became = " << castorplus[1][j]/castorplus[0][j] << endl;
303 
304  double depth_mean = 0.;
305  double fhot_mean = 0.;
306  double sum_energy = 0.;
307 
308  //cout << "energyplus size = " << energyplus[j].size()<< endl;
309  for (size_t p = 0; p < energyplus[j].size(); p++) {
310  depth_mean = depth_mean + depthplus[j][p] * energyplus[j][p];
311  fhot_mean = fhot_mean + fhotplus[j][p] * energyplus[j][p];
312  sum_energy = sum_energy + energyplus[j][p];
313  }
314  depth_mean = depth_mean / sum_energy;
315  fhot_mean = fhot_mean / sum_energy;
316  //cout << "computed depth/fhot" << endl;
317 
318  //std::vector<CastorCell> usedCells;
320  CastorTowers->push_back(reco::CastorTower(
321  castorplus[0][j], pt2, castorplus[1][j], castorplus[2][j], fem, depth_mean, fhot_mean, refvector));
322  }
323  }
324  // eta = -5.9
325  for (int j = 0; j < 16; j++) {
326  if (castormin[0][j] != 0.) {
327  double fem = 0.;
328  fem = castormin[1][j] / castormin[0][j];
329  TowerPoint pt1(88.5, -5.9, castormin[3][j]);
330  Point pt2(pt1);
331 
332  double depth_mean = 0.;
333  double fhot_mean = 0.;
334  double sum_energy = 0.;
335 
336  for (size_t p = 0; p < energymin[j].size(); p++) {
337  depth_mean = depth_mean + depthmin[j][p] * energymin[j][p];
338  fhot_mean = fhot_mean + fhotmin[j][p] * energymin[j][p];
339  sum_energy = sum_energy + energymin[j][p];
340  }
341  depth_mean = depth_mean / sum_energy;
342  fhot_mean = fhot_mean / sum_energy;
343 
344  //std::vector<CastorCell> usedCells;
346  CastorTowers->push_back(reco::CastorTower(
347  castormin[0][j], pt2, castormin[1][j], castormin[2][j], fem, depth_mean, fhot_mean, refvector));
348  }
349  }
350 
351  iEvent.put(std::move(CastorTowers));
352 }
ROOT::Math::RhoEtaPhiPoint TowerPoint
int iEvent
Definition: GenABIO.cc:224
const edm::EDGetTokenT< reco::GenParticleCollection > tokGenPart_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define M_PI
static const double tmax[3]
static constexpr float d0
double b
Definition: hdecay.h:120
fixed size matrix
HLT enums.
Structure Point Contains parameters of Gaussian fits to DMRs.
double a
Definition: hdecay.h:121
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
static constexpr float d1
std::vector< CastorTower > CastorTowerCollection
collection of CastorTower objects
Definition: CastorTower.h:137
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
def move(src, dest)
Definition: eostools.py:511

Member Data Documentation

◆ tokGenPart_

const edm::EDGetTokenT<reco::GenParticleCollection> CastorFastTowerProducer::tokGenPart_
private

Definition at line 30 of file CastorFastTowerProducer.h.

Referenced by produce().