CMS 3D CMS Logo

CastorFastClusterProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: CastorFastClusterProducer
4 // Class: CastorFastClusterProducer
5 //
13 //
14 // Original Author: Hans Van Haevermaet
15 // Created: Thu Mar 13 12:00:56 CET 2008
16 // $Id: CastorFastClusterProducer.cc,v 1.2 2011/03/04 11:00:55 hvanhaev Exp $
17 //
18 //
19 
20 // system include files
21 #include <memory>
22 #include <vector>
23 #include <iostream>
24 #include <sstream>
25 #include <TMath.h>
26 #include <TRandom3.h>
27 #include <TF1.h>
28 
29 // user include files
32 
35 
37 
39 
40 // Castorobject includes
44 
45 // genCandidate particle includes
48 
50 
51 //
52 // constructors and destructor
53 //
55  //register your products
56  produces<CastorClusterCollection>();
57 
58  //now do what ever other initialization is needed
59 }
60 
62  // do anything here that needs to be done at desctruction time
63  // (e.g. close files, deallocate resources etc.)
64 }
65 
66 //
67 // member functions
68 //
69 
70 // ------------ method called to produce the data ------------
72  using namespace edm;
73  using namespace reco;
74  using namespace std;
75  using namespace TMath;
76 
77  //
78  // Make CastorCluster objects
79  //
80 
81  //cout << "entering event" << endl;
82 
84  iEvent.getByLabel("genParticles", genParticles);
85 
86  // make pointer to towers that will be made
87  unique_ptr<CastorClusterCollection> CastorClusters(new CastorClusterCollection);
88 
89  /*
90  // declare castor array
91  double castorplus [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = 5.9
92  double castormin [4][16]; // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position - eta = -5.9
93  // set phi values of array sectors and everything else to zero
94  for (int j = 0; j < 16; j++) {
95  castorplus[3][j] = -2.94524 + j*0.3927;
96  castormin[3][j] = -2.94524 + j*0.3927;
97  castorplus[0][j] = 0.;
98  castormin[0][j] = 0.;
99  castorplus[1][j] = 0.;
100  castormin[1][j] = 0.;
101  castorplus[2][j] = 0.;
102  castormin[2][j] = 0.;
103  //castorplus[4][j] = 0.;
104  //castormin[4][j] = 0.;
105  }
106 
107  // declare properties vectors
108  vector<double> depthplus[16];
109  vector<double> depthmin[16];
110  vector<double> fhotplus [16];
111  vector<double> fhotmin [16];
112  vector<double> energyplus [16];
113  vector<double> energymin [16];
114 
115  for (int i=0;i<16;i++) {
116  depthplus[i].clear();
117  depthmin[i].clear();
118  fhotplus[i].clear();
119  fhotmin[i].clear();
120  energyplus[i].clear();
121  energymin[i].clear();
122  }
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 emEnergy = 0.;
138  double hadEnergy = 0.;
139  //double fhot = 0.;
140  //double depth = 0.;
141 
142  // add energies - em: if particle is e- or gamma
143  if (p.pdgId() == 11 || p.pdgId() == 22) {
144  while (gaus_E < 0.) {
145  // apply energy smearing with gaussian random generator
146  TRandom3 r(0);
147  // use sigma/E parametrization for the EM sections of CASTOR TB 2007 results
148  double sigma = p.energy() * (sqrt(pow(0.044, 2) + pow(0.513 / sqrt(p.energy()), 2)));
149  gaus_E = r.Gaus(p.energy(), sigma);
150  }
151 
152  // calculate electromagnetic electron/photon energy leakage
153  double tmax;
154  double a;
155  double cte;
156  if (p.pdgId() == 11) {
157  cte = -0.5;
158  } else {
159  cte = 0.5;
160  }
161  tmax = 1.0 * (log(gaus_E / 0.0015) + cte);
162  a = tmax * 0.5 + 1;
163  double leakage;
164  double x = 0.5 * 19.38;
165  leakage = gaus_E - gaus_E * Gamma(a, x);
166 
167  // add emEnergy
168  emEnergy = gaus_E - leakage;
169  // add hadEnergy leakage
170  hadEnergy = leakage;
171 
172  // make cluster
174  if (p.eta() > 0.) {
175  ClusterPoint temp(88.5, 5.9, p.phi());
176  pt1 = temp;
177  }
178  if (p.eta() < 0.) {
179  ClusterPoint temp(88.5, -5.9, p.phi());
180  pt1 = temp;
181  }
182  Point pt2(pt1);
183  CastorTowerRefVector refvector;
184  CastorClusters->push_back(
185  reco::CastorCluster(gaus_E, pt2, emEnergy, hadEnergy, emEnergy / gaus_E, 0., 0., 0., 0., refvector));
186 
187  } else {
188  while (gaus_E < 0.) {
189  // apply energy smearing with gaussian random generator
190  TRandom3 r(0);
191  // use sigma/E parametrization for the HAD sections of CASTOR TB 2007 results
192  double sigma = p.energy() * (sqrt(pow(0.121, 2) + pow(1.684 / sqrt(p.energy()), 2)));
193  gaus_E = r.Gaus(p.energy(), sigma);
194  }
195 
196  // add hadEnergy
197  hadEnergy = gaus_E;
198 
199  // make cluster
201  if (p.eta() > 0.) {
202  ClusterPoint temp(88.5, 5.9, p.phi());
203  pt1 = temp;
204  }
205  if (p.eta() < 0.) {
206  ClusterPoint temp(88.5, -5.9, p.phi());
207  pt1 = temp;
208  }
209  Point pt2(pt1);
210  CastorTowerRefVector refvector;
211  CastorClusters->push_back(reco::CastorCluster(gaus_E, pt2, 0., hadEnergy, 0., 0., 0., 0., 0., refvector));
212  }
213 
214  /*
215  // make tower
216 
217  // set sector
218  int sector = -1;
219  for (int j = 0; j < 16; j++) {
220  double a = -M_PI + j*0.3927;
221  double b = -M_PI + (j+1)*0.3927;
222  if ( (p.phi() > a) && (p.phi() < b)) {
223  sector = j;
224  }
225  }
226 
227  // set eta
228  if (p.eta() > 0) {
229  castorplus[0][sector] = castorplus[0][sector] + gaus_E;
230  castorplus[1][sector] = castorplus[1][sector] + emEnergy;
231  castorplus[2][sector] = castorplus[2][sector] + hadEnergy;
232 
233  depthplus[sector].push_back(depth);
234  fhotplus[sector].push_back(fhot);
235  energyplus[sector].push_back(gaus_E);
236  //cout << "filled vectors" << endl;
237  //cout << "energyplus size = " << energyplus[sector].size() << endl;
238  //cout << "depthplus size = " << depthplus[sector].size() << endl;
239  //cout << "fhotplus size = " << fhotplus[sector].size() << endl;
240 
241  } else {
242  castormin[0][sector] = castormin[0][sector] + gaus_E;
243  castormin[1][sector] = castormin[1][sector] + emEnergy;
244  castormin[2][sector] = castormin[2][sector] + hadEnergy;
245 
246 
247  depthmin[sector].push_back(depth);
248  fhotmin[sector].push_back(fhot);
249  energymin[sector].push_back(gaus_E);
250  //cout << "filled vectors" << endl;
251 
252  }
253  */
254  }
255  }
256 
257  /*
258  // substract pedestals/noise
259  for (int j = 0; j < 16; j++) {
260  double hadnoise = 0.;
261  for (int i=0;i<12;i++) {
262  hadnoise = hadnoise + make_noise();
263  }
264  castorplus[0][j] = castorplus[0][j] - hadnoise - make_noise() - make_noise();
265  castormin[0][j] = castormin[0][j] - hadnoise - make_noise() - make_noise();
266  castorplus[1][j] = castorplus[1][j] - make_noise() - make_noise();
267  castormin[1][j] = castormin[1][j] - make_noise() - make_noise();
268  castorplus[2][j] = castorplus[2][j] - hadnoise;
269  castormin[2][j] = castormin[2][j] - hadnoise;
270 
271  // set possible negative values to zero
272  if (castorplus[0][j] < 0.) castorplus[0][j] = 0.;
273  if (castormin[0][j] < 0.) castormin[0][j] = 0.;
274  if (castorplus[1][j] < 0.) castorplus[1][j] = 0.;
275  if (castormin[1][j] < 0.) castormin[1][j] = 0.;
276  if (castorplus[2][j] < 0.) castorplus[2][j] = 0.;
277  if (castormin[2][j] < 0.) castormin[2][j] = 0.;
278  }
279  */
280 
281  /*
282  // store towers from castor arrays
283  // eta = 5.9
284  for (int j=0;j<16;j++) {
285  if (castorplus[0][j] > 0.) {
286 
287  double fem = 0.;
288  fem = castorplus[1][j]/castorplus[0][j];
289  ClusterPoint pt1(88.5,5.9,castorplus[3][j]);
290  Point pt2(pt1);
291 
292  // parametrize depth and fhot from full sim
293  // get fit parameters from energy
294  // get random number according to distribution with fit parameters
295  double depth_mean = 0.;
296  double fhot_mean = 0.;
297  double sum_energy = 0.;
298 
299  //cout << "energyplus size = " << energyplus[j].size()<< endl;
300  for (size_t p = 0; p<energyplus[j].size();p++) {
301  depth_mean = depth_mean + depthplus[j][p]*energyplus[j][p];
302  fhot_mean = fhot_mean + fhotplus[j][p]*energyplus[j][p];
303  sum_energy = sum_energy + energyplus[j][p];
304  }
305  depth_mean = depth_mean/sum_energy;
306  fhot_mean = fhot_mean/sum_energy;
307  cout << "computed depth/fhot" << endl;
308 
309 
310  edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
311  CastorClusters->push_back(reco::CastorCluster(castorplus[0][j],pt2,castorplus[1][j],castorplus[2][j],fem,depth_mean,fhot_mean,refvector));
312  }
313  }
314  // eta = -5.9
315  for (int j=0;j<16;j++) {
316  if (castormin[0][j] > 0.) {
317  double fem = 0.;
318  fem = castormin[1][j]/castormin[0][j];
319  ClusterPoint pt1(88.5,-5.9,castormin[3][j]);
320  Point pt2(pt1);
321 
322  // parametrize depth and fhot from full sim
323  // get fit parameters from energy
324  // get random number according to distribution with fit parameters
325  double depth_mean = 0.;
326  double fhot_mean = 0.;
327  double sum_energy = 0.;
328 
329 
330  for (size_t p = 0; p<energymin[j].size();p++) {
331  depth_mean = depth_mean + depthmin[j][p]*energymin[j][p];
332  fhot_mean = fhot_mean + fhotmin[j][p]*energymin[j][p];
333  sum_energy = sum_energy + energymin[j][p];
334  }
335  depth_mean = depth_mean/sum_energy;
336  fhot_mean = fhot_mean/sum_energy;
337 
338 
339  edm::RefVector<edm::SortedCollection<CastorRecHit> > refvector;
340  CastorClusters->push_back(reco::CastorCluster(castormin[0][j],pt2,castormin[1][j],castormin[2][j],fem,depth_mean,fhot_mean,refvector));
341  }
342  }
343  */
344 
345  iEvent.put(std::move(CastorClusters));
346 }
347 
349  double result = -1.;
350  TRandom3 r2(0);
351  double mu_noise = 0.053; // GeV (from 1.214 ADC) per channel
352  double sigma_noise = 0.027; // GeV (from 0.6168 ADC) per channel
353 
354  while (result < 0.) {
355  result = r2.Gaus(mu_noise, sigma_noise);
356  }
357 
358  return result;
359 }
360 
361 //define this as a plug-in
mps_fire.i
i
Definition: mps_fire.py:428
genParticles2HepMC_cfi.genParticles
genParticles
Definition: genParticles2HepMC_cfi.py:4
CastorFastClusterProducer::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: CastorFastClusterProducer.cc:71
EDProducer.h
CastorFastClusterProducer
Definition: CastorFastClusterProducer.h:17
edm
HLT enums.
Definition: AlignableModifier.h:19
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
CastorFastClusterProducer::make_noise
double make_noise()
Definition: CastorFastClusterProducer.cc:348
CastorCluster.h
DDAxes::x
edm::RefVector< CastorTowerCollection >
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm::Handle
Definition: AssociativeIterator.h:50
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
GenParticle.h
MakerMacros.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
reco::CastorCluster
Definition: CastorCluster.h:24
CastorFastClusterProducer::CastorFastClusterProducer
CastorFastClusterProducer(const edm::ParameterSet &)
Definition: CastorFastClusterProducer.cc:54
HLT_FULL_cff.pt1
pt1
Definition: HLT_FULL_cff.py:9870
CastorFastClusterProducer::~CastorFastClusterProducer
~CastorFastClusterProducer() override
Definition: CastorFastClusterProducer.cc:61
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
edm::ParameterSet
Definition: ParameterSet.h:47
a
double a
Definition: hdecay.h:119
Event.h
CastorFastClusterProducer::CastorClusterCollection
std::vector< reco::CastorCluster > CastorClusterCollection
Definition: CastorFastClusterProducer.h:29
diffTwoXMLs.r2
r2
Definition: diffTwoXMLs.py:73
CastorFastClusterProducer.h
iEvent
int iEvent
Definition: GenABIO.cc:224
edm::EventSetup
Definition: EventSetup.h:57
HLT_FULL_cff.pt2
pt2
Definition: HLT_FULL_cff.py:9872
edm::RefVector::push_back
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
reco::Candidate
Definition: Candidate.h:27
alignCSCRings.r
r
Definition: alignCSCRings.py:93
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
Frameworkfwd.h
Point3D.h
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
mps_fire.result
result
Definition: mps_fire.py:311
CastorRecHit.h
Candidate.h
ParameterSet.h
HcalCastorDetId.h
edm::Event
Definition: Event.h:73
CastorFastClusterProducer::ClusterPoint
ROOT::Math::RhoEtaPhiPoint ClusterPoint
Definition: CastorFastClusterProducer.h:28