Fork me on GitHub

source: git/modules/Calorimeter.cc@ cdeea24

ImprovedOutputFile Timing dual_readout llp
Last change on this file since cdeea24 was cdeea24, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

fix track momentum scaling in Calorimeter
  • Property mode set to 100644
File size: 19.3 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
435 continue;
436 }
437
438 // check for photon and electron hits in current tower
439 if(flags & 2) ++fTowerPhotonHits;
440
441 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
442 momentum = particle->Momentum;
443 position = particle->Position;
444
445 // fill current tower
446 ecalEnergy = momentum.E() * fECalTowerFractions[number];
447 hcalEnergy = momentum.E() * fHCalTowerFractions[number];
448
449 fECalTowerEnergy += ecalEnergy;
450 fHCalTowerEnergy += hcalEnergy;
451
452 if(ecalEnergy > fTimingEnergyMin && fTower)
453 {
454 if (abs(particle->PID) != 11 || !fElectronsFromTrack)
455 {
456 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
457 }
458 }
459
460 fTower->AddCandidate(particle);
461 }
462
463 // finalize last tower
464 FinalizeTower();
465}
466
467//------------------------------------------------------------------------------
468
469void Calorimeter::FinalizeTower()
470{
471 Candidate *track, *tower, *mother;
472 Double_t energy, pt, eta, phi;
473 Double_t ecalEnergy, hcalEnergy;
474 Double_t ecalSigma, hcalSigma;
475
476 TLorentzVector momentum;
477 TFractionMap::iterator itFractionMap;
478
479 Float_t weight, sumWeightedTime, sumWeight;
480
481 if(!fTower) return;
482
483
484 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
485 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
486
487 ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
488 hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
489
490 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
491 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
492
493 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
494 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
495
496 energy = ecalEnergy + hcalEnergy;
497
498 if(fSmearTowerCenter)
499 {
500 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
501 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
502 }
503 else
504 {
505 eta = fTowerEta;
506 phi = fTowerPhi;
507 }
508
509 pt = energy / TMath::CosH(eta);
510
511 // Time calculation for tower
512 fTower->NTimeHits = 0;
513 sumWeightedTime = 0.0;
514 sumWeight = 0.0;
515
516 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
517 {
518 weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
519 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
520 sumWeight += weight;
521 fTower->NTimeHits++;
522 }
523
524 if(sumWeight > 0.0)
525 {
526 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
527 }
528 else
529 {
530 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
531 }
532
533
534 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
535 fTower->Eem = ecalEnergy;
536 fTower->Ehad = hcalEnergy;
537
538 fTower->Edges[0] = fTowerEdges[0];
539 fTower->Edges[1] = fTowerEdges[1];
540 fTower->Edges[2] = fTowerEdges[2];
541 fTower->Edges[3] = fTowerEdges[3];
542
543 if(energy > 0.0)
544 {
545 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
546 {
547 fPhotonOutputArray->Add(fTower);
548 }
549
550 fTowerOutputArray->Add(fTower);
551 }
552
553 // fill energy flow candidates
554
555 ecalEnergy -= fECalTrackEnergy[1];
556 hcalEnergy -= fHCalTrackEnergy[1];
557
558 fItECalTowerTrackArray[0]->Reset();
559 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
560 {
561 mother = track;
562 track = static_cast<Candidate*>(track->Clone());
563 track->AddCandidate(mother);
564
565 track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
566 ecalEnergy = 0.0;
567
568 fEFlowTrackOutputArray->Add(track);
569 }
570
571 fItECalTowerTrackArray[1]->Reset();
572 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
573 {
574 mother = track;
575 track = static_cast<Candidate*>(track->Clone());
576 track->AddCandidate(mother);
577
578 fEFlowTrackOutputArray->Add(track);
579 }
580
581 fItHCalTowerTrackArray[0]->Reset();
582 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
583 {
584 mother = track;
585 track = static_cast<Candidate*>(track->Clone());
586 track->AddCandidate(mother);
587
588 track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
589 hcalEnergy = 0.0;
590
591 fEFlowTrackOutputArray->Add(track);
592 }
593
594 fItHCalTowerTrackArray[1]->Reset();
595 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
596 {
597 mother = track;
598 track = static_cast<Candidate*>(track->Clone());
599 track->AddCandidate(mother);
600
601 fEFlowTrackOutputArray->Add(track);
602 }
603
604 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
605 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
606
607 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
608 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
609
610 energy = ecalEnergy + hcalEnergy;
611
612 if(ecalEnergy > 0.0)
613 {
614 // create new photon tower
615 tower = static_cast<Candidate*>(fTower->Clone());
616
617 pt = ecalEnergy / TMath::CosH(eta);
618
619 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
620 tower->Eem = ecalEnergy;
621 tower->Ehad = 0.0;
622
623 fEFlowPhotonOutputArray->Add(tower);
624 }
625 if(hcalEnergy > 0.0)
626 {
627 // create new neutral hadron tower
628 tower = static_cast<Candidate*>(fTower->Clone());
629
630 pt = hcalEnergy / TMath::CosH(eta);
631
632 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
633 tower->Eem = 0.0;
634 tower->Ehad = hcalEnergy;
635
636 fEFlowNeutralHadronOutputArray->Add(tower);
637 }
638}
639
640//------------------------------------------------------------------------------
641
642Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
643{
644 Double_t a, b;
645
646 if(mean > 0.0)
647 {
648 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
649 a = TMath::Log(mean) - 0.5*b*b;
650
651 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
652 }
653 else
654 {
655 return 0.0;
656 }
657}
Note: See TracBrowser for help on using the repository browser.