[608] | 1 |
|
---|
[814] | 2 | /** \class Calorimeter
|
---|
| 3 | *
|
---|
[894] | 4 | * Fills calorimeter towers, performs calorimeter resolution smearing,
|
---|
| 5 | * preselects towers hit by photons and creates energy flow objects.
|
---|
[814] | 6 | *
|
---|
| 7 | * $Date: 2013-09-03 19:41:25 +0000 (Tue, 03 Sep 2013) $
|
---|
| 8 | * $Revision: 1274 $
|
---|
| 9 | *
|
---|
| 10 | *
|
---|
| 11 | * \author P. Demin - UCL, Louvain-la-Neuve
|
---|
| 12 | *
|
---|
| 13 | */
|
---|
| 14 |
|
---|
[608] | 15 | #include "modules/Calorimeter.h"
|
---|
| 16 |
|
---|
[687] | 17 | #include "classes/DelphesClasses.h"
|
---|
| 18 | #include "classes/DelphesFactory.h"
|
---|
[766] | 19 | #include "classes/DelphesFormula.h"
|
---|
[608] | 20 |
|
---|
[703] | 21 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
| 22 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
| 23 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
| 24 |
|
---|
[608] | 25 | #include "TMath.h"
|
---|
| 26 | #include "TString.h"
|
---|
[629] | 27 | #include "TFormula.h"
|
---|
[703] | 28 | #include "TRandom3.h"
|
---|
| 29 | #include "TObjArray.h"
|
---|
| 30 | #include "TDatabasePDG.h"
|
---|
| 31 | #include "TLorentzVector.h"
|
---|
[608] | 32 |
|
---|
[935] | 33 | #include <algorithm>
|
---|
[703] | 34 | #include <stdexcept>
|
---|
[608] | 35 | #include <iostream>
|
---|
[703] | 36 | #include <sstream>
|
---|
[608] | 37 |
|
---|
| 38 | using namespace std;
|
---|
| 39 |
|
---|
| 40 | //------------------------------------------------------------------------------
|
---|
| 41 |
|
---|
| 42 | Calorimeter::Calorimeter() :
|
---|
[886] | 43 | fECalResolutionFormula(0), fHCalResolutionFormula(0),
|
---|
| 44 | fItParticleInputArray(0), fItTrackInputArray(0),
|
---|
[1273] | 45 | fTowerTrackArray(0), fItTowerTrackArray(0)
|
---|
[608] | 46 | {
|
---|
[766] | 47 | fECalResolutionFormula = new DelphesFormula;
|
---|
| 48 | fHCalResolutionFormula = new DelphesFormula;
|
---|
[1233] | 49 |
|
---|
[1235] | 50 | fTowerTrackArray = new TObjArray;
|
---|
| 51 | fItTowerTrackArray = fTowerTrackArray->MakeIterator();
|
---|
[608] | 52 | }
|
---|
| 53 |
|
---|
| 54 | //------------------------------------------------------------------------------
|
---|
| 55 |
|
---|
| 56 | Calorimeter::~Calorimeter()
|
---|
| 57 | {
|
---|
[766] | 58 | if(fECalResolutionFormula) delete fECalResolutionFormula;
|
---|
| 59 | if(fHCalResolutionFormula) delete fHCalResolutionFormula;
|
---|
[608] | 60 |
|
---|
[1235] | 61 | if(fTowerTrackArray) delete fTowerTrackArray;
|
---|
| 62 | if(fItTowerTrackArray) delete fItTowerTrackArray;
|
---|
[1273] | 63 | }
|
---|
[1233] | 64 |
|
---|
[608] | 65 | //------------------------------------------------------------------------------
|
---|
| 66 |
|
---|
| 67 | void Calorimeter::Init()
|
---|
| 68 | {
|
---|
[704] | 69 | ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
|
---|
| 70 | Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
|
---|
[608] | 71 | Double_t ecalFraction, hcalFraction;
|
---|
[622] | 72 | TBinMap::iterator itEtaBin;
|
---|
[659] | 73 | set< Double_t >::iterator itPhiBin;
|
---|
| 74 | vector< Double_t > *phiBins;
|
---|
[608] | 75 |
|
---|
[613] | 76 | // read eta and phi bins
|
---|
[646] | 77 | param = GetParam("EtaPhiBins");
|
---|
[608] | 78 | size = param.GetSize();
|
---|
[622] | 79 | fBinMap.clear();
|
---|
[663] | 80 | fEtaBins.clear();
|
---|
| 81 | fPhiBins.clear();
|
---|
[613] | 82 | for(i = 0; i < size/2; ++i)
|
---|
[608] | 83 | {
|
---|
[646] | 84 | paramEtaBins = param[i*2];
|
---|
| 85 | sizeEtaBins = paramEtaBins.GetSize();
|
---|
[613] | 86 | paramPhiBins = param[i*2 + 1];
|
---|
| 87 | sizePhiBins = paramPhiBins.GetSize();
|
---|
[646] | 88 |
|
---|
| 89 | for(j = 0; j < sizeEtaBins; ++j)
|
---|
[613] | 90 | {
|
---|
[646] | 91 | for(k = 0; k < sizePhiBins; ++k)
|
---|
| 92 | {
|
---|
[659] | 93 | fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
|
---|
[646] | 94 | }
|
---|
[613] | 95 | }
|
---|
[608] | 96 | }
|
---|
| 97 |
|
---|
[661] | 98 | // for better performance we transform map of sets to parallel vectors:
|
---|
[659] | 99 | // vector< double > and vector< vector< double >* >
|
---|
[622] | 100 | for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
|
---|
[813] | 101 | {
|
---|
[659] | 102 | fEtaBins.push_back(itEtaBin->first);
|
---|
| 103 | phiBins = new vector< double >(itEtaBin->second.size());
|
---|
| 104 | fPhiBins.push_back(phiBins);
|
---|
[666] | 105 | phiBins->clear();
|
---|
[659] | 106 | for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
|
---|
[622] | 107 | {
|
---|
[659] | 108 | phiBins->push_back(*itPhiBin);
|
---|
[622] | 109 | }
|
---|
| 110 | }
|
---|
[935] | 111 |
|
---|
[608] | 112 | // read energy fractions for different particles
|
---|
| 113 | param = GetParam("EnergyFraction");
|
---|
| 114 | size = param.GetSize();
|
---|
[935] | 115 |
|
---|
[608] | 116 | // set default energy fractions values
|
---|
| 117 | fFractionMap.clear();
|
---|
| 118 | fFractionMap[0] = make_pair(0.0, 1.0);
|
---|
[935] | 119 |
|
---|
[608] | 120 | for(i = 0; i < size/2; ++i)
|
---|
| 121 | {
|
---|
[703] | 122 | paramFractions = param[i*2 + 1];
|
---|
[608] | 123 | sizeFractions = paramFractions.GetSize();
|
---|
[935] | 124 |
|
---|
[608] | 125 | ecalFraction = paramFractions[0].GetDouble();
|
---|
| 126 | hcalFraction = paramFractions[1].GetDouble();
|
---|
| 127 |
|
---|
[703] | 128 | fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
|
---|
[608] | 129 | }
|
---|
[629] | 130 | /*
|
---|
[626] | 131 | TFractionMap::iterator itFractionMap;
|
---|
| 132 | for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
|
---|
| 133 | {
|
---|
| 134 | cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
|
---|
| 135 | }
|
---|
[629] | 136 | */
|
---|
[935] | 137 | // read resolution formulas
|
---|
[766] | 138 | fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
|
---|
| 139 | fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
|
---|
| 140 |
|
---|
[608] | 141 | // import array with output from other modules
|
---|
[886] | 142 | fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
|
---|
| 143 | fItParticleInputArray = fParticleInputArray->MakeIterator();
|
---|
[608] | 144 |
|
---|
[886] | 145 | fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
|
---|
| 146 | fItTrackInputArray = fTrackInputArray->MakeIterator();
|
---|
| 147 |
|
---|
[766] | 148 | // create output arrays
|
---|
[629] | 149 | fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
|
---|
| 150 | fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
|
---|
[938] | 151 |
|
---|
| 152 | fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
|
---|
| 153 | fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
|
---|
[608] | 154 | }
|
---|
| 155 |
|
---|
| 156 | //------------------------------------------------------------------------------
|
---|
| 157 |
|
---|
| 158 | void Calorimeter::Finish()
|
---|
| 159 | {
|
---|
[1273] | 160 | vector< vector< Double_t >* >::iterator itPhiBin;
|
---|
[886] | 161 | if(fItParticleInputArray) delete fItParticleInputArray;
|
---|
[938] | 162 | if(fItTrackInputArray) delete fItTrackInputArray;
|
---|
[660] | 163 | for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
|
---|
| 164 | {
|
---|
| 165 | delete *itPhiBin;
|
---|
| 166 | }
|
---|
[608] | 167 | }
|
---|
| 168 |
|
---|
| 169 | //------------------------------------------------------------------------------
|
---|
| 170 |
|
---|
| 171 | void Calorimeter::Process()
|
---|
| 172 | {
|
---|
[886] | 173 | Candidate *particle, *track;
|
---|
[680] | 174 | TLorentzVector position, momentum;
|
---|
[1002] | 175 | Short_t etaBin, phiBin, flags;
|
---|
| 176 | Int_t number;
|
---|
[745] | 177 | Long64_t towerHit, towerEtaPhi, hitEtaPhi;
|
---|
[730] | 178 | Double_t ecalFraction, hcalFraction;
|
---|
[735] | 179 | Double_t ecalEnergy, hcalEnergy;
|
---|
[608] | 180 | Int_t pdgCode;
|
---|
| 181 |
|
---|
| 182 | TFractionMap::iterator itFractionMap;
|
---|
| 183 |
|
---|
[659] | 184 | vector< Double_t >::iterator itEtaBin;
|
---|
| 185 | vector< Double_t >::iterator itPhiBin;
|
---|
| 186 | vector< Double_t > *phiBins;
|
---|
[622] | 187 |
|
---|
[745] | 188 | vector< Long64_t >::iterator itTowerHits;
|
---|
[935] | 189 |
|
---|
[687] | 190 | DelphesFactory *factory = GetFactory();
|
---|
[661] | 191 | fTowerHits.clear();
|
---|
[1273] | 192 | fTowerECalFractions.clear();
|
---|
| 193 | fTowerHCalFractions.clear();
|
---|
| 194 | fTrackECalFractions.clear();
|
---|
| 195 | fTrackHCalFractions.clear();
|
---|
[730] | 196 |
|
---|
[608] | 197 | // loop over all particles
|
---|
[886] | 198 | fItParticleInputArray->Reset();
|
---|
[666] | 199 | number = -1;
|
---|
[886] | 200 | while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
|
---|
[608] | 201 | {
|
---|
[886] | 202 | const TLorentzVector &particlePosition = particle->Position;
|
---|
[666] | 203 | ++number;
|
---|
[608] | 204 |
|
---|
[886] | 205 | pdgCode = TMath::Abs(particle->PID);
|
---|
[626] | 206 |
|
---|
| 207 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 208 | if(itFractionMap == fFractionMap.end())
|
---|
| 209 | {
|
---|
| 210 | itFractionMap = fFractionMap.find(0);
|
---|
| 211 | }
|
---|
| 212 |
|
---|
| 213 | ecalFraction = itFractionMap->second.first;
|
---|
| 214 | hcalFraction = itFractionMap->second.second;
|
---|
[666] | 215 |
|
---|
[1273] | 216 | fTowerECalFractions.push_back(ecalFraction);
|
---|
| 217 | fTowerHCalFractions.push_back(hcalFraction);
|
---|
[935] | 218 |
|
---|
[626] | 219 | if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
|
---|
[935] | 220 |
|
---|
[730] | 221 | // find eta bin [1, fEtaBins.size - 1]
|
---|
[886] | 222 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
|
---|
[659] | 223 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
[730] | 224 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
[935] | 225 |
|
---|
[730] | 226 | // phi bins for given eta bin
|
---|
[659] | 227 | phiBins = fPhiBins[etaBin];
|
---|
[730] | 228 |
|
---|
| 229 | // find phi bin [1, phiBins.size - 1]
|
---|
[886] | 230 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
|
---|
[659] | 231 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
[666] | 232 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
[608] | 233 |
|
---|
[1233] | 234 | flags = 0;
|
---|
[1273] | 235 | flags |= (pdgCode == 11 || pdgCode == 22) << 1;
|
---|
[745] | 236 |
|
---|
[1004] | 237 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
|
---|
| 238 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
[935] | 239 |
|
---|
[661] | 240 | fTowerHits.push_back(towerHit);
|
---|
| 241 | }
|
---|
| 242 |
|
---|
[886] | 243 | // loop over all tracks
|
---|
| 244 | fItTrackInputArray->Reset();
|
---|
| 245 | number = -1;
|
---|
| 246 | while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
|
---|
| 247 | {
|
---|
| 248 | const TLorentzVector &trackPosition = track->Position;
|
---|
| 249 | ++number;
|
---|
| 250 |
|
---|
[1233] | 251 | pdgCode = TMath::Abs(track->PID);
|
---|
| 252 |
|
---|
| 253 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 254 | if(itFractionMap == fFractionMap.end())
|
---|
| 255 | {
|
---|
| 256 | itFractionMap = fFractionMap.find(0);
|
---|
| 257 | }
|
---|
| 258 |
|
---|
| 259 | ecalFraction = itFractionMap->second.first;
|
---|
| 260 | hcalFraction = itFractionMap->second.second;
|
---|
| 261 |
|
---|
[1273] | 262 | fTrackECalFractions.push_back(ecalFraction);
|
---|
| 263 | fTrackHCalFractions.push_back(hcalFraction);
|
---|
| 264 |
|
---|
[886] | 265 | // find eta bin [1, fEtaBins.size - 1]
|
---|
| 266 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
|
---|
| 267 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
| 268 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
[935] | 269 |
|
---|
[886] | 270 | // phi bins for given eta bin
|
---|
| 271 | phiBins = fPhiBins[etaBin];
|
---|
| 272 |
|
---|
| 273 | // find phi bin [1, phiBins.size - 1]
|
---|
| 274 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
|
---|
| 275 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
| 276 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
| 277 |
|
---|
[1233] | 278 | flags = 1;
|
---|
| 279 |
|
---|
[1004] | 280 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
|
---|
[1233] | 281 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
[935] | 282 |
|
---|
[886] | 283 | fTowerHits.push_back(towerHit);
|
---|
| 284 | }
|
---|
| 285 |
|
---|
[894] | 286 | // all hits are sorted first by eta bin number, then by phi bin number,
|
---|
| 287 | // then by flags and then by particle or track number
|
---|
[661] | 288 | sort(fTowerHits.begin(), fTowerHits.end());
|
---|
| 289 |
|
---|
[730] | 290 | // loop over all hits
|
---|
[661] | 291 | towerEtaPhi = 0;
|
---|
[732] | 292 | fTower = 0;
|
---|
[661] | 293 | for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
|
---|
| 294 | {
|
---|
| 295 | towerHit = (*itTowerHits);
|
---|
[1004] | 296 | flags = (towerHit >> 24) & 0x00000000000000FFLL;
|
---|
| 297 | number = (towerHit) & 0x0000000000FFFFFFLL;
|
---|
| 298 | hitEtaPhi = towerHit >> 32;
|
---|
[730] | 299 |
|
---|
[661] | 300 | if(towerEtaPhi != hitEtaPhi)
|
---|
[608] | 301 | {
|
---|
[730] | 302 | // switch to next tower
|
---|
[661] | 303 | towerEtaPhi = hitEtaPhi;
|
---|
[730] | 304 |
|
---|
[661] | 305 | // finalize previous tower
|
---|
[730] | 306 | FinalizeTower();
|
---|
[661] | 307 |
|
---|
| 308 | // create new tower
|
---|
[730] | 309 | fTower = factory->NewCandidate();
|
---|
[661] | 310 |
|
---|
[1004] | 311 | phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
|
---|
| 312 | etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
|
---|
[730] | 313 |
|
---|
| 314 | // phi bins for given eta bin
|
---|
[661] | 315 | phiBins = fPhiBins[etaBin];
|
---|
[608] | 316 |
|
---|
[730] | 317 | // calculate eta and phi of the tower's center
|
---|
[731] | 318 | fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
|
---|
| 319 | fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
|
---|
[626] | 320 |
|
---|
[898] | 321 | fTowerEdges[0] = fEtaBins[etaBin - 1];
|
---|
| 322 | fTowerEdges[1] = fEtaBins[etaBin];
|
---|
| 323 | fTowerEdges[2] = (*phiBins)[phiBin - 1];
|
---|
| 324 | fTowerEdges[3] = (*phiBins)[phiBin];
|
---|
| 325 |
|
---|
[735] | 326 | fTowerECalEnergy = 0.0;
|
---|
| 327 | fTowerHCalEnergy = 0.0;
|
---|
| 328 |
|
---|
[1273] | 329 | fTrackECalEnergy = 0.0;
|
---|
| 330 | fTrackHCalEnergy = 0.0;
|
---|
| 331 |
|
---|
| 332 | fTowerTrackHits = 0;
|
---|
[1233] | 333 | fTowerPhotonHits = 0;
|
---|
[935] | 334 |
|
---|
[1235] | 335 | fTowerTrackArray->Clear();
|
---|
[661] | 336 | }
|
---|
| 337 |
|
---|
[886] | 338 | // check for track hits
|
---|
[1233] | 339 | if(flags & 1)
|
---|
[886] | 340 | {
|
---|
[1273] | 341 | ++fTowerTrackHits;
|
---|
| 342 |
|
---|
[886] | 343 | track = static_cast<Candidate*>(fTrackInputArray->At(number));
|
---|
[1273] | 344 | momentum = track->Momentum;
|
---|
[1233] | 345 |
|
---|
[1273] | 346 | ecalEnergy = momentum.E() * fTrackECalFractions[number];
|
---|
| 347 | hcalEnergy = momentum.E() * fTrackHCalFractions[number];
|
---|
| 348 |
|
---|
| 349 | fTrackECalEnergy += ecalEnergy;
|
---|
| 350 | fTrackHCalEnergy += hcalEnergy;
|
---|
| 351 |
|
---|
[886] | 352 | fTowerTrackArray->Add(track);
|
---|
[1235] | 353 |
|
---|
[886] | 354 | continue;
|
---|
| 355 | }
|
---|
| 356 |
|
---|
[1233] | 357 | // check for photon and electron hits in current tower
|
---|
[1273] | 358 | if(flags & 2) ++fTowerPhotonHits;
|
---|
[1233] | 359 |
|
---|
[886] | 360 | particle = static_cast<Candidate*>(fParticleInputArray->At(number));
|
---|
[888] | 361 | momentum = particle->Momentum;
|
---|
[886] | 362 |
|
---|
[735] | 363 | // fill current tower
|
---|
[1273] | 364 | ecalEnergy = momentum.E() * fTowerECalFractions[number];
|
---|
| 365 | hcalEnergy = momentum.E() * fTowerHCalFractions[number];
|
---|
[931] | 366 |
|
---|
[735] | 367 | fTowerECalEnergy += ecalEnergy;
|
---|
| 368 | fTowerHCalEnergy += hcalEnergy;
|
---|
[1242] | 369 |
|
---|
| 370 | fTower->AddCandidate(particle);
|
---|
[608] | 371 | }
|
---|
| 372 |
|
---|
[661] | 373 | // finalize last tower
|
---|
[730] | 374 | FinalizeTower();
|
---|
| 375 | }
|
---|
| 376 |
|
---|
| 377 | //------------------------------------------------------------------------------
|
---|
[931] | 378 |
|
---|
[730] | 379 | void Calorimeter::FinalizeTower()
|
---|
| 380 | {
|
---|
[1273] | 381 | Candidate *track, *tower;
|
---|
[1078] | 382 | Double_t energy, pt, eta, phi;
|
---|
[735] | 383 | Double_t ecalEnergy, hcalEnergy;
|
---|
[1145] | 384 |
|
---|
[734] | 385 | if(!fTower) return;
|
---|
[733] | 386 |
|
---|
[1142] | 387 | // ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
[1145] | 388 | // if(ecalEnergy < 0.0) ecalEnergy = 0.0;
|
---|
[733] | 389 |
|
---|
[1145] | 390 | ecalEnergy = LogNormal(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
[733] | 391 |
|
---|
[1145] | 392 | // hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
| 393 | // if(hcalEnergy < 0.0) hcalEnergy = 0.0;
|
---|
| 394 |
|
---|
| 395 | hcalEnergy = LogNormal(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
| 396 |
|
---|
[735] | 397 | energy = ecalEnergy + hcalEnergy;
|
---|
[629] | 398 |
|
---|
[1145] | 399 | // eta = fTowerEta;
|
---|
| 400 | // phi = fTowerPhi;
|
---|
[1078] | 401 |
|
---|
| 402 | eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
|
---|
| 403 | phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
|
---|
| 404 |
|
---|
[1086] | 405 | pt = energy / TMath::CosH(eta);
|
---|
| 406 |
|
---|
[1078] | 407 | fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
| 408 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
[735] | 409 | fTower->Eem = ecalEnergy;
|
---|
| 410 | fTower->Ehad = hcalEnergy;
|
---|
[734] | 411 |
|
---|
[898] | 412 | fTower->Edges[0] = fTowerEdges[0];
|
---|
| 413 | fTower->Edges[1] = fTowerEdges[1];
|
---|
| 414 | fTower->Edges[2] = fTowerEdges[2];
|
---|
| 415 | fTower->Edges[3] = fTowerEdges[3];
|
---|
[886] | 416 |
|
---|
[931] | 417 | // fill calorimeter towers and photon candidates
|
---|
| 418 | if(energy > 0.0)
|
---|
| 419 | {
|
---|
[1273] | 420 | if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
|
---|
[935] | 421 | {
|
---|
[931] | 422 | fPhotonOutputArray->Add(fTower);
|
---|
| 423 | }
|
---|
| 424 |
|
---|
| 425 | fTowerOutputArray->Add(fTower);
|
---|
| 426 | }
|
---|
[935] | 427 |
|
---|
[931] | 428 | // fill energy flow candidates
|
---|
[1273] | 429 |
|
---|
| 430 | // save all the tracks as energy flow tracks
|
---|
| 431 | fItTowerTrackArray->Reset();
|
---|
| 432 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
|
---|
[735] | 433 | {
|
---|
[1273] | 434 | fEFlowTrackOutputArray->Add(track);
|
---|
[886] | 435 | }
|
---|
[1241] | 436 |
|
---|
[1273] | 437 | ecalEnergy -= fTrackECalEnergy;
|
---|
| 438 | if(ecalEnergy < 0.0) ecalEnergy = 0.0;
|
---|
[1241] | 439 |
|
---|
[1273] | 440 | hcalEnergy -= fTrackHCalEnergy;
|
---|
| 441 | if(hcalEnergy < 0.0) hcalEnergy = 0.0;
|
---|
[1241] | 442 |
|
---|
[1273] | 443 | energy = ecalEnergy + hcalEnergy;
|
---|
[1241] | 444 |
|
---|
[1273] | 445 | // save ECAL and/or HCAL energy excess as an energy flow tower
|
---|
| 446 | if(energy > 0.0)
|
---|
| 447 | {
|
---|
| 448 | // create new tower
|
---|
| 449 | tower = static_cast<Candidate*>(fTower->Clone());
|
---|
[1241] | 450 |
|
---|
[1273] | 451 | pt = energy / TMath::CosH(eta);
|
---|
[1241] | 452 |
|
---|
[1273] | 453 | tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
| 454 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
| 455 | tower->Eem = ecalEnergy;
|
---|
| 456 | tower->Ehad = hcalEnergy;
|
---|
[1241] | 457 |
|
---|
[1273] | 458 | tower->Edges[0] = fTowerEdges[0];
|
---|
| 459 | tower->Edges[1] = fTowerEdges[1];
|
---|
| 460 | tower->Edges[2] = fTowerEdges[2];
|
---|
| 461 | tower->Edges[3] = fTowerEdges[3];
|
---|
[1241] | 462 |
|
---|
[1273] | 463 | fEFlowTowerOutputArray->Add(tower);
|
---|
[931] | 464 | }
|
---|
[608] | 465 | }
|
---|
| 466 |
|
---|
| 467 | //------------------------------------------------------------------------------
|
---|
[1142] | 468 |
|
---|
[1145] | 469 | Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
|
---|
| 470 | {
|
---|
| 471 | Double_t a, b;
|
---|
| 472 |
|
---|
| 473 | if(mean > 0.0)
|
---|
[1142] | 474 | {
|
---|
[1145] | 475 | b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
|
---|
| 476 | a = TMath::Log(mean) - 0.5*b*b;
|
---|
| 477 |
|
---|
| 478 | return TMath::Exp(a + b*gRandom->Gaus(0, 1));
|
---|
[1142] | 479 | }
|
---|
[1145] | 480 | else
|
---|
| 481 | {
|
---|
| 482 | return 0.0;
|
---|
| 483 | }
|
---|
| 484 | }
|
---|