Fork me on GitHub

source: git/modules/DualReadoutCalorimeter.cc@ 952bbbc

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

fix calo tower position

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