1 |
|
---|
2 | /** \class Calorimeter
|
---|
3 | *
|
---|
4 | * Fills calorimeter towers, performs calorimeter resolution smearing,
|
---|
5 | * preselects towers hit by photons and creates energy flow objects.
|
---|
6 | *
|
---|
7 | * $Date: 2013-07-22 22:48:31 +0000 (Mon, 22 Jul 2013) $
|
---|
8 | * $Revision: 1239 $
|
---|
9 | *
|
---|
10 | *
|
---|
11 | * \author P. Demin - UCL, Louvain-la-Neuve
|
---|
12 | *
|
---|
13 | */
|
---|
14 |
|
---|
15 | #include "modules/Calorimeter.h"
|
---|
16 |
|
---|
17 | #include "classes/DelphesClasses.h"
|
---|
18 | #include "classes/DelphesFactory.h"
|
---|
19 | #include "classes/DelphesFormula.h"
|
---|
20 |
|
---|
21 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
22 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
23 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
24 |
|
---|
25 | #include "TMath.h"
|
---|
26 | #include "TString.h"
|
---|
27 | #include "TFormula.h"
|
---|
28 | #include "TRandom3.h"
|
---|
29 | #include "TObjArray.h"
|
---|
30 | #include "TDatabasePDG.h"
|
---|
31 | #include "TLorentzVector.h"
|
---|
32 |
|
---|
33 | #include <algorithm>
|
---|
34 | #include <stdexcept>
|
---|
35 | #include <iostream>
|
---|
36 | #include <sstream>
|
---|
37 |
|
---|
38 | using namespace std;
|
---|
39 |
|
---|
40 | //------------------------------------------------------------------------------
|
---|
41 |
|
---|
42 | Calorimeter::Calorimeter() :
|
---|
43 | fECalResolutionFormula(0), fHCalResolutionFormula(0),
|
---|
44 | fItParticleInputArray(0), fItTrackInputArray(0),
|
---|
45 | fTowerECalArray(0), fItTowerECalArray(0),
|
---|
46 | fTowerHCalArray(0), fItTowerHCalArray(0),
|
---|
47 | fTowerTrackArray(0), fItTowerTrackArray(0),
|
---|
48 | fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
|
---|
49 | fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
|
---|
50 | {
|
---|
51 | fECalResolutionFormula = new DelphesFormula;
|
---|
52 | fHCalResolutionFormula = new DelphesFormula;
|
---|
53 |
|
---|
54 | fTowerECalArray = new TObjArray;
|
---|
55 | fItTowerECalArray = fTowerECalArray->MakeIterator();
|
---|
56 | fTowerHCalArray = new TObjArray;
|
---|
57 | fItTowerHCalArray = fTowerHCalArray->MakeIterator();
|
---|
58 |
|
---|
59 | fTowerTrackArray = new TObjArray;
|
---|
60 | fItTowerTrackArray = fTowerTrackArray->MakeIterator();
|
---|
61 | fTowerECalTrackArray = new TObjArray;
|
---|
62 | fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
|
---|
63 | fTowerHCalTrackArray = new TObjArray;
|
---|
64 | fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
|
---|
65 | }
|
---|
66 |
|
---|
67 | //------------------------------------------------------------------------------
|
---|
68 |
|
---|
69 | Calorimeter::~Calorimeter()
|
---|
70 | {
|
---|
71 | if(fECalResolutionFormula) delete fECalResolutionFormula;
|
---|
72 | if(fHCalResolutionFormula) delete fHCalResolutionFormula;
|
---|
73 |
|
---|
74 | if(fTowerECalArray) delete fTowerECalArray;
|
---|
75 | if(fItTowerECalArray) delete fItTowerECalArray;
|
---|
76 | if(fTowerHCalArray) delete fTowerHCalArray;
|
---|
77 | if(fItTowerHCalArray) delete fItTowerHCalArray;
|
---|
78 |
|
---|
79 | if(fTowerTrackArray) delete fTowerTrackArray;
|
---|
80 | if(fItTowerTrackArray) delete fItTowerTrackArray;
|
---|
81 | if(fTowerECalTrackArray) delete fTowerECalTrackArray;
|
---|
82 | if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
|
---|
83 | if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
|
---|
84 | if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
|
---|
85 |
|
---|
86 | //------------------------------------------------------------------------------
|
---|
87 |
|
---|
88 | void Calorimeter::Init()
|
---|
89 | {
|
---|
90 | ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
|
---|
91 | Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
|
---|
92 | Double_t ecalFraction, hcalFraction;
|
---|
93 | TBinMap::iterator itEtaBin;
|
---|
94 | set< Double_t >::iterator itPhiBin;
|
---|
95 | vector< Double_t > *phiBins;
|
---|
96 |
|
---|
97 | // read eta and phi bins
|
---|
98 | param = GetParam("EtaPhiBins");
|
---|
99 | size = param.GetSize();
|
---|
100 | fBinMap.clear();
|
---|
101 | fEtaBins.clear();
|
---|
102 | fPhiBins.clear();
|
---|
103 | for(i = 0; i < size/2; ++i)
|
---|
104 | {
|
---|
105 | paramEtaBins = param[i*2];
|
---|
106 | sizeEtaBins = paramEtaBins.GetSize();
|
---|
107 | paramPhiBins = param[i*2 + 1];
|
---|
108 | sizePhiBins = paramPhiBins.GetSize();
|
---|
109 |
|
---|
110 | for(j = 0; j < sizeEtaBins; ++j)
|
---|
111 | {
|
---|
112 | for(k = 0; k < sizePhiBins; ++k)
|
---|
113 | {
|
---|
114 | fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
|
---|
115 | }
|
---|
116 | }
|
---|
117 | }
|
---|
118 |
|
---|
119 | // for better performance we transform map of sets to parallel vectors:
|
---|
120 | // vector< double > and vector< vector< double >* >
|
---|
121 | for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
|
---|
122 | {
|
---|
123 | fEtaBins.push_back(itEtaBin->first);
|
---|
124 | phiBins = new vector< double >(itEtaBin->second.size());
|
---|
125 | fPhiBins.push_back(phiBins);
|
---|
126 | phiBins->clear();
|
---|
127 | for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
|
---|
128 | {
|
---|
129 | phiBins->push_back(*itPhiBin);
|
---|
130 | }
|
---|
131 | }
|
---|
132 |
|
---|
133 | // read energy fractions for different particles
|
---|
134 | param = GetParam("EnergyFraction");
|
---|
135 | size = param.GetSize();
|
---|
136 |
|
---|
137 | // set default energy fractions values
|
---|
138 | fFractionMap.clear();
|
---|
139 | fFractionMap[0] = make_pair(0.0, 1.0);
|
---|
140 |
|
---|
141 | for(i = 0; i < size/2; ++i)
|
---|
142 | {
|
---|
143 | paramFractions = param[i*2 + 1];
|
---|
144 | sizeFractions = paramFractions.GetSize();
|
---|
145 |
|
---|
146 | ecalFraction = paramFractions[0].GetDouble();
|
---|
147 | hcalFraction = paramFractions[1].GetDouble();
|
---|
148 |
|
---|
149 | fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
|
---|
150 | }
|
---|
151 | /*
|
---|
152 | TFractionMap::iterator itFractionMap;
|
---|
153 | for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
|
---|
154 | {
|
---|
155 | cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
|
---|
156 | }
|
---|
157 | */
|
---|
158 | // read resolution formulas
|
---|
159 | fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
|
---|
160 | fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
|
---|
161 |
|
---|
162 | // import array with output from other modules
|
---|
163 | fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
|
---|
164 | fItParticleInputArray = fParticleInputArray->MakeIterator();
|
---|
165 |
|
---|
166 | fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
|
---|
167 | fItTrackInputArray = fTrackInputArray->MakeIterator();
|
---|
168 |
|
---|
169 | // create output arrays
|
---|
170 | fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
|
---|
171 | fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
|
---|
172 |
|
---|
173 | fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
|
---|
174 | fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
|
---|
175 | }
|
---|
176 |
|
---|
177 | //------------------------------------------------------------------------------
|
---|
178 |
|
---|
179 | void Calorimeter::Finish()
|
---|
180 | {
|
---|
181 | vector< vector< Double_t>* >::iterator itPhiBin;
|
---|
182 | if(fItParticleInputArray) delete fItParticleInputArray;
|
---|
183 | if(fItTrackInputArray) delete fItTrackInputArray;
|
---|
184 | for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
|
---|
185 | {
|
---|
186 | delete *itPhiBin;
|
---|
187 | }
|
---|
188 | }
|
---|
189 |
|
---|
190 | //------------------------------------------------------------------------------
|
---|
191 |
|
---|
192 | void Calorimeter::Process()
|
---|
193 | {
|
---|
194 | Candidate *particle, *track;
|
---|
195 | TLorentzVector position, momentum;
|
---|
196 | Short_t etaBin, phiBin, flags;
|
---|
197 | Int_t number;
|
---|
198 | Long64_t towerHit, towerEtaPhi, hitEtaPhi;
|
---|
199 | Double_t ecalFraction, hcalFraction;
|
---|
200 | Double_t ecalEnergy, hcalEnergy;
|
---|
201 | Int_t pdgCode;
|
---|
202 |
|
---|
203 | TFractionMap::iterator itFractionMap;
|
---|
204 |
|
---|
205 | vector< Double_t >::iterator itEtaBin;
|
---|
206 | vector< Double_t >::iterator itPhiBin;
|
---|
207 | vector< Double_t > *phiBins;
|
---|
208 |
|
---|
209 | vector< Long64_t >::iterator itTowerHits;
|
---|
210 |
|
---|
211 | DelphesFactory *factory = GetFactory();
|
---|
212 | fTowerHits.clear();
|
---|
213 | fECalFractions.clear();
|
---|
214 | fHCalFractions.clear();
|
---|
215 |
|
---|
216 | // loop over all particles
|
---|
217 | fItParticleInputArray->Reset();
|
---|
218 | number = -1;
|
---|
219 | while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
|
---|
220 | {
|
---|
221 | const TLorentzVector &particlePosition = particle->Position;
|
---|
222 | ++number;
|
---|
223 |
|
---|
224 | pdgCode = TMath::Abs(particle->PID);
|
---|
225 |
|
---|
226 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
227 | if(itFractionMap == fFractionMap.end())
|
---|
228 | {
|
---|
229 | itFractionMap = fFractionMap.find(0);
|
---|
230 | }
|
---|
231 |
|
---|
232 | ecalFraction = itFractionMap->second.first;
|
---|
233 | hcalFraction = itFractionMap->second.second;
|
---|
234 |
|
---|
235 | fECalFractions.push_back(ecalFraction);
|
---|
236 | fHCalFractions.push_back(hcalFraction);
|
---|
237 |
|
---|
238 | if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
|
---|
239 |
|
---|
240 | // find eta bin [1, fEtaBins.size - 1]
|
---|
241 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
|
---|
242 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
243 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
244 |
|
---|
245 | // phi bins for given eta bin
|
---|
246 | phiBins = fPhiBins[etaBin];
|
---|
247 |
|
---|
248 | // find phi bin [1, phiBins.size - 1]
|
---|
249 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
|
---|
250 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
251 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
252 |
|
---|
253 | flags = 0;
|
---|
254 | flags |= (ecalFraction >= 1.0E-9) << 1;
|
---|
255 | flags |= (hcalFraction >= 1.0E-9) << 2;
|
---|
256 | flags |= (pdgCode == 11 || pdgCode == 22) << 3;
|
---|
257 |
|
---|
258 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
|
---|
259 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
260 |
|
---|
261 | fTowerHits.push_back(towerHit);
|
---|
262 | }
|
---|
263 |
|
---|
264 | // loop over all tracks
|
---|
265 | fItTrackInputArray->Reset();
|
---|
266 | number = -1;
|
---|
267 | while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
|
---|
268 | {
|
---|
269 | const TLorentzVector &trackPosition = track->Position;
|
---|
270 | ++number;
|
---|
271 |
|
---|
272 | pdgCode = TMath::Abs(track->PID);
|
---|
273 |
|
---|
274 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
275 | if(itFractionMap == fFractionMap.end())
|
---|
276 | {
|
---|
277 | itFractionMap = fFractionMap.find(0);
|
---|
278 | }
|
---|
279 |
|
---|
280 | ecalFraction = itFractionMap->second.first;
|
---|
281 | hcalFraction = itFractionMap->second.second;
|
---|
282 |
|
---|
283 | // find eta bin [1, fEtaBins.size - 1]
|
---|
284 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
|
---|
285 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
286 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
287 |
|
---|
288 | // phi bins for given eta bin
|
---|
289 | phiBins = fPhiBins[etaBin];
|
---|
290 |
|
---|
291 | // find phi bin [1, phiBins.size - 1]
|
---|
292 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
|
---|
293 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
294 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
295 |
|
---|
296 | flags = 1;
|
---|
297 | flags |= (ecalFraction >= 1.0E-9) << 1;
|
---|
298 | flags |= (hcalFraction >= 1.0E-9) << 2;
|
---|
299 |
|
---|
300 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
|
---|
301 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
302 |
|
---|
303 | fTowerHits.push_back(towerHit);
|
---|
304 | }
|
---|
305 |
|
---|
306 | // all hits are sorted first by eta bin number, then by phi bin number,
|
---|
307 | // then by flags and then by particle or track number
|
---|
308 | sort(fTowerHits.begin(), fTowerHits.end());
|
---|
309 |
|
---|
310 | // loop over all hits
|
---|
311 | towerEtaPhi = 0;
|
---|
312 | fTower = 0;
|
---|
313 | for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
|
---|
314 | {
|
---|
315 | towerHit = (*itTowerHits);
|
---|
316 | flags = (towerHit >> 24) & 0x00000000000000FFLL;
|
---|
317 | number = (towerHit) & 0x0000000000FFFFFFLL;
|
---|
318 | hitEtaPhi = towerHit >> 32;
|
---|
319 |
|
---|
320 | if(towerEtaPhi != hitEtaPhi)
|
---|
321 | {
|
---|
322 | // switch to next tower
|
---|
323 | towerEtaPhi = hitEtaPhi;
|
---|
324 |
|
---|
325 | // finalize previous tower
|
---|
326 | FinalizeTower();
|
---|
327 |
|
---|
328 | // create new tower
|
---|
329 | fTower = factory->NewCandidate();
|
---|
330 |
|
---|
331 | phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
|
---|
332 | etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
|
---|
333 |
|
---|
334 | // phi bins for given eta bin
|
---|
335 | phiBins = fPhiBins[etaBin];
|
---|
336 |
|
---|
337 | // calculate eta and phi of the tower's center
|
---|
338 | fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
|
---|
339 | fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
|
---|
340 |
|
---|
341 | fTowerEdges[0] = fEtaBins[etaBin - 1];
|
---|
342 | fTowerEdges[1] = fEtaBins[etaBin];
|
---|
343 | fTowerEdges[2] = (*phiBins)[phiBin - 1];
|
---|
344 | fTowerEdges[3] = (*phiBins)[phiBin];
|
---|
345 |
|
---|
346 | fTowerECalEnergy = 0.0;
|
---|
347 | fTowerHCalEnergy = 0.0;
|
---|
348 |
|
---|
349 | fTowerPhotonHits = 0;
|
---|
350 |
|
---|
351 | fTowerAllHits = 0;
|
---|
352 | fTowerECalHits = 0;
|
---|
353 | fTowerHCalHits = 0;
|
---|
354 |
|
---|
355 | fTowerTrackAllHits = 0;
|
---|
356 | fTowerECalTrackHits = 0;
|
---|
357 | fTowerHCalTrackHits = 0;
|
---|
358 |
|
---|
359 | fTowerECalArray->Clear();
|
---|
360 | fTowerHCalArray->Clear();
|
---|
361 |
|
---|
362 | fTowerTrackArray->Clear();
|
---|
363 | fTowerECalTrackArray->Clear();
|
---|
364 | fTowerHCalTrackArray->Clear();
|
---|
365 | }
|
---|
366 |
|
---|
367 | // check for track hits
|
---|
368 | if(flags & 1)
|
---|
369 | {
|
---|
370 | track = static_cast<Candidate*>(fTrackInputArray->At(number));
|
---|
371 |
|
---|
372 | ++fTowerTrackAllHits;
|
---|
373 | fTowerTrackArray->Add(track);
|
---|
374 |
|
---|
375 | // check for track ECAL hits
|
---|
376 | if(flags & 2)
|
---|
377 | {
|
---|
378 | ++fTowerECalTrackHits;
|
---|
379 | fTowerECalTrackArray->Add(track);
|
---|
380 | }
|
---|
381 |
|
---|
382 | // check for track HCAL hits
|
---|
383 | if(flags & 4)
|
---|
384 | {
|
---|
385 | ++fTowerHCalTrackHits;
|
---|
386 | fTowerHCalTrackArray->Add(track);
|
---|
387 | }
|
---|
388 | continue;
|
---|
389 | }
|
---|
390 |
|
---|
391 | ++fTowerAllHits;
|
---|
392 | fTower->AddCandidate(particle);
|
---|
393 |
|
---|
394 | // check for ECAL hits
|
---|
395 | if(flags & 2)
|
---|
396 | {
|
---|
397 | ++fTowerECalHits;
|
---|
398 | fTowerECalArray->Add(particle);
|
---|
399 | }
|
---|
400 |
|
---|
401 | // check for HCAL hits
|
---|
402 | if(flags & 4)
|
---|
403 | {
|
---|
404 | ++fTowerHCalHits;
|
---|
405 | fTowerHCalArray->Add(particle);
|
---|
406 | }
|
---|
407 |
|
---|
408 | // check for photon and electron hits in current tower
|
---|
409 | if(flags & 8) ++fTowerPhotonHits;
|
---|
410 |
|
---|
411 | particle = static_cast<Candidate*>(fParticleInputArray->At(number));
|
---|
412 | momentum = particle->Momentum;
|
---|
413 |
|
---|
414 | // fill current tower
|
---|
415 | ecalEnergy = momentum.E() * fECalFractions[number];
|
---|
416 | hcalEnergy = momentum.E() * fHCalFractions[number];
|
---|
417 |
|
---|
418 | fTowerECalEnergy += ecalEnergy;
|
---|
419 | fTowerHCalEnergy += hcalEnergy;
|
---|
420 | }
|
---|
421 |
|
---|
422 | // finalize last tower
|
---|
423 | FinalizeTower();
|
---|
424 | }
|
---|
425 |
|
---|
426 | //------------------------------------------------------------------------------
|
---|
427 |
|
---|
428 | void Calorimeter::FinalizeTower()
|
---|
429 | {
|
---|
430 | Candidate *particle, *track, *tower;
|
---|
431 | Double_t energy, pt, eta, phi;
|
---|
432 | Double_t ecalEnergy, hcalEnergy;
|
---|
433 | TIterator *iterator;
|
---|
434 |
|
---|
435 | if(!fTower) return;
|
---|
436 |
|
---|
437 | // ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
438 | // if(ecalEnergy < 0.0) ecalEnergy = 0.0;
|
---|
439 |
|
---|
440 | ecalEnergy = LogNormal(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
441 |
|
---|
442 | // hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
443 | // if(hcalEnergy < 0.0) hcalEnergy = 0.0;
|
---|
444 |
|
---|
445 | hcalEnergy = LogNormal(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
446 |
|
---|
447 | energy = ecalEnergy + hcalEnergy;
|
---|
448 |
|
---|
449 | // eta = fTowerEta;
|
---|
450 | // phi = fTowerPhi;
|
---|
451 |
|
---|
452 | eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
|
---|
453 | phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
|
---|
454 |
|
---|
455 | pt = energy / TMath::CosH(eta);
|
---|
456 |
|
---|
457 | fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
458 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
459 | fTower->Eem = ecalEnergy;
|
---|
460 | fTower->Ehad = hcalEnergy;
|
---|
461 |
|
---|
462 | fTower->Edges[0] = fTowerEdges[0];
|
---|
463 | fTower->Edges[1] = fTowerEdges[1];
|
---|
464 | fTower->Edges[2] = fTowerEdges[2];
|
---|
465 | fTower->Edges[3] = fTowerEdges[3];
|
---|
466 |
|
---|
467 | // fill calorimeter towers and photon candidates
|
---|
468 | if(energy > 0.0)
|
---|
469 | {
|
---|
470 | if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
|
---|
471 | {
|
---|
472 | fPhotonOutputArray->Add(fTower);
|
---|
473 | }
|
---|
474 |
|
---|
475 | fTowerOutputArray->Add(fTower);
|
---|
476 | }
|
---|
477 |
|
---|
478 | // fill energy flow candidates
|
---|
479 | if(fTowerTrackAllHits == fTowerAllHits)
|
---|
480 | {
|
---|
481 | fItTowerTrackArray->Reset();
|
---|
482 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
|
---|
483 | {
|
---|
484 | fEFlowTrackOutputArray->Add(track);
|
---|
485 | }
|
---|
486 | }
|
---|
487 | else if(fTowerTrackAllHits > 0 &&
|
---|
488 | fTowerECalHits + fTowerHCalHits == fTowerAllHits &&
|
---|
489 | (fTowerECalHits == fTowerECalTrackHits || fTowerHCalHits == fTowerHCalTrackHits))
|
---|
490 | {
|
---|
491 | if(fTowerECalHits == fTowerECalTrackHits)
|
---|
492 | {
|
---|
493 | fItTowerECalTrackArray->Reset();
|
---|
494 | while((track = static_cast<Candidate*>(fItTowerECalTrackArray->Next())))
|
---|
495 | {
|
---|
496 | fEFlowTrackOutputArray->Add(track);
|
---|
497 | }
|
---|
498 | energy = hcalEnergy;
|
---|
499 | iterator = fItTowerHCalArray;
|
---|
500 | }
|
---|
501 |
|
---|
502 | if(fTowerHCalHits == fTowerHCalTrackHits)
|
---|
503 | {
|
---|
504 | fItTowerHCalTrackArray->Reset();
|
---|
505 | while((track = static_cast<Candidate*>(fItTowerHCalTrackArray->Next())))
|
---|
506 | {
|
---|
507 | fEFlowTrackOutputArray->Add(track);
|
---|
508 | }
|
---|
509 | energy = ecalEnergy;
|
---|
510 | iterator = fItTowerECalArray;
|
---|
511 | }
|
---|
512 |
|
---|
513 | if(fTowerECalHits == fTowerECalTrackHits &&
|
---|
514 | fTowerHCalHits == fTowerHCalTrackHits)
|
---|
515 | {
|
---|
516 | energy = 0.0;
|
---|
517 | iterator = 0;
|
---|
518 | }
|
---|
519 |
|
---|
520 | if(energy > 0.0 && iterator)
|
---|
521 | {
|
---|
522 | DelphesFactory *factory = GetFactory();
|
---|
523 |
|
---|
524 | // create new tower
|
---|
525 | tower = factory->NewCandidate();
|
---|
526 |
|
---|
527 | iterator->Reset();
|
---|
528 | while((particle = static_cast<Candidate*>(iterator->Next())))
|
---|
529 | {
|
---|
530 | tower->AddCandidate(particle);
|
---|
531 | }
|
---|
532 |
|
---|
533 | pt = energy / TMath::CosH(eta);
|
---|
534 |
|
---|
535 | tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
536 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
537 | tower->Eem = energy;
|
---|
538 | tower->Ehad = 0.0;
|
---|
539 |
|
---|
540 | tower->Edges[0] = fTowerEdges[0];
|
---|
541 | tower->Edges[1] = fTowerEdges[1];
|
---|
542 | tower->Edges[2] = fTowerEdges[2];
|
---|
543 | tower->Edges[3] = fTowerEdges[3];
|
---|
544 |
|
---|
545 | fEFlowTowerOutputArray->Add(tower);
|
---|
546 | }
|
---|
547 | }
|
---|
548 | else if(energy > 0.0)
|
---|
549 | {
|
---|
550 | fEFlowTowerOutputArray->Add(fTower);
|
---|
551 | }
|
---|
552 | }
|
---|
553 |
|
---|
554 | //------------------------------------------------------------------------------
|
---|
555 |
|
---|
556 | Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
|
---|
557 | {
|
---|
558 | Double_t a, b;
|
---|
559 |
|
---|
560 | if(mean > 0.0)
|
---|
561 | {
|
---|
562 | b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
|
---|
563 | a = TMath::Log(mean) - 0.5*b*b;
|
---|
564 |
|
---|
565 | return TMath::Exp(a + b*gRandom->Gaus(0, 1));
|
---|
566 | }
|
---|
567 | else
|
---|
568 | {
|
---|
569 | return 0.0;
|
---|
570 | }
|
---|
571 | }
|
---|