Fork me on GitHub

source: git/modules/Calorimeter.cc@ 40e890c

Last change on this file since 40e890c was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

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