Fork me on GitHub

source: git/modules/Calorimeter.cc@ 78e55cf

Last change on this file since 78e55cf was 5eda6767, checked in by Michele Selvaggi <michele@…>, 4 years ago

fix calo tower position

  • Property mode set to 100644
File size: 21.6 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 fTowerRmax=0.;
234 while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
235 {
236 const TLorentzVector &particlePosition = particle->Position;
237 ++number;
238
239 // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
240 if (particlePosition.Perp() > fTowerRmax)
241 fTowerRmax=particlePosition.Perp();
242
243 pdgCode = TMath::Abs(particle->PID);
244
245 itFractionMap = fFractionMap.find(pdgCode);
246 if(itFractionMap == fFractionMap.end())
247 {
248 itFractionMap = fFractionMap.find(0);
249 }
250
251 ecalFraction = itFractionMap->second.first;
252 hcalFraction = itFractionMap->second.second;
253
254 fECalTowerFractions.push_back(ecalFraction);
255 fHCalTowerFractions.push_back(hcalFraction);
256
257 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
258
259 // find eta bin [1, fEtaBins.size - 1]
260 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
261 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
262 etaBin = distance(fEtaBins.begin(), itEtaBin);
263
264 // phi bins for given eta bin
265 phiBins = fPhiBins[etaBin];
266
267 // find phi bin [1, phiBins.size - 1]
268 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
269 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
270 phiBin = distance(phiBins->begin(), itPhiBin);
271
272 flags = 0;
273 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
274
275 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
276 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
277
278 fTowerHits.push_back(towerHit);
279 }
280
281 // loop over all tracks
282 fItTrackInputArray->Reset();
283 number = -1;
284 while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
285 {
286 const TLorentzVector &trackPosition = track->Position;
287 ++number;
288
289 pdgCode = TMath::Abs(track->PID);
290
291 itFractionMap = fFractionMap.find(pdgCode);
292 if(itFractionMap == fFractionMap.end())
293 {
294 itFractionMap = fFractionMap.find(0);
295 }
296
297 ecalFraction = itFractionMap->second.first;
298 hcalFraction = itFractionMap->second.second;
299
300 fECalTrackFractions.push_back(ecalFraction);
301 fHCalTrackFractions.push_back(hcalFraction);
302
303 // find eta bin [1, fEtaBins.size - 1]
304 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
305 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
306 etaBin = distance(fEtaBins.begin(), itEtaBin);
307
308 // phi bins for given eta bin
309 phiBins = fPhiBins[etaBin];
310
311 // find phi bin [1, phiBins.size - 1]
312 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
313 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
314 phiBin = distance(phiBins->begin(), itPhiBin);
315
316 flags = 1;
317
318 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
319 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
320
321 fTowerHits.push_back(towerHit);
322 }
323
324 // all hits are sorted first by eta bin number, then by phi bin number,
325 // then by flags and then by particle or track number
326 sort(fTowerHits.begin(), fTowerHits.end());
327
328 // loop over all hits
329 towerEtaPhi = 0;
330 fTower = 0;
331 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
332 {
333 towerHit = (*itTowerHits);
334 flags = (towerHit >> 24) & 0x00000000000000FFLL;
335 number = (towerHit)&0x0000000000FFFFFFLL;
336 hitEtaPhi = towerHit >> 32;
337
338 if(towerEtaPhi != hitEtaPhi)
339 {
340 // switch to next tower
341 towerEtaPhi = hitEtaPhi;
342
343 // finalize previous tower
344 FinalizeTower();
345
346 // create new tower
347 fTower = factory->NewCandidate();
348
349 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
350 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
351
352 // phi bins for given eta bin
353 phiBins = fPhiBins[etaBin];
354
355 // calculate eta and phi of the tower's center
356 fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
357 fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
358
359 fTowerEdges[0] = fEtaBins[etaBin - 1];
360 fTowerEdges[1] = fEtaBins[etaBin];
361 fTowerEdges[2] = (*phiBins)[phiBin - 1];
362 fTowerEdges[3] = (*phiBins)[phiBin];
363
364 fECalTowerEnergy = 0.0;
365 fHCalTowerEnergy = 0.0;
366
367 fECalTrackEnergy = 0.0;
368 fHCalTrackEnergy = 0.0;
369
370 fECalTrackSigma = 0.0;
371 fHCalTrackSigma = 0.0;
372
373 fTowerTrackHits = 0;
374 fTowerPhotonHits = 0;
375
376 fECalTowerTrackArray->Clear();
377 fHCalTowerTrackArray->Clear();
378 }
379
380 // check for track hits
381 if(flags & 1)
382 {
383 ++fTowerTrackHits;
384
385 track = static_cast<Candidate *>(fTrackInputArray->At(number));
386 momentum = track->Momentum;
387 position = track->Position;
388
389 ecalEnergy = momentum.E() * fECalTrackFractions[number];
390 hcalEnergy = momentum.E() * fHCalTrackFractions[number];
391
392 if(ecalEnergy > fTimingEnergyMin && fTower)
393 {
394 if(fElectronsFromTrack)
395 {
396 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
397 }
398 }
399
400 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
401 {
402 fECalTrackEnergy += ecalEnergy;
403 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
404 if(ecalSigma / momentum.E() < track->TrackResolution)
405 energyGuess = ecalEnergy;
406 else
407 energyGuess = momentum.E();
408
409 fECalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
410 fECalTowerTrackArray->Add(track);
411 }
412
413 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
414 {
415 fHCalTrackEnergy += hcalEnergy;
416 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
417 if(hcalSigma / momentum.E() < track->TrackResolution)
418 energyGuess = hcalEnergy;
419 else
420 energyGuess = momentum.E();
421
422 fHCalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
423 fHCalTowerTrackArray->Add(track);
424 }
425
426 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
427 {
428 fEFlowTrackOutputArray->Add(track);
429 }
430
431 continue;
432 }
433
434 // check for photon and electron hits in current tower
435 if(flags & 2) ++fTowerPhotonHits;
436
437 particle = static_cast<Candidate *>(fParticleInputArray->At(number));
438 momentum = particle->Momentum;
439 position = particle->Position;
440
441 // fill current tower
442 ecalEnergy = momentum.E() * fECalTowerFractions[number];
443 hcalEnergy = momentum.E() * fHCalTowerFractions[number];
444
445 fECalTowerEnergy += ecalEnergy;
446 fHCalTowerEnergy += hcalEnergy;
447
448 if(ecalEnergy > fTimingEnergyMin && fTower)
449 {
450 if(abs(particle->PID) != 11 || !fElectronsFromTrack)
451 {
452 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
453 }
454 }
455
456 fTower->AddCandidate(particle);
457 fTower->Position = position;
458 }
459
460 // finalize last tower
461 FinalizeTower();
462}
463
464//------------------------------------------------------------------------------
465
466void Calorimeter::FinalizeTower()
467{
468 Candidate *track, *tower, *mother;
469 Double_t energy, pt, eta, phi, r;
470 Double_t ecalEnergy, hcalEnergy;
471 Double_t ecalNeutralEnergy, hcalNeutralEnergy;
472
473 Double_t ecalSigma, hcalSigma;
474 Double_t ecalNeutralSigma, hcalNeutralSigma;
475
476 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
477
478 TLorentzVector momentum;
479 TFractionMap::iterator itFractionMap;
480
481 Float_t weight, sumWeightedTime, sumWeight;
482
483 if(!fTower) return;
484
485 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
486 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
487
488 ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
489 hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
490
491 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
492 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
493
494 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin * ecalSigma) ecalEnergy = 0.0;
495 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin * hcalSigma) hcalEnergy = 0.0;
496
497 energy = ecalEnergy + hcalEnergy;
498
499 if(fSmearTowerCenter)
500 {
501 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
502 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
503 }
504 else
505 {
506 eta = fTowerEta;
507 phi = fTowerPhi;
508 }
509
510 pt = energy / TMath::CosH(eta);
511
512 // Time calculation for tower
513 fTower->NTimeHits = 0;
514 sumWeightedTime = 0.0;
515 sumWeight = 0.0;
516
517 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
518 {
519 weight = TMath::Power((fTower->ECalEnergyTimePairs[i].first),2);
520 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
521 sumWeight += weight;
522 fTower->NTimeHits++;
523 }
524
525 // check whether barrel or endcap tower
526 if (fTower->Position.Perp() < fTowerRmax && TMath::Abs(eta) > 0.)
527 r = fTower->Position.Z()/TMath::SinH(eta);
528 else
529 r = fTower->Position.Pt();
530
531 if(sumWeight > 0.0)
532 {
533 fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime / sumWeight);
534 }
535 else
536 {
537 fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9);
538 }
539
540 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
541 fTower->Eem = ecalEnergy;
542 fTower->Ehad = hcalEnergy;
543
544 fTower->Edges[0] = fTowerEdges[0];
545 fTower->Edges[1] = fTowerEdges[1];
546 fTower->Edges[2] = fTowerEdges[2];
547 fTower->Edges[3] = fTowerEdges[3];
548
549 if(energy > 0.0)
550 {
551 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
552 {
553 fPhotonOutputArray->Add(fTower);
554 }
555
556 fTowerOutputArray->Add(fTower);
557 }
558
559 // fill energy flow candidates
560 fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
561 fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
562
563 //compute neutral excesses
564 ecalNeutralEnergy = max((ecalEnergy - fECalTrackEnergy), 0.0);
565 hcalNeutralEnergy = max((hcalEnergy - fHCalTrackEnergy), 0.0);
566
567 ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma * fECalTrackSigma + ecalSigma * ecalSigma);
568 hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma * fHCalTrackSigma + hcalSigma * hcalSigma);
569
570 // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
571 if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
572 {
573 // create new photon tower assuming null mass
574 tower = static_cast<Candidate *>(fTower->Clone());
575 pt = ecalNeutralEnergy / TMath::CosH(eta);
576
577 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
578 tower->Eem = ecalNeutralEnergy;
579 tower->Ehad = 0.0;
580 tower->PID = 22;
581
582 fEFlowPhotonOutputArray->Add(tower);
583
584 //clone tracks
585 fItECalTowerTrackArray->Reset();
586 while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
587 {
588 mother = track;
589 track = static_cast<Candidate *>(track->Clone());
590 track->AddCandidate(mother);
591
592 fEFlowTrackOutputArray->Add(track);
593 }
594 }
595
596 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
597 else if(fECalTrackEnergy > 0.0)
598 {
599 weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma * fECalTrackSigma) : 0.0;
600 weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma * ecalSigma) : 0.0;
601
602 bestEnergyEstimate = (weightTrack * fECalTrackEnergy + weightCalo * ecalEnergy) / (weightTrack + weightCalo);
603 rescaleFactor = bestEnergyEstimate / fECalTrackEnergy;
604
605 //rescale tracks
606 fItECalTowerTrackArray->Reset();
607 while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
608 {
609 mother = track;
610 track = static_cast<Candidate *>(track->Clone());
611 track->AddCandidate(mother);
612
613 track->Momentum *= rescaleFactor;
614
615 fEFlowTrackOutputArray->Add(track);
616 }
617 }
618
619 // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
620 if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
621 {
622 // create new photon tower
623 tower = static_cast<Candidate *>(fTower->Clone());
624 pt = hcalNeutralEnergy / TMath::CosH(eta);
625
626 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
627 tower->Ehad = hcalNeutralEnergy;
628 tower->Eem = 0.0;
629
630 fEFlowNeutralHadronOutputArray->Add(tower);
631
632 //clone tracks
633 fItHCalTowerTrackArray->Reset();
634 while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
635 {
636 mother = track;
637 track = static_cast<Candidate *>(track->Clone());
638 track->AddCandidate(mother);
639
640 fEFlowTrackOutputArray->Add(track);
641 }
642 }
643
644 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
645 else if(fHCalTrackEnergy > 0.0)
646 {
647 weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma * fHCalTrackSigma) : 0.0;
648 weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma * hcalSigma) : 0.0;
649
650 bestEnergyEstimate = (weightTrack * fHCalTrackEnergy + weightCalo * hcalEnergy) / (weightTrack + weightCalo);
651 rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
652
653 //rescale tracks
654 fItHCalTowerTrackArray->Reset();
655 while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
656 {
657 mother = track;
658 track = static_cast<Candidate *>(track->Clone());
659 track->AddCandidate(mother);
660 track->Momentum *= rescaleFactor;
661 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
662
663 fEFlowTrackOutputArray->Add(track);
664 }
665 }
666}
667
668//------------------------------------------------------------------------------
669
670Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
671{
672 Double_t a, b;
673
674 if(mean > 0.0)
675 {
676 b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
677 a = TMath::Log(mean) - 0.5 * b * b;
678
679 return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0));
680 }
681 else
682 {
683 return 0.0;
684 }
685}
Note: See TracBrowser for help on using the repository browser.