Fork me on GitHub

source: git/modules/Calorimeter.cc@ 2628d2b

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 2628d2b was f8299bc, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

added PID to EFlowPhotons
  • Property mode set to 100644
File size: 19.5 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class Calorimeter
21 *
22 * Fills calorimeter towers, performs calorimeter resolution smearing,
23 * and creates energy flow objects (tracks, photons, and neutral hadrons).
24 *
25 * \author P. Demin - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/Calorimeter.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56Calorimeter::Calorimeter() :
57 fECalResolutionFormula(0), fHCalResolutionFormula(0),
58 fItParticleInputArray(0), fItTrackInputArray(0)
59{
60 Int_t i;
61
62 fECalResolutionFormula = new DelphesFormula;
63 fHCalResolutionFormula = new DelphesFormula;
64
65 for(i = 0; i < 2; ++i)
66 {
67 fECalTowerTrackArray[i] = new TObjArray;
68 fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
69
70 fHCalTowerTrackArray[i] = new TObjArray;
71 fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
72 }
73}
74
75//------------------------------------------------------------------------------
76
77Calorimeter::~Calorimeter()
78{
79 Int_t i;
80
81 if(fECalResolutionFormula) delete fECalResolutionFormula;
82 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
83
84 for(i = 0; i < 2; ++i)
85 {
86 if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
87 if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
88
89 if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
90 if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
91 }
92}
93
94//------------------------------------------------------------------------------
95
96void Calorimeter::Init()
97{
98 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
99 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
100 Double_t ecalFraction, hcalFraction;
101 TBinMap::iterator itEtaBin;
102 set< Double_t >::iterator itPhiBin;
103 vector< Double_t > *phiBins;
104
105 // read eta and phi bins
106 param = GetParam("EtaPhiBins");
107 size = param.GetSize();
108 fBinMap.clear();
109 fEtaBins.clear();
110 fPhiBins.clear();
111 for(i = 0; i < size/2; ++i)
112 {
113 paramEtaBins = param[i*2];
114 sizeEtaBins = paramEtaBins.GetSize();
115 paramPhiBins = param[i*2 + 1];
116 sizePhiBins = paramPhiBins.GetSize();
117
118 for(j = 0; j < sizeEtaBins; ++j)
119 {
120 for(k = 0; k < sizePhiBins; ++k)
121 {
122 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
123 }
124 }
125 }
126
127 // for better performance we transform map of sets to parallel vectors:
128 // vector< double > and vector< vector< double >* >
129 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
130 {
131 fEtaBins.push_back(itEtaBin->first);
132 phiBins = new vector< double >(itEtaBin->second.size());
133 fPhiBins.push_back(phiBins);
134 phiBins->clear();
135 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
136 {
137 phiBins->push_back(*itPhiBin);
138 }
139 }
140
141 // read energy fractions for different particles
142 param = GetParam("EnergyFraction");
143 size = param.GetSize();
144
145 // set default energy fractions values
146 fFractionMap.clear();
147 fFractionMap[0] = make_pair(0.0, 1.0);
148
149 for(i = 0; i < size/2; ++i)
150 {
151 paramFractions = param[i*2 + 1];
152
153 ecalFraction = paramFractions[0].GetDouble();
154 hcalFraction = paramFractions[1].GetDouble();
155
156 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
157 }
158
159 // read min E value for timing measurement in ECAL
160 fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
161 // For timing
162 // So far this flag needs to be false
163 // Curved extrapolation not supported
164 fElectronsFromTrack = false;
165
166 // read min E value for towers to be saved
167 fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
168 fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
169
170 fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
171 fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
172
173 // switch on or off the dithering of the center of calorimeter towers
174 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
175
176 // read resolution formulas
177 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
178 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
179
180 // import array with output from other modules
181 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
182 fItParticleInputArray = fParticleInputArray->MakeIterator();
183
184 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
185 fItTrackInputArray = fTrackInputArray->MakeIterator();
186
187 // create output arrays
188 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
189 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
190
191 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
192 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
193 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
194}
195
196//------------------------------------------------------------------------------
197
198void Calorimeter::Finish()
199{
200 vector< vector< Double_t >* >::iterator itPhiBin;
201 if(fItParticleInputArray) delete fItParticleInputArray;
202 if(fItTrackInputArray) delete fItTrackInputArray;
203 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
204 {
205 delete *itPhiBin;
206 }
207}
208
209//------------------------------------------------------------------------------
210
211void Calorimeter::Process()
212{
213 Candidate *particle, *track;
214 TLorentzVector position, momentum;
215 Short_t etaBin, phiBin, flags;
216 Int_t number;
217 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
218 Double_t ecalFraction, hcalFraction;
219 Double_t ecalEnergy, hcalEnergy;
220 Double_t ecalSigma, hcalSigma;
221 Int_t pdgCode;
222
223 TFractionMap::iterator itFractionMap;
224
225 vector< Double_t >::iterator itEtaBin;
226 vector< Double_t >::iterator itPhiBin;
227 vector< Double_t > *phiBins;
228
229 vector< Long64_t >::iterator itTowerHits;
230
231 DelphesFactory *factory = GetFactory();
232 fTowerHits.clear();
233 fECalTowerFractions.clear();
234 fHCalTowerFractions.clear();
235 fECalTrackFractions.clear();
236 fHCalTrackFractions.clear();
237
238 // loop over all particles
239 fItParticleInputArray->Reset();
240 number = -1;
241 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
242 {
243 const TLorentzVector &particlePosition = particle->Position;
244 ++number;
245
246 pdgCode = TMath::Abs(particle->PID);
247
248 itFractionMap = fFractionMap.find(pdgCode);
249 if(itFractionMap == fFractionMap.end())
250 {
251 itFractionMap = fFractionMap.find(0);
252 }
253
254 ecalFraction = itFractionMap->second.first;
255 hcalFraction = itFractionMap->second.second;
256
257 fECalTowerFractions.push_back(ecalFraction);
258 fHCalTowerFractions.push_back(hcalFraction);
259
260 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
261
262 // find eta bin [1, fEtaBins.size - 1]
263 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
264 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
265 etaBin = distance(fEtaBins.begin(), itEtaBin);
266
267 // phi bins for given eta bin
268 phiBins = fPhiBins[etaBin];
269
270 // find phi bin [1, phiBins.size - 1]
271 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
272 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
273 phiBin = distance(phiBins->begin(), itPhiBin);
274
275 flags = 0;
276 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
277
278 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
279 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
280
281 fTowerHits.push_back(towerHit);
282 }
283
284 // loop over all tracks
285 fItTrackInputArray->Reset();
286 number = -1;
287 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
288 {
289 const TLorentzVector &trackPosition = track->Position;
290 ++number;
291
292 pdgCode = TMath::Abs(track->PID);
293
294 itFractionMap = fFractionMap.find(pdgCode);
295 if(itFractionMap == fFractionMap.end())
296 {
297 itFractionMap = fFractionMap.find(0);
298 }
299
300 ecalFraction = itFractionMap->second.first;
301 hcalFraction = itFractionMap->second.second;
302
303 fECalTrackFractions.push_back(ecalFraction);
304 fHCalTrackFractions.push_back(hcalFraction);
305
306 // find eta bin [1, fEtaBins.size - 1]
307 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
308 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
309 etaBin = distance(fEtaBins.begin(), itEtaBin);
310
311 // phi bins for given eta bin
312 phiBins = fPhiBins[etaBin];
313
314 // find phi bin [1, phiBins.size - 1]
315 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
316 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
317 phiBin = distance(phiBins->begin(), itPhiBin);
318
319 flags = 1;
320
321 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
322 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
323
324 fTowerHits.push_back(towerHit);
325 }
326
327 // all hits are sorted first by eta bin number, then by phi bin number,
328 // then by flags and then by particle or track number
329 sort(fTowerHits.begin(), fTowerHits.end());
330
331 // loop over all hits
332 towerEtaPhi = 0;
333 fTower = 0;
334 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
335 {
336 towerHit = (*itTowerHits);
337 flags = (towerHit >> 24) & 0x00000000000000FFLL;
338 number = (towerHit) & 0x0000000000FFFFFFLL;
339 hitEtaPhi = towerHit >> 32;
340
341 if(towerEtaPhi != hitEtaPhi)
342 {
343 // switch to next tower
344 towerEtaPhi = hitEtaPhi;
345
346 // finalize previous tower
347 FinalizeTower();
348
349 // create new tower
350 fTower = factory->NewCandidate();
351
352 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
353 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
354
355 // phi bins for given eta bin
356 phiBins = fPhiBins[etaBin];
357
358 // calculate eta and phi of the tower's center
359 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
360 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
361
362 fTowerEdges[0] = fEtaBins[etaBin - 1];
363 fTowerEdges[1] = fEtaBins[etaBin];
364 fTowerEdges[2] = (*phiBins)[phiBin - 1];
365 fTowerEdges[3] = (*phiBins)[phiBin];
366
367 fECalTowerEnergy = 0.0;
368 fHCalTowerEnergy = 0.0;
369
370 fECalTrackEnergy[0] = 0.0;
371 fECalTrackEnergy[1] = 0.0;
372
373 fHCalTrackEnergy[0] = 0.0;
374 fHCalTrackEnergy[1] = 0.0;
375
376 fTowerTrackHits = 0;
377 fTowerPhotonHits = 0;
378
379 fECalTowerTrackArray[0]->Clear();
380 fECalTowerTrackArray[1]->Clear();
381
382 fHCalTowerTrackArray[0]->Clear();
383 fHCalTowerTrackArray[1]->Clear();
384 }
385
386 // check for track hits
387 if(flags & 1)
388 {
389 ++fTowerTrackHits;
390
391 track = static_cast<Candidate*>(fTrackInputArray->At(number));
392 momentum = track->Momentum;
393 position = track->Position;
394
395 ecalEnergy = momentum.E() * fECalTrackFractions[number];
396 hcalEnergy = momentum.E() * fHCalTrackFractions[number];
397
398 if(ecalEnergy > fTimingEnergyMin && fTower)
399 {
400 if(fElectronsFromTrack)
401 {
402 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
403 }
404 }
405
406 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
407 {
408 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
409 if(ecalSigma/momentum.E() < track->TrackResolution)
410 {
411 fECalTrackEnergy[0] += ecalEnergy;
412 fECalTowerTrackArray[0]->Add(track);
413 }
414 else
415 {
416 fECalTrackEnergy[1] += ecalEnergy;
417 fECalTowerTrackArray[1]->Add(track);
418 }
419 }
420 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
421 {
422 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
423 if(hcalSigma/momentum.E() < track->TrackResolution)
424 {
425 fHCalTrackEnergy[0] += hcalEnergy;
426 fHCalTowerTrackArray[0]->Add(track);
427 }
428 else
429 {
430 fHCalTrackEnergy[1] += hcalEnergy;
431 fHCalTowerTrackArray[1]->Add(track);
432 }
433 }
434 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
435 {
436 fEFlowTrackOutputArray->Add(track);
437 }
438
439 continue;
440 }
441
442 // check for photon and electron hits in current tower
443 if(flags & 2) ++fTowerPhotonHits;
444
445 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
446 momentum = particle->Momentum;
447 position = particle->Position;
448
449 // fill current tower
450 ecalEnergy = momentum.E() * fECalTowerFractions[number];
451 hcalEnergy = momentum.E() * fHCalTowerFractions[number];
452
453 fECalTowerEnergy += ecalEnergy;
454 fHCalTowerEnergy += hcalEnergy;
455
456 if(ecalEnergy > fTimingEnergyMin && fTower)
457 {
458 if (abs(particle->PID) != 11 || !fElectronsFromTrack)
459 {
460 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
461 }
462 }
463
464 fTower->AddCandidate(particle);
465 }
466
467 // finalize last tower
468 FinalizeTower();
469}
470
471//------------------------------------------------------------------------------
472
473void Calorimeter::FinalizeTower()
474{
475 Candidate *track, *tower, *mother;
476 Double_t energy, pt, eta, phi;
477 Double_t ecalEnergy, hcalEnergy;
478 Double_t ecalSigma, hcalSigma;
479
480 TLorentzVector momentum;
481 TFractionMap::iterator itFractionMap;
482
483 Float_t weight, sumWeightedTime, sumWeight;
484
485 if(!fTower) return;
486
487
488 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
489 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
490
491 ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
492 hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
493
494 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
495 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
496
497 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
498 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
499
500 energy = ecalEnergy + hcalEnergy;
501
502 if(fSmearTowerCenter)
503 {
504 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
505 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
506 }
507 else
508 {
509 eta = fTowerEta;
510 phi = fTowerPhi;
511 }
512
513 pt = energy / TMath::CosH(eta);
514
515 // Time calculation for tower
516 fTower->NTimeHits = 0;
517 sumWeightedTime = 0.0;
518 sumWeight = 0.0;
519
520 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
521 {
522 weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
523 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
524 sumWeight += weight;
525 fTower->NTimeHits++;
526 }
527
528 if(sumWeight > 0.0)
529 {
530 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
531 }
532 else
533 {
534 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
535 }
536
537
538 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
539 fTower->Eem = ecalEnergy;
540 fTower->Ehad = hcalEnergy;
541
542 fTower->Edges[0] = fTowerEdges[0];
543 fTower->Edges[1] = fTowerEdges[1];
544 fTower->Edges[2] = fTowerEdges[2];
545 fTower->Edges[3] = fTowerEdges[3];
546
547 if(energy > 0.0)
548 {
549 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
550 {
551 fPhotonOutputArray->Add(fTower);
552 }
553
554 fTowerOutputArray->Add(fTower);
555 }
556
557 // fill energy flow candidates
558
559 ecalEnergy -= fECalTrackEnergy[1];
560 hcalEnergy -= fHCalTrackEnergy[1];
561
562 fItECalTowerTrackArray[0]->Reset();
563 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
564 {
565 mother = track;
566 track = static_cast<Candidate*>(track->Clone());
567 track->AddCandidate(mother);
568
569 track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
570
571 fEFlowTrackOutputArray->Add(track);
572 }
573
574 fItECalTowerTrackArray[1]->Reset();
575 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
576 {
577 mother = track;
578 track = static_cast<Candidate*>(track->Clone());
579 track->AddCandidate(mother);
580
581 fEFlowTrackOutputArray->Add(track);
582 }
583
584 fItHCalTowerTrackArray[0]->Reset();
585 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
586 {
587 mother = track;
588 track = static_cast<Candidate*>(track->Clone());
589 track->AddCandidate(mother);
590
591 track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
592
593 fEFlowTrackOutputArray->Add(track);
594 }
595
596 fItHCalTowerTrackArray[1]->Reset();
597 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
598 {
599 mother = track;
600 track = static_cast<Candidate*>(track->Clone());
601 track->AddCandidate(mother);
602
603 fEFlowTrackOutputArray->Add(track);
604 }
605
606 if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
607 if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0;
608
609 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
610 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
611
612 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
613 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
614
615 energy = ecalEnergy + hcalEnergy;
616
617 if(ecalEnergy > 0.0)
618 {
619 // create new photon tower
620 tower = static_cast<Candidate*>(fTower->Clone());
621
622 pt = ecalEnergy / TMath::CosH(eta);
623
624 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
625 tower->Eem = ecalEnergy;
626 tower->Ehad = 0.0;
627 tower->PID = 22;
628
629 fEFlowPhotonOutputArray->Add(tower);
630 }
631 if(hcalEnergy > 0.0)
632 {
633 // create new neutral hadron tower
634 tower = static_cast<Candidate*>(fTower->Clone());
635
636 pt = hcalEnergy / TMath::CosH(eta);
637
638 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
639 tower->Eem = 0.0;
640 tower->Ehad = hcalEnergy;
641
642 fEFlowNeutralHadronOutputArray->Add(tower);
643 }
644}
645
646//------------------------------------------------------------------------------
647
648Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
649{
650 Double_t a, b;
651
652 if(mean > 0.0)
653 {
654 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
655 a = TMath::Log(mean) - 0.5*b*b;
656
657 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
658 }
659 else
660 {
661 return 0.0;
662 }
663}
Note: See TracBrowser for help on using the repository browser.