CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CastorFastTowerProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: CastorFastTowerProducer
4 // Class: CastorFastTowerProducer
5 //
13 //
14 // Original Author: Hans Van Haevermaet
15 // Created: Thu Mar 13 12:00:56 CET 2008
16 // $Id: CastorFastTowerProducer.cc,v 1.1 2009/03/27 21:59:06 hvanhaev Exp $
17 //
18 //
19 
20 
21 // system include files
22 #include <memory>
23 #include <vector>
24 #include <iostream>
25 #include <sstream>
26 #include <TMath.h>
27 #include <TRandom3.h>
28 #include <TF1.h>
29 
30 // user include files
33 
36 
38 
40 
41 // Castorobject includes
44 
45 // genCandidate particle includes
48 
50 
51 
52 //
53 // constructors and destructor
54 //
56 {
57  //register your products
58  produces<CastorTowerCollection>();
59 
60  //now do what ever other initialization is needed
61 
62 }
63 
64 
66 {
67 
68  // do anything here that needs to be done at desctruction time
69  // (e.g. close files, deallocate resources etc.)
70 
71 }
72 
73 
74 //
75 // member functions
76 //
77 
78 // ------------ method called to produce the data ------------
79 void
81 {
82  using namespace edm;
83  using namespace reco;
84  using namespace std;
85  using namespace TMath;
86 
87  //
88  // Make CastorTower objects
89  //
90 
91  //cout << "entering event" << endl;
92 
94  iEvent.getByLabel("genParticles", genParticles);
95 
96  // make pointer to towers that will be made
97  auto_ptr<CastorTowerCollection> CastorTowers (new CastorTowerCollection);
98 
99  // declare castor array
100  double castorplus [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
101  double castormin [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
102  // set phi values of array sectors and everything else to zero
103  for (int j = 0; j < 16; j++) {
104  castorplus[3][j] = -2.94524 + j*0.3927;
105  castormin[3][j] = -2.94524 + j*0.3927;
106  castorplus[0][j] = 0.;
107  castormin[0][j] = 0.;
108  castorplus[1][j] = 0.;
109  castormin[1][j] = 0.;
110  castorplus[2][j] = 0.;
111  castormin[2][j] = 0.;
112  //castorplus[4][j] = 0.;
113  //castormin[4][j] = 0.;
114  }
115 
116  // declare properties vectors
117  vector<double> depthplus[16];
118  vector<double> depthmin[16];
119  vector<double> fhotplus [16];
120  vector<double> fhotmin [16];
121  vector<double> energyplus [16];
122  vector<double> energymin [16];
123  //vector<double> femplus [16];
124  //vector<double> femmin [16];
125 
126  for (int i=0;i<16;i++) {
127  depthplus[i].clear();
128  depthmin[i].clear();
129  fhotplus[i].clear();
130  fhotmin[i].clear();
131  energyplus[i].clear();
132  energymin[i].clear();
133  //femplus[i].clear();
134  //femmin[i].clear();
135  }
136 
137  //cout << "declared everything" << endl;
138 
139  // start particle loop
140  for (size_t i = 0; i < genParticles->size(); ++i) {
141  const Candidate & p = (*genParticles)[i];
142 
143  // select particles in castor
144  if ( fabs(p.eta()) > 5.2 && fabs(p.eta()) < 6.6) {
145 
146  //cout << "found particle in castor, start calculating" << endl;
147 
148  // declare energies
149  //double gaus_E = -1.;
150  double energy = -1.;
151  double emEnergy = 0.;
152  double hadEnergy = 0.;
153  double fhot = 0.;
154  double depth = 0.;
155  //double fem = 0.;
156 
157  // add energies - em: if particle is e- or gamma
158  if (p.pdgId() == 11 || p.pdgId() == 22) {
159 
160  // calculate primary tower energy for electrons
161  while ( energy < 0.) {
162  // apply energy smearing with gaussian random generator
163  TRandom3 r(0);
164  // use sigma/E parametrization from the full simulation
165  double mean = 1.0024*p.energy() - 0.3859;
166  double sigma = 0.0228*p.energy() + 2.1061;
167  energy = r.Gaus(mean,sigma);
168  }
169 
170  // calculate electromagnetic electron/photon energy leakage
171  double tmax;
172  double a;
173  double cte;
174  if ( p.pdgId() == 11) { cte = -0.5; } else { cte = 0.5; }
175  tmax = 1.0*(log(energy/0.0015)+cte);
176  a = tmax*0.5 + 1;
177  double leakage;
178  double x = 0.5*19.38;
179  leakage = energy - energy*Gamma(a,x);
180 
181  // add emEnergy
182  emEnergy = energy - leakage;
183  // add hadEnergy leakage
184  hadEnergy = leakage;
185 
186  // calculate EM depth from parametrization
187  double d0 = 0.2338 * pow(p.energy(),-0.1634);
188  double d1 = 5.4336 * pow(p.energy(),0.2410) + 14408.1025;
189  double d2 = 1.4692 * pow(p.energy(),0.1307) - 0.5216;
190  if (d0 < 0.) d0 = abs(d0);
191 
192  TF1 *fdepth = new TF1("fdepth","[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2])))",14400.,14460.);
193  fdepth->SetParameters(d0,d1,d2);
194  depth = fdepth->GetRandom();
195  fdepth->Delete();
196  if (p.eta() < 0.) depth = -1*depth;
197 
198  } else {
199 
200  // calculate primary tower energy for hadrons
201  while (energy < 0.) {
202  // apply energy smearing with gaussian random generator
203  TRandom3 r(0);
204  // use sigma/E parametrization from the full simulation
205  double mean = 0.8340*p.energy() - 8.5054;
206  double sigma = 0.1595*p.energy() + 3.1183;
207  energy = r.Gaus(mean,sigma);
208  }
209 
210  // add hadEnergy
211  hadEnergy = energy;
212 
213  // in the near future add fem parametrization
214 
215  // calculate depth for HAD particle from parametrization
216  double d0 = -0.000012 * p.energy() + 0.0661;
217  double d1 = 785.7524 * pow(p.energy(),0.0262) + 13663.4262;
218  double d2 = 9.8748 * pow(p.energy(),0.1720) + 37.0187;
219  if (d0 < 0.) d0 = abs(d0);
220 
221  TF1 *fdepth = new TF1("fdepth","[0] * TMath::Exp(-0.5*( (x-[1])/[2] + TMath::Exp(-(x-[1])/[2]) ))",14400.,15500.);
222  fdepth->SetParameters(d0,d1,d2);
223  depth = fdepth->GetRandom();
224  fdepth->Delete();
225  if (p.eta() < 0.) depth = -1*depth;
226 
227 
228  }
229 
230  // make tower
231 
232  // set sector
233  int sector = -1;
234  for (int j = 0; j < 16; j++) {
235  double a = -M_PI + j*0.3927;
236  double b = -M_PI + (j+1)*0.3927;
237  if ( (p.phi() > a) && (p.phi() < b)) {
238  sector = j;
239  }
240  }
241 
242  // set eta
243  if (p.eta() > 0) {
244  castorplus[0][sector] = castorplus[0][sector] + energy;
245  castorplus[1][sector] = castorplus[1][sector] + emEnergy;
246  castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
247 
248  depthplus[sector].push_back(depth);
249  fhotplus[sector].push_back(fhot);
250  energyplus[sector].push_back(energy);
251  //cout << "filled vectors" << endl;
252  //cout << "energyplus size = " << energyplus[sector].size() << endl;
253  //cout << "depthplus size = " << depthplus[sector].size() << endl;
254  //cout << "fhotplus size = " << fhotplus[sector].size() << endl;
255 
256  } else {
257  castormin[0][sector] = castormin[0][sector] + energy;
258  castormin[1][sector] = castormin[1][sector] + emEnergy;
259  castormin[2][sector] = castormin[2][sector] + hadEnergy;
260 
261 
262  depthmin[sector].push_back(depth);
263  fhotmin[sector].push_back(fhot);
264  energymin[sector].push_back(energy);
265  //cout << "filled vectors" << endl;
266 
267  }
268 
269  }
270 
271  }
272 
273 
274  // add and substract pedestals/noise
275  for (int j = 0; j < 16; j++) {
276  double hadnoise = 0.;
277  double emnoise = 0.;
278  for (int i=0;i<12;i++) {
279  hadnoise = hadnoise + make_noise();
280  if (i<2) emnoise = emnoise + make_noise();
281  }
282 
283  hadnoise = hadnoise - 12*0.053;
284  emnoise = emnoise - 2*0.053;
285  if ( hadnoise < 0.) hadnoise = 0.;
286  if ( emnoise < 0.) emnoise = 0.;
287  double totnoise = hadnoise + emnoise;
288 
289  // add random noise
290  castorplus[0][j] = castorplus[0][j] + totnoise;
291  castormin[0][j] = castormin[0][j] + totnoise;
292  castorplus[1][j] = castorplus[1][j] + emnoise;
293  castormin[1][j] = castormin[1][j] + emnoise;
294  castorplus[2][j] = castorplus[2][j] + hadnoise;
295  castormin[2][j] = castormin[2][j] + hadnoise;
296 
297  //cout << "after constant substraction" << endl;
298  //cout << "total noise = " << castorplus[0][j] << " em noise = " << castorplus[1][j] << " had noise = " << castorplus[2][j] << endl;
299  //cout << "fem should be = " << castorplus[1][j]/castorplus[0][j] << endl;
300 
301  }
302 
303 
304  // store towers from castor arrays
305  // eta = 5.9
306  for (int j=0;j<16;j++) {
307  if (castorplus[0][j] != 0.) {
308 
309  double fem = 0.;
310  fem = castorplus[1][j]/castorplus[0][j];
311  TowerPoint pt1(88.5,5.9,castorplus[3][j]);
312  Point pt2(pt1);
313 
314  //cout << "fem became = " << castorplus[1][j]/castorplus[0][j] << endl;
315 
316  double depth_mean = 0.;
317  double fhot_mean = 0.;
318  double sum_energy = 0.;
319 
320  //cout << "energyplus size = " << energyplus[j].size()<< endl;
321  for (size_t p = 0; p<energyplus[j].size();p++) {
322  depth_mean = depth_mean + depthplus[j][p]*energyplus[j][p];
323  fhot_mean = fhot_mean + fhotplus[j][p]*energyplus[j][p];
324  sum_energy = sum_energy + energyplus[j][p];
325  }
326  depth_mean = depth_mean/sum_energy;
327  fhot_mean = fhot_mean/sum_energy;
328  //cout << "computed depth/fhot" << endl;
329 
330 
331  //std::vector<CastorCell> usedCells;
332  CastorCellRefVector refvector;
333  CastorTowers->push_back(reco::CastorTower(castorplus[0][j],pt2,castorplus[1][j],castorplus[2][j],fem,depth_mean,fhot_mean,refvector));
334  }
335  }
336  // eta = -5.9
337  for (int j=0;j<16;j++) {
338  if (castormin[0][j] != 0.) {
339  double fem = 0.;
340  fem = castormin[1][j]/castormin[0][j];
341  TowerPoint pt1(88.5,-5.9,castormin[3][j]);
342  Point pt2(pt1);
343 
344  double depth_mean = 0.;
345  double fhot_mean = 0.;
346  double sum_energy = 0.;
347 
348 
349  for (size_t p = 0; p<energymin[j].size();p++) {
350  depth_mean = depth_mean + depthmin[j][p]*energymin[j][p];
351  fhot_mean = fhot_mean + fhotmin[j][p]*energymin[j][p];
352  sum_energy = sum_energy + energymin[j][p];
353  }
354  depth_mean = depth_mean/sum_energy;
355  fhot_mean = fhot_mean/sum_energy;
356 
357 
358  //std::vector<CastorCell> usedCells;
359  CastorCellRefVector refvector;
360  CastorTowers->push_back(reco::CastorTower(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
361  }
362  }
363 
364  iEvent.put(CastorTowers);
365 
366 }
367 
369  double result = -1.;
370  TRandom3 r2(0);
371  double mu_noise = 0.053; // GeV (from 1.214 ADC) per channel
372  double sigma_noise = 0.027; // GeV (from 0.6168 ADC) per channel
373 
374  while (result < 0.) {
375  result = r2.Gaus(mu_noise,sigma_noise);
376  }
377 
378  return result;
379 }
380 
381 
382 // ------------ method called once each job just before starting event loop ------------
383 void
385 {
386 }
387 
388 // ------------ method called once each job just after ending the event loop ------------
389 void
391 }
392 
393 //define this as a plug-in
int i
Definition: DBlmapReader.cc:9
virtual double energy() const =0
energy
virtual void beginRun(edm::Run &, edm::EventSetup const &)
CastorFastTowerProducer(const edm::ParameterSet &)
std::vector< reco::CastorTower > CastorTowerCollection
tuple d0
Definition: debug_cff.py:3
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
#define abs(x)
Definition: mlp_lapack.h:159
ROOT::Math::RhoEtaPhiPoint TowerPoint
tuple d1
Definition: debug_cff.py:7
virtual void produce(edm::Event &, const edm::EventSetup &)
int iEvent
Definition: GenABIO.cc:243
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
tuple result
Definition: query.py:137
int j
Definition: DBlmapReader.cc:9
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
virtual int pdgId() const =0
PDG identifier.
static const double tmax[3]
#define M_PI
Definition: BFit3D.cc:3
Log< T >::type log(const T &t)
Definition: Log.h:22
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
dbl * Gamma
Definition: mlp_gen.cc:38
void push_back(value_type const &ref)
Add a Ref&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:61
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: Run.h:31
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity