[b906e2b] | 1 |
|
---|
| 2 | /** \class OldCalorimeter
|
---|
| 3 | *
|
---|
| 4 | * Fills calorimeter towers, performs calorimeter resolution smearing,
|
---|
| 5 | * preselects towers hit by photons and creates energy flow objects.
|
---|
| 6 | *
|
---|
| 7 | * $Date$
|
---|
| 8 | * $Revision$
|
---|
| 9 | *
|
---|
| 10 | *
|
---|
| 11 | * \author P. Demin - UCL, Louvain-la-Neuve
|
---|
| 12 | *
|
---|
| 13 | */
|
---|
| 14 |
|
---|
| 15 | #include "modules/OldCalorimeter.h"
|
---|
| 16 |
|
---|
| 17 | #include "classes/DelphesClasses.h"
|
---|
| 18 | #include "classes/DelphesFactory.h"
|
---|
| 19 | #include "classes/DelphesFormula.h"
|
---|
| 20 |
|
---|
| 21 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
| 22 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
| 23 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
| 24 |
|
---|
| 25 | #include "TMath.h"
|
---|
| 26 | #include "TString.h"
|
---|
| 27 | #include "TFormula.h"
|
---|
| 28 | #include "TRandom3.h"
|
---|
| 29 | #include "TObjArray.h"
|
---|
| 30 | #include "TDatabasePDG.h"
|
---|
| 31 | #include "TLorentzVector.h"
|
---|
| 32 |
|
---|
| 33 | #include <algorithm>
|
---|
| 34 | #include <stdexcept>
|
---|
| 35 | #include <iostream>
|
---|
| 36 | #include <sstream>
|
---|
| 37 |
|
---|
| 38 | using namespace std;
|
---|
| 39 |
|
---|
| 40 | //------------------------------------------------------------------------------
|
---|
| 41 |
|
---|
| 42 | OldCalorimeter::OldCalorimeter() :
|
---|
| 43 | fECalResolutionFormula(0), fHCalResolutionFormula(0),
|
---|
| 44 | fItParticleInputArray(0), fItTrackInputArray(0),
|
---|
| 45 | fTowerECalArray(0), fItTowerECalArray(0),
|
---|
| 46 | fTowerHCalArray(0), fItTowerHCalArray(0),
|
---|
| 47 | fTowerTrackArray(0), fItTowerTrackArray(0),
|
---|
| 48 | fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
|
---|
| 49 | fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
|
---|
| 50 | {
|
---|
| 51 | fECalResolutionFormula = new DelphesFormula;
|
---|
| 52 | fHCalResolutionFormula = new DelphesFormula;
|
---|
| 53 |
|
---|
| 54 | fTowerECalArray = new TObjArray;
|
---|
| 55 | fItTowerECalArray = fTowerECalArray->MakeIterator();
|
---|
| 56 | fTowerHCalArray = new TObjArray;
|
---|
| 57 | fItTowerHCalArray = fTowerHCalArray->MakeIterator();
|
---|
| 58 |
|
---|
| 59 | fTowerTrackArray = new TObjArray;
|
---|
| 60 | fItTowerTrackArray = fTowerTrackArray->MakeIterator();
|
---|
| 61 | fTowerECalTrackArray = new TObjArray;
|
---|
| 62 | fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
|
---|
| 63 | fTowerHCalTrackArray = new TObjArray;
|
---|
| 64 | fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
|
---|
| 65 | }
|
---|
| 66 |
|
---|
| 67 | //------------------------------------------------------------------------------
|
---|
| 68 |
|
---|
| 69 | OldCalorimeter::~OldCalorimeter()
|
---|
| 70 | {
|
---|
| 71 | if(fECalResolutionFormula) delete fECalResolutionFormula;
|
---|
| 72 | if(fHCalResolutionFormula) delete fHCalResolutionFormula;
|
---|
| 73 |
|
---|
| 74 | if(fTowerECalArray) delete fTowerECalArray;
|
---|
| 75 | if(fItTowerECalArray) delete fItTowerECalArray;
|
---|
| 76 | if(fTowerHCalArray) delete fTowerHCalArray;
|
---|
| 77 | if(fItTowerHCalArray) delete fItTowerHCalArray;
|
---|
| 78 |
|
---|
| 79 | if(fTowerTrackArray) delete fTowerTrackArray;
|
---|
| 80 | if(fItTowerTrackArray) delete fItTowerTrackArray;
|
---|
| 81 | if(fTowerECalTrackArray) delete fTowerECalTrackArray;
|
---|
| 82 | if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
|
---|
| 83 | if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
|
---|
| 84 | if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
|
---|
| 85 |
|
---|
| 86 | //------------------------------------------------------------------------------
|
---|
| 87 |
|
---|
| 88 | void OldCalorimeter::Init()
|
---|
| 89 | {
|
---|
| 90 | ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
|
---|
| 91 | Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
|
---|
| 92 | Double_t ecalFraction, hcalFraction;
|
---|
| 93 | TBinMap::iterator itEtaBin;
|
---|
| 94 | set< Double_t >::iterator itPhiBin;
|
---|
| 95 | vector< Double_t > *phiBins;
|
---|
| 96 |
|
---|
| 97 | // read eta and phi bins
|
---|
| 98 | param = GetParam("EtaPhiBins");
|
---|
| 99 | size = param.GetSize();
|
---|
| 100 | fBinMap.clear();
|
---|
| 101 | fEtaBins.clear();
|
---|
| 102 | fPhiBins.clear();
|
---|
| 103 | for(i = 0; i < size/2; ++i)
|
---|
| 104 | {
|
---|
| 105 | paramEtaBins = param[i*2];
|
---|
| 106 | sizeEtaBins = paramEtaBins.GetSize();
|
---|
| 107 | paramPhiBins = param[i*2 + 1];
|
---|
| 108 | sizePhiBins = paramPhiBins.GetSize();
|
---|
| 109 |
|
---|
| 110 | for(j = 0; j < sizeEtaBins; ++j)
|
---|
| 111 | {
|
---|
| 112 | for(k = 0; k < sizePhiBins; ++k)
|
---|
| 113 | {
|
---|
| 114 | fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
|
---|
| 115 | }
|
---|
| 116 | }
|
---|
| 117 | }
|
---|
| 118 |
|
---|
| 119 | // for better performance we transform map of sets to parallel vectors:
|
---|
| 120 | // vector< double > and vector< vector< double >* >
|
---|
| 121 | for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
|
---|
| 122 | {
|
---|
| 123 | fEtaBins.push_back(itEtaBin->first);
|
---|
| 124 | phiBins = new vector< double >(itEtaBin->second.size());
|
---|
| 125 | fPhiBins.push_back(phiBins);
|
---|
| 126 | phiBins->clear();
|
---|
| 127 | for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
|
---|
| 128 | {
|
---|
| 129 | phiBins->push_back(*itPhiBin);
|
---|
| 130 | }
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | // read energy fractions for different particles
|
---|
| 134 | param = GetParam("EnergyFraction");
|
---|
| 135 | size = param.GetSize();
|
---|
| 136 |
|
---|
| 137 | // set default energy fractions values
|
---|
| 138 | fFractionMap.clear();
|
---|
| 139 | fFractionMap[0] = make_pair(0.0, 1.0);
|
---|
| 140 |
|
---|
| 141 | for(i = 0; i < size/2; ++i)
|
---|
| 142 | {
|
---|
| 143 | paramFractions = param[i*2 + 1];
|
---|
| 144 | sizeFractions = paramFractions.GetSize();
|
---|
| 145 |
|
---|
| 146 | ecalFraction = paramFractions[0].GetDouble();
|
---|
| 147 | hcalFraction = paramFractions[1].GetDouble();
|
---|
| 148 |
|
---|
| 149 | fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
|
---|
| 150 | }
|
---|
| 151 | /*
|
---|
| 152 | TFractionMap::iterator itFractionMap;
|
---|
| 153 | for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
|
---|
| 154 | {
|
---|
| 155 | cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
|
---|
| 156 | }
|
---|
| 157 | */
|
---|
| 158 | // read resolution formulas
|
---|
| 159 | fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
|
---|
| 160 | fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
|
---|
| 161 |
|
---|
| 162 | // import array with output from other modules
|
---|
| 163 | fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
|
---|
| 164 | fItParticleInputArray = fParticleInputArray->MakeIterator();
|
---|
| 165 |
|
---|
| 166 | fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
|
---|
| 167 | fItTrackInputArray = fTrackInputArray->MakeIterator();
|
---|
| 168 |
|
---|
| 169 | // create output arrays
|
---|
| 170 | fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
|
---|
| 171 | fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
|
---|
| 172 |
|
---|
| 173 | fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
|
---|
| 174 | fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
|
---|
| 175 | }
|
---|
| 176 |
|
---|
| 177 | //------------------------------------------------------------------------------
|
---|
| 178 |
|
---|
| 179 | void OldCalorimeter::Finish()
|
---|
| 180 | {
|
---|
| 181 | vector< vector< Double_t>* >::iterator itPhiBin;
|
---|
| 182 | if(fItParticleInputArray) delete fItParticleInputArray;
|
---|
| 183 | if(fItTrackInputArray) delete fItTrackInputArray;
|
---|
| 184 | for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
|
---|
| 185 | {
|
---|
| 186 | delete *itPhiBin;
|
---|
| 187 | }
|
---|
| 188 | }
|
---|
| 189 |
|
---|
| 190 | //------------------------------------------------------------------------------
|
---|
| 191 |
|
---|
| 192 | void OldCalorimeter::Process()
|
---|
| 193 | {
|
---|
| 194 | Candidate *particle, *track;
|
---|
| 195 | TLorentzVector position, momentum;
|
---|
| 196 | Short_t etaBin, phiBin, flags;
|
---|
| 197 | Int_t number;
|
---|
| 198 | Long64_t towerHit, towerEtaPhi, hitEtaPhi;
|
---|
| 199 | Double_t ecalFraction, hcalFraction;
|
---|
| 200 | Double_t ecalEnergy, hcalEnergy;
|
---|
| 201 | Int_t pdgCode;
|
---|
| 202 |
|
---|
| 203 | TFractionMap::iterator itFractionMap;
|
---|
| 204 |
|
---|
| 205 | vector< Double_t >::iterator itEtaBin;
|
---|
| 206 | vector< Double_t >::iterator itPhiBin;
|
---|
| 207 | vector< Double_t > *phiBins;
|
---|
| 208 |
|
---|
| 209 | vector< Long64_t >::iterator itTowerHits;
|
---|
| 210 |
|
---|
| 211 | DelphesFactory *factory = GetFactory();
|
---|
| 212 | fTowerHits.clear();
|
---|
| 213 | fECalFractions.clear();
|
---|
| 214 | fHCalFractions.clear();
|
---|
| 215 |
|
---|
| 216 | // loop over all particles
|
---|
| 217 | fItParticleInputArray->Reset();
|
---|
| 218 | number = -1;
|
---|
| 219 | while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
|
---|
| 220 | {
|
---|
| 221 | const TLorentzVector &particlePosition = particle->Position;
|
---|
| 222 | ++number;
|
---|
| 223 |
|
---|
| 224 | pdgCode = TMath::Abs(particle->PID);
|
---|
| 225 |
|
---|
| 226 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 227 | if(itFractionMap == fFractionMap.end())
|
---|
| 228 | {
|
---|
| 229 | itFractionMap = fFractionMap.find(0);
|
---|
| 230 | }
|
---|
| 231 |
|
---|
| 232 | ecalFraction = itFractionMap->second.first;
|
---|
| 233 | hcalFraction = itFractionMap->second.second;
|
---|
| 234 |
|
---|
| 235 | fECalFractions.push_back(ecalFraction);
|
---|
| 236 | fHCalFractions.push_back(hcalFraction);
|
---|
| 237 |
|
---|
| 238 | if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
|
---|
| 239 |
|
---|
| 240 | // find eta bin [1, fEtaBins.size - 1]
|
---|
| 241 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
|
---|
| 242 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
| 243 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
| 244 |
|
---|
| 245 | // phi bins for given eta bin
|
---|
| 246 | phiBins = fPhiBins[etaBin];
|
---|
| 247 |
|
---|
| 248 | // find phi bin [1, phiBins.size - 1]
|
---|
| 249 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
|
---|
| 250 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
| 251 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
| 252 |
|
---|
| 253 | flags = 0;
|
---|
| 254 | flags |= (ecalFraction >= 1.0E-9) << 1;
|
---|
| 255 | flags |= (hcalFraction >= 1.0E-9) << 2;
|
---|
| 256 | flags |= (pdgCode == 11 || pdgCode == 22) << 3;
|
---|
| 257 |
|
---|
| 258 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
|
---|
| 259 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
| 260 |
|
---|
| 261 | fTowerHits.push_back(towerHit);
|
---|
| 262 | }
|
---|
| 263 |
|
---|
| 264 | // loop over all tracks
|
---|
| 265 | fItTrackInputArray->Reset();
|
---|
| 266 | number = -1;
|
---|
| 267 | while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
|
---|
| 268 | {
|
---|
| 269 | const TLorentzVector &trackPosition = track->Position;
|
---|
| 270 | ++number;
|
---|
| 271 |
|
---|
| 272 | pdgCode = TMath::Abs(track->PID);
|
---|
| 273 |
|
---|
| 274 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 275 | if(itFractionMap == fFractionMap.end())
|
---|
| 276 | {
|
---|
| 277 | itFractionMap = fFractionMap.find(0);
|
---|
| 278 | }
|
---|
| 279 |
|
---|
| 280 | ecalFraction = itFractionMap->second.first;
|
---|
| 281 | hcalFraction = itFractionMap->second.second;
|
---|
| 282 |
|
---|
| 283 | // find eta bin [1, fEtaBins.size - 1]
|
---|
| 284 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
|
---|
| 285 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
| 286 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
| 287 |
|
---|
| 288 | // phi bins for given eta bin
|
---|
| 289 | phiBins = fPhiBins[etaBin];
|
---|
| 290 |
|
---|
| 291 | // find phi bin [1, phiBins.size - 1]
|
---|
| 292 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
|
---|
| 293 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
| 294 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
| 295 |
|
---|
| 296 | flags = 1;
|
---|
| 297 | flags |= (ecalFraction >= 1.0E-9) << 1;
|
---|
| 298 | flags |= (hcalFraction >= 1.0E-9) << 2;
|
---|
| 299 |
|
---|
| 300 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
|
---|
| 301 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
| 302 |
|
---|
| 303 | fTowerHits.push_back(towerHit);
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | // all hits are sorted first by eta bin number, then by phi bin number,
|
---|
| 307 | // then by flags and then by particle or track number
|
---|
| 308 | sort(fTowerHits.begin(), fTowerHits.end());
|
---|
| 309 |
|
---|
| 310 | // loop over all hits
|
---|
| 311 | towerEtaPhi = 0;
|
---|
| 312 | fTower = 0;
|
---|
| 313 | for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
|
---|
| 314 | {
|
---|
| 315 | towerHit = (*itTowerHits);
|
---|
| 316 | flags = (towerHit >> 24) & 0x00000000000000FFLL;
|
---|
| 317 | number = (towerHit) & 0x0000000000FFFFFFLL;
|
---|
| 318 | hitEtaPhi = towerHit >> 32;
|
---|
| 319 |
|
---|
| 320 | if(towerEtaPhi != hitEtaPhi)
|
---|
| 321 | {
|
---|
| 322 | // switch to next tower
|
---|
| 323 | towerEtaPhi = hitEtaPhi;
|
---|
| 324 |
|
---|
| 325 | // finalize previous tower
|
---|
| 326 | FinalizeTower();
|
---|
| 327 |
|
---|
| 328 | // create new tower
|
---|
| 329 | fTower = factory->NewCandidate();
|
---|
| 330 |
|
---|
| 331 | phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
|
---|
| 332 | etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
|
---|
| 333 |
|
---|
| 334 | // phi bins for given eta bin
|
---|
| 335 | phiBins = fPhiBins[etaBin];
|
---|
| 336 |
|
---|
| 337 | // calculate eta and phi of the tower's center
|
---|
| 338 | fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
|
---|
| 339 | fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
|
---|
| 340 |
|
---|
| 341 | fTowerEdges[0] = fEtaBins[etaBin - 1];
|
---|
| 342 | fTowerEdges[1] = fEtaBins[etaBin];
|
---|
| 343 | fTowerEdges[2] = (*phiBins)[phiBin - 1];
|
---|
| 344 | fTowerEdges[3] = (*phiBins)[phiBin];
|
---|
| 345 |
|
---|
| 346 | fTowerECalEnergy = 0.0;
|
---|
| 347 | fTowerHCalEnergy = 0.0;
|
---|
| 348 |
|
---|
| 349 | fTowerPhotonHits = 0;
|
---|
| 350 |
|
---|
| 351 | fTowerAllHits = 0;
|
---|
| 352 | fTowerECalHits = 0;
|
---|
| 353 | fTowerHCalHits = 0;
|
---|
| 354 |
|
---|
| 355 | fTowerTrackAllHits = 0;
|
---|
| 356 | fTowerECalTrackHits = 0;
|
---|
| 357 | fTowerHCalTrackHits = 0;
|
---|
| 358 |
|
---|
| 359 | fTowerECalArray->Clear();
|
---|
| 360 | fTowerHCalArray->Clear();
|
---|
| 361 |
|
---|
| 362 | fTowerTrackArray->Clear();
|
---|
| 363 | fTowerECalTrackArray->Clear();
|
---|
| 364 | fTowerHCalTrackArray->Clear();
|
---|
| 365 | }
|
---|
| 366 |
|
---|
| 367 | // check for track hits
|
---|
| 368 | if(flags & 1)
|
---|
| 369 | {
|
---|
| 370 | track = static_cast<Candidate*>(fTrackInputArray->At(number));
|
---|
| 371 |
|
---|
| 372 | ++fTowerTrackAllHits;
|
---|
| 373 | fTowerTrackArray->Add(track);
|
---|
| 374 |
|
---|
| 375 | // check for track ECAL hits
|
---|
| 376 | if(flags & 2)
|
---|
| 377 | {
|
---|
| 378 | ++fTowerECalTrackHits;
|
---|
| 379 | fTowerECalTrackArray->Add(track);
|
---|
| 380 | }
|
---|
| 381 |
|
---|
| 382 | // check for track HCAL hits
|
---|
| 383 | if(flags & 4)
|
---|
| 384 | {
|
---|
| 385 | ++fTowerHCalTrackHits;
|
---|
| 386 | fTowerHCalTrackArray->Add(track);
|
---|
| 387 | }
|
---|
| 388 | continue;
|
---|
| 389 | }
|
---|
| 390 |
|
---|
| 391 | ++fTowerAllHits;
|
---|
| 392 |
|
---|
| 393 | // check for ECAL hits
|
---|
| 394 | if(flags & 2)
|
---|
| 395 | {
|
---|
| 396 | ++fTowerECalHits;
|
---|
| 397 | fTowerECalArray->Add(particle);
|
---|
| 398 | }
|
---|
| 399 |
|
---|
| 400 | // check for HCAL hits
|
---|
| 401 | if(flags & 4)
|
---|
| 402 | {
|
---|
| 403 | ++fTowerHCalHits;
|
---|
| 404 | fTowerHCalArray->Add(particle);
|
---|
| 405 | }
|
---|
| 406 |
|
---|
| 407 | // check for photon and electron hits in current tower
|
---|
| 408 | if(flags & 8) ++fTowerPhotonHits;
|
---|
| 409 |
|
---|
| 410 | particle = static_cast<Candidate*>(fParticleInputArray->At(number));
|
---|
| 411 | momentum = particle->Momentum;
|
---|
| 412 |
|
---|
| 413 | // fill current tower
|
---|
| 414 | ecalEnergy = momentum.E() * fECalFractions[number];
|
---|
| 415 | hcalEnergy = momentum.E() * fHCalFractions[number];
|
---|
| 416 |
|
---|
| 417 | fTowerECalEnergy += ecalEnergy;
|
---|
| 418 | fTowerHCalEnergy += hcalEnergy;
|
---|
| 419 |
|
---|
| 420 | fTower->AddCandidate(particle);
|
---|
| 421 | }
|
---|
| 422 |
|
---|
| 423 | // finalize last tower
|
---|
| 424 | FinalizeTower();
|
---|
| 425 | }
|
---|
| 426 |
|
---|
| 427 | //------------------------------------------------------------------------------
|
---|
| 428 |
|
---|
| 429 | void OldCalorimeter::FinalizeTower()
|
---|
| 430 | {
|
---|
| 431 | Candidate *particle, *track, *tower;
|
---|
| 432 | Double_t energy, pt, eta, phi;
|
---|
| 433 | Double_t ecalEnergy, hcalEnergy;
|
---|
| 434 | TIterator *itTowerTrackArray;
|
---|
| 435 |
|
---|
| 436 | if(!fTower) return;
|
---|
| 437 |
|
---|
| 438 | // ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
| 439 | // if(ecalEnergy < 0.0) ecalEnergy = 0.0;
|
---|
| 440 |
|
---|
| 441 | ecalEnergy = LogNormal(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
|
---|
| 442 |
|
---|
| 443 | // hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
| 444 | // if(hcalEnergy < 0.0) hcalEnergy = 0.0;
|
---|
| 445 |
|
---|
| 446 | hcalEnergy = LogNormal(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
|
---|
| 447 |
|
---|
| 448 | energy = ecalEnergy + hcalEnergy;
|
---|
| 449 |
|
---|
| 450 | // eta = fTowerEta;
|
---|
| 451 | // phi = fTowerPhi;
|
---|
| 452 |
|
---|
| 453 | eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
|
---|
| 454 | phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
|
---|
| 455 |
|
---|
| 456 | pt = energy / TMath::CosH(eta);
|
---|
| 457 |
|
---|
| 458 | fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
| 459 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
| 460 | fTower->Eem = ecalEnergy;
|
---|
| 461 | fTower->Ehad = hcalEnergy;
|
---|
| 462 |
|
---|
| 463 | fTower->Edges[0] = fTowerEdges[0];
|
---|
| 464 | fTower->Edges[1] = fTowerEdges[1];
|
---|
| 465 | fTower->Edges[2] = fTowerEdges[2];
|
---|
| 466 | fTower->Edges[3] = fTowerEdges[3];
|
---|
| 467 |
|
---|
| 468 | // fill calorimeter towers and photon candidates
|
---|
| 469 | if(energy > 0.0)
|
---|
| 470 | {
|
---|
| 471 | if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
|
---|
| 472 | {
|
---|
| 473 | fPhotonOutputArray->Add(fTower);
|
---|
| 474 | }
|
---|
| 475 |
|
---|
| 476 | fTowerOutputArray->Add(fTower);
|
---|
| 477 | }
|
---|
| 478 |
|
---|
| 479 | // fill energy flow candidates
|
---|
| 480 | if(fTowerTrackAllHits == fTowerAllHits)
|
---|
| 481 | {
|
---|
| 482 | fItTowerTrackArray->Reset();
|
---|
| 483 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
|
---|
| 484 | {
|
---|
| 485 | fEFlowTrackOutputArray->Add(track);
|
---|
| 486 | }
|
---|
| 487 | }
|
---|
| 488 | else if(fTowerTrackAllHits > 0 &&
|
---|
| 489 | fTowerECalHits + fTowerHCalHits == fTowerAllHits)
|
---|
| 490 | {
|
---|
| 491 | if(fTowerECalHits == fTowerECalTrackHits &&
|
---|
| 492 | fTowerHCalHits == fTowerHCalTrackHits)
|
---|
| 493 | {
|
---|
| 494 | itTowerTrackArray = fItTowerTrackArray;
|
---|
| 495 | }
|
---|
| 496 | else if(fTowerECalHits == fTowerECalTrackHits)
|
---|
| 497 | {
|
---|
| 498 | itTowerTrackArray = fItTowerECalTrackArray;
|
---|
| 499 |
|
---|
| 500 | if(hcalEnergy > 0.0)
|
---|
| 501 | {
|
---|
| 502 | DelphesFactory *factory = GetFactory();
|
---|
| 503 |
|
---|
| 504 | // create new tower
|
---|
| 505 | tower = factory->NewCandidate();
|
---|
| 506 |
|
---|
| 507 | fItTowerHCalArray->Reset();
|
---|
| 508 | while((particle = static_cast<Candidate*>(fItTowerHCalArray->Next())))
|
---|
| 509 | {
|
---|
| 510 | tower->AddCandidate(particle);
|
---|
| 511 | }
|
---|
| 512 |
|
---|
| 513 | pt = hcalEnergy / TMath::CosH(eta);
|
---|
| 514 |
|
---|
| 515 | tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
| 516 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
|
---|
| 517 | tower->Eem = 0.0;
|
---|
| 518 | tower->Ehad = hcalEnergy;
|
---|
| 519 |
|
---|
| 520 | tower->Edges[0] = fTowerEdges[0];
|
---|
| 521 | tower->Edges[1] = fTowerEdges[1];
|
---|
| 522 | tower->Edges[2] = fTowerEdges[2];
|
---|
| 523 | tower->Edges[3] = fTowerEdges[3];
|
---|
| 524 |
|
---|
| 525 | fEFlowTowerOutputArray->Add(tower);
|
---|
| 526 | }
|
---|
| 527 | }
|
---|
| 528 | else if(fTowerHCalHits == fTowerHCalTrackHits)
|
---|
| 529 | {
|
---|
| 530 | itTowerTrackArray = fItTowerHCalTrackArray;
|
---|
| 531 |
|
---|
| 532 | if(ecalEnergy > 0.0)
|
---|
| 533 | {
|
---|
| 534 | DelphesFactory *factory = GetFactory();
|
---|
| 535 |
|
---|
| 536 | // create new tower
|
---|
| 537 | tower = factory->NewCandidate();
|
---|
| 538 |
|
---|
| 539 | fItTowerECalArray->Reset();
|
---|
| 540 | while((particle = static_cast<Candidate*>(fItTowerECalArray->Next())))
|
---|
| 541 | {
|
---|
| 542 | tower->AddCandidate(particle);
|
---|
| 543 | }
|
---|
| 544 |
|
---|
| 545 | pt = ecalEnergy / TMath::CosH(eta);
|
---|
| 546 |
|
---|
| 547 | tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
|
---|
| 548 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
|
---|
| 549 | tower->Eem = ecalEnergy;
|
---|
| 550 | tower->Ehad = 0.0;
|
---|
| 551 |
|
---|
| 552 | tower->Edges[0] = fTowerEdges[0];
|
---|
| 553 | tower->Edges[1] = fTowerEdges[1];
|
---|
| 554 | tower->Edges[2] = fTowerEdges[2];
|
---|
| 555 | tower->Edges[3] = fTowerEdges[3];
|
---|
| 556 |
|
---|
| 557 | fEFlowTowerOutputArray->Add(tower);
|
---|
| 558 | }
|
---|
| 559 | }
|
---|
| 560 | else
|
---|
| 561 | {
|
---|
| 562 | itTowerTrackArray = 0;
|
---|
| 563 | fEFlowTowerOutputArray->Add(fTower);
|
---|
| 564 | }
|
---|
| 565 |
|
---|
| 566 | if(itTowerTrackArray)
|
---|
| 567 | {
|
---|
| 568 | itTowerTrackArray->Reset();
|
---|
| 569 | while((track = static_cast<Candidate*>(itTowerTrackArray->Next())))
|
---|
| 570 | {
|
---|
| 571 | fEFlowTrackOutputArray->Add(track);
|
---|
| 572 | }
|
---|
| 573 | }
|
---|
| 574 | }
|
---|
| 575 | else if(energy > 0.0)
|
---|
| 576 | {
|
---|
| 577 | fEFlowTowerOutputArray->Add(fTower);
|
---|
| 578 | }
|
---|
| 579 | }
|
---|
| 580 |
|
---|
| 581 | //------------------------------------------------------------------------------
|
---|
| 582 |
|
---|
| 583 | Double_t OldCalorimeter::LogNormal(Double_t mean, Double_t sigma)
|
---|
| 584 | {
|
---|
| 585 | Double_t a, b;
|
---|
| 586 |
|
---|
| 587 | if(mean > 0.0)
|
---|
| 588 | {
|
---|
| 589 | b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
|
---|
| 590 | a = TMath::Log(mean) - 0.5*b*b;
|
---|
| 591 |
|
---|
| 592 | return TMath::Exp(a + b*gRandom->Gaus(0, 1));
|
---|
| 593 | }
|
---|
| 594 | else
|
---|
| 595 | {
|
---|
| 596 | return 0.0;
|
---|
| 597 | }
|
---|
| 598 | }
|
---|