[ef8a06d] | 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 | *
|
---|
[a1b19ea] | 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))
|
---|
[fd4b326] | 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.
|
---|
[a1b19ea] | 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 |
|
---|
[ef8a06d] | 29 | *
|
---|
[a1b19ea] | 30 | * \author M. Selvaggi - CERN
|
---|
[ef8a06d] | 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 |
|
---|
| 57 | using namespace std;
|
---|
| 58 |
|
---|
| 59 | //------------------------------------------------------------------------------
|
---|
| 60 |
|
---|
| 61 | DualReadoutCalorimeter::DualReadoutCalorimeter() :
|
---|
| 62 | fECalResolutionFormula(0), fHCalResolutionFormula(0),
|
---|
| 63 | fItParticleInputArray(0), fItTrackInputArray(0)
|
---|
| 64 | {
|
---|
[fd4b326] | 65 |
|
---|
[ef8a06d] | 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();
|
---|
[a1b19ea] | 74 |
|
---|
| 75 | fTowerTrackArray = new TObjArray;
|
---|
| 76 | fItTowerTrackArray = fTowerTrackArray->MakeIterator();
|
---|
[fd4b326] | 77 |
|
---|
[ef8a06d] | 78 | }
|
---|
| 79 |
|
---|
| 80 | //------------------------------------------------------------------------------
|
---|
| 81 |
|
---|
| 82 | DualReadoutCalorimeter::~DualReadoutCalorimeter()
|
---|
| 83 | {
|
---|
[fd4b326] | 84 |
|
---|
[ef8a06d] | 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;
|
---|
[a1b19ea] | 93 |
|
---|
| 94 | if(fTowerTrackArray) delete fTowerTrackArray;
|
---|
| 95 | if(fItTowerTrackArray) delete fItTowerTrackArray;
|
---|
[fd4b326] | 96 |
|
---|
[ef8a06d] | 97 | }
|
---|
| 98 |
|
---|
| 99 | //------------------------------------------------------------------------------
|
---|
| 100 |
|
---|
| 101 | void 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);
|
---|
[de2e39d] | 174 | fEnergyMin = GetDouble("EnergyMin", 0.0);
|
---|
[ef8a06d] | 175 |
|
---|
| 176 | fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
|
---|
| 177 | fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
|
---|
[de2e39d] | 178 | fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
|
---|
[ef8a06d] | 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 |
|
---|
| 205 | void 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 |
|
---|
| 218 | void 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;
|
---|
[a1b19ea] | 227 | Double_t ecalSigma, hcalSigma, sigma;
|
---|
| 228 | Double_t energyGuess, energy;
|
---|
[ef8a06d] | 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;
|
---|
[5eda6767] | 249 | fTowerRmax=0.;
|
---|
[61dccd3] | 250 |
|
---|
| 251 | //cout<<"--------- new event ---------- "<<endl;
|
---|
| 252 |
|
---|
[ef8a06d] | 253 | while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
|
---|
| 254 | {
|
---|
| 255 | const TLorentzVector &particlePosition = particle->Position;
|
---|
| 256 | ++number;
|
---|
| 257 |
|
---|
[5eda6767] | 258 | // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
|
---|
| 259 | if (particlePosition.Perp() > fTowerRmax)
|
---|
| 260 | fTowerRmax=particlePosition.Perp();
|
---|
| 261 |
|
---|
[ef8a06d] | 262 | pdgCode = TMath::Abs(particle->PID);
|
---|
| 263 |
|
---|
| 264 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 265 | if(itFractionMap == fFractionMap.end())
|
---|
| 266 | {
|
---|
| 267 | itFractionMap = fFractionMap.find(0);
|
---|
| 268 | }
|
---|
| 269 |
|
---|
| 270 | ecalFraction = itFractionMap->second.first;
|
---|
| 271 | hcalFraction = itFractionMap->second.second;
|
---|
| 272 |
|
---|
| 273 | fECalTowerFractions.push_back(ecalFraction);
|
---|
| 274 | fHCalTowerFractions.push_back(hcalFraction);
|
---|
| 275 |
|
---|
| 276 | if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
|
---|
| 277 |
|
---|
| 278 | // find eta bin [1, fEtaBins.size - 1]
|
---|
| 279 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
|
---|
| 280 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
| 281 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
| 282 |
|
---|
| 283 | // phi bins for given eta bin
|
---|
| 284 | phiBins = fPhiBins[etaBin];
|
---|
| 285 |
|
---|
| 286 | // find phi bin [1, phiBins.size - 1]
|
---|
| 287 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
|
---|
| 288 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
| 289 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
| 290 |
|
---|
| 291 | flags = 0;
|
---|
| 292 | flags |= (pdgCode == 11 || pdgCode == 22) << 1;
|
---|
| 293 |
|
---|
| 294 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
|
---|
| 295 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
| 296 |
|
---|
| 297 | fTowerHits.push_back(towerHit);
|
---|
| 298 | }
|
---|
| 299 |
|
---|
| 300 | // loop over all tracks
|
---|
| 301 | fItTrackInputArray->Reset();
|
---|
| 302 | number = -1;
|
---|
| 303 | while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
|
---|
| 304 | {
|
---|
| 305 | const TLorentzVector &trackPosition = track->Position;
|
---|
| 306 | ++number;
|
---|
| 307 |
|
---|
| 308 | pdgCode = TMath::Abs(track->PID);
|
---|
| 309 |
|
---|
| 310 | itFractionMap = fFractionMap.find(pdgCode);
|
---|
| 311 | if(itFractionMap == fFractionMap.end())
|
---|
| 312 | {
|
---|
| 313 | itFractionMap = fFractionMap.find(0);
|
---|
| 314 | }
|
---|
| 315 |
|
---|
| 316 | ecalFraction = itFractionMap->second.first;
|
---|
| 317 | hcalFraction = itFractionMap->second.second;
|
---|
| 318 |
|
---|
| 319 | fECalTrackFractions.push_back(ecalFraction);
|
---|
| 320 | fHCalTrackFractions.push_back(hcalFraction);
|
---|
| 321 |
|
---|
| 322 | // find eta bin [1, fEtaBins.size - 1]
|
---|
| 323 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
|
---|
| 324 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
|
---|
| 325 | etaBin = distance(fEtaBins.begin(), itEtaBin);
|
---|
| 326 |
|
---|
| 327 | // phi bins for given eta bin
|
---|
| 328 | phiBins = fPhiBins[etaBin];
|
---|
| 329 |
|
---|
| 330 | // find phi bin [1, phiBins.size - 1]
|
---|
| 331 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
|
---|
| 332 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
|
---|
| 333 | phiBin = distance(phiBins->begin(), itPhiBin);
|
---|
| 334 |
|
---|
| 335 | flags = 1;
|
---|
| 336 |
|
---|
| 337 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
|
---|
| 338 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
|
---|
| 339 |
|
---|
| 340 | fTowerHits.push_back(towerHit);
|
---|
| 341 | }
|
---|
| 342 |
|
---|
| 343 | // all hits are sorted first by eta bin number, then by phi bin number,
|
---|
| 344 | // then by flags and then by particle or track number
|
---|
| 345 | sort(fTowerHits.begin(), fTowerHits.end());
|
---|
| 346 |
|
---|
| 347 | // loop over all hits
|
---|
| 348 | towerEtaPhi = 0;
|
---|
| 349 | fTower = 0;
|
---|
| 350 | for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
|
---|
| 351 | {
|
---|
| 352 | towerHit = (*itTowerHits);
|
---|
| 353 | flags = (towerHit >> 24) & 0x00000000000000FFLL;
|
---|
| 354 | number = (towerHit) & 0x0000000000FFFFFFLL;
|
---|
| 355 | hitEtaPhi = towerHit >> 32;
|
---|
| 356 |
|
---|
| 357 | if(towerEtaPhi != hitEtaPhi)
|
---|
| 358 | {
|
---|
| 359 | // switch to next tower
|
---|
| 360 | towerEtaPhi = hitEtaPhi;
|
---|
| 361 |
|
---|
| 362 | // finalize previous tower
|
---|
| 363 | FinalizeTower();
|
---|
| 364 |
|
---|
| 365 | // create new tower
|
---|
| 366 | fTower = factory->NewCandidate();
|
---|
| 367 |
|
---|
| 368 | phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
|
---|
| 369 | etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
|
---|
| 370 |
|
---|
| 371 | // phi bins for given eta bin
|
---|
| 372 | phiBins = fPhiBins[etaBin];
|
---|
| 373 |
|
---|
| 374 | // calculate eta and phi of the tower's center
|
---|
| 375 | fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
|
---|
| 376 | fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
|
---|
| 377 |
|
---|
| 378 | fTowerEdges[0] = fEtaBins[etaBin - 1];
|
---|
| 379 | fTowerEdges[1] = fEtaBins[etaBin];
|
---|
| 380 | fTowerEdges[2] = (*phiBins)[phiBin - 1];
|
---|
| 381 | fTowerEdges[3] = (*phiBins)[phiBin];
|
---|
| 382 |
|
---|
| 383 | fECalTowerEnergy = 0.0;
|
---|
| 384 | fHCalTowerEnergy = 0.0;
|
---|
| 385 |
|
---|
| 386 | fECalTrackEnergy = 0.0;
|
---|
| 387 | fHCalTrackEnergy = 0.0;
|
---|
[a1b19ea] | 388 | fTrackEnergy = 0.0;
|
---|
[fd4b326] | 389 |
|
---|
[ef8a06d] | 390 | fECalTrackSigma = 0.0;
|
---|
| 391 | fHCalTrackSigma = 0.0;
|
---|
[a1b19ea] | 392 | fTrackSigma = 0.0;
|
---|
[fd4b326] | 393 |
|
---|
[ef8a06d] | 394 | fTowerTrackHits = 0;
|
---|
| 395 | fTowerPhotonHits = 0;
|
---|
| 396 |
|
---|
[61dccd3] | 397 | fTowerTime = 0.0;
|
---|
| 398 | fTowerTimeWeight = 0.0;
|
---|
| 399 |
|
---|
[ef8a06d] | 400 | fECalTowerTrackArray->Clear();
|
---|
| 401 | fHCalTowerTrackArray->Clear();
|
---|
[a1b19ea] | 402 | fTowerTrackArray->Clear();
|
---|
[fd4b326] | 403 |
|
---|
[ef8a06d] | 404 | }
|
---|
| 405 |
|
---|
| 406 | // check for track hits
|
---|
| 407 | if(flags & 1)
|
---|
| 408 | {
|
---|
| 409 | ++fTowerTrackHits;
|
---|
| 410 |
|
---|
| 411 | track = static_cast<Candidate*>(fTrackInputArray->At(number));
|
---|
| 412 | momentum = track->Momentum;
|
---|
| 413 | position = track->Position;
|
---|
| 414 |
|
---|
| 415 | ecalEnergy = momentum.E() * fECalTrackFractions[number];
|
---|
| 416 | hcalEnergy = momentum.E() * fHCalTrackFractions[number];
|
---|
[a1b19ea] | 417 | energy = ecalEnergy + hcalEnergy;
|
---|
[ef8a06d] | 418 |
|
---|
| 419 | if(ecalEnergy > fTimingEnergyMin && fTower)
|
---|
| 420 | {
|
---|
| 421 | if(fElectronsFromTrack)
|
---|
| 422 | {
|
---|
| 423 | fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
|
---|
| 424 | }
|
---|
| 425 | }
|
---|
| 426 |
|
---|
| 427 |
|
---|
[a1b19ea] | 428 | // in Dual Readout we do not care if tracks are ECAL of HCAL
|
---|
| 429 | if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
|
---|
[fd4b326] | 430 | {
|
---|
[a1b19ea] | 431 | fTrackEnergy += energy;
|
---|
[fd4b326] | 432 | // 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)
|
---|
[a1b19ea] | 433 | sigma = 0.0;
|
---|
| 434 | if(fHCalTrackFractions[number] > 0)
|
---|
| 435 | sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
|
---|
| 436 | else
|
---|
| 437 | sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
|
---|
[fd4b326] | 438 |
|
---|
| 439 | if(sigma/momentum.E() < track->TrackResolution)
|
---|
| 440 | energyGuess = ecalEnergy + hcalEnergy;
|
---|
[a1b19ea] | 441 | else
|
---|
| 442 | energyGuess = momentum.E();
|
---|
[fd4b326] | 443 |
|
---|
[a1b19ea] | 444 | fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
|
---|
| 445 | fTowerTrackArray->Add(track);
|
---|
[fd4b326] | 446 |
|
---|
[a1b19ea] | 447 | }
|
---|
| 448 | else
|
---|
| 449 | {
|
---|
| 450 | fEFlowTrackOutputArray->Add(track);
|
---|
| 451 | }
|
---|
[ef8a06d] | 452 |
|
---|
| 453 | continue;
|
---|
| 454 | }
|
---|
| 455 |
|
---|
| 456 | // check for photon and electron hits in current tower
|
---|
| 457 | if(flags & 2) ++fTowerPhotonHits;
|
---|
| 458 |
|
---|
| 459 | particle = static_cast<Candidate*>(fParticleInputArray->At(number));
|
---|
| 460 | momentum = particle->Momentum;
|
---|
| 461 | position = particle->Position;
|
---|
| 462 |
|
---|
[61dccd3] | 463 |
|
---|
[ef8a06d] | 464 | // fill current tower
|
---|
| 465 | ecalEnergy = momentum.E() * fECalTowerFractions[number];
|
---|
| 466 | hcalEnergy = momentum.E() * fHCalTowerFractions[number];
|
---|
| 467 |
|
---|
| 468 | fECalTowerEnergy += ecalEnergy;
|
---|
| 469 | fHCalTowerEnergy += hcalEnergy;
|
---|
| 470 |
|
---|
[61dccd3] | 471 | // assume combined timing measurements in ECAL/HCAL sections
|
---|
| 472 | fTowerTime += (ecalEnergy + hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
|
---|
| 473 | fTowerTimeWeight += ecalEnergy + hcalEnergy;
|
---|
| 474 | //fTowerTime += (hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
|
---|
| 475 | //fTowerTimeWeight += hcalEnergy;
|
---|
| 476 | //fTowerTime += position.T(); //sigma_t ~ 1/sqrt(E)
|
---|
| 477 | //fTowerTimeWeight += 1;
|
---|
| 478 |
|
---|
| 479 | //cout<<" tower particle PID, pt, eta, phi, l, tof: "<<particle->PID<<", "<<momentum.E()<<", "<<momentum.Eta()<<", "<<momentum.Phi()<<", "<<position.Vect().Mag()<<", "<<position.T()/2.99792458E2<<endl;
|
---|
| 480 | //cout<<" tower particle time, weight: "<<fTowerTime/2.99792458E2<<", "<<fTowerTimeWeight<<endl;
|
---|
[ef8a06d] | 481 |
|
---|
| 482 | fTower->AddCandidate(particle);
|
---|
[d1678fd] | 483 | fTower->Position = position;
|
---|
[ef8a06d] | 484 | }
|
---|
| 485 |
|
---|
| 486 | // finalize last tower
|
---|
| 487 | FinalizeTower();
|
---|
| 488 | }
|
---|
| 489 |
|
---|
| 490 | //------------------------------------------------------------------------------
|
---|
| 491 |
|
---|
| 492 | void DualReadoutCalorimeter::FinalizeTower()
|
---|
| 493 | {
|
---|
| 494 |
|
---|
| 495 | Candidate *track, *tower, *mother;
|
---|
[61dccd3] | 496 | Double_t energy, pt, eta, phi, r, time;
|
---|
[ef8a06d] | 497 | Double_t ecalEnergy, hcalEnergy;
|
---|
[a1b19ea] | 498 | Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
|
---|
[fd4b326] | 499 |
|
---|
[ef8a06d] | 500 | Double_t ecalSigma, hcalSigma, sigma;
|
---|
[a1b19ea] | 501 | Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
|
---|
[ef8a06d] | 502 |
|
---|
| 503 | Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
|
---|
[fd4b326] | 504 |
|
---|
[ef8a06d] | 505 | TLorentzVector momentum;
|
---|
| 506 | TFractionMap::iterator itFractionMap;
|
---|
| 507 |
|
---|
| 508 | Float_t weight, sumWeightedTime, sumWeight;
|
---|
| 509 |
|
---|
| 510 | if(!fTower) return;
|
---|
| 511 |
|
---|
| 512 |
|
---|
[fd4b326] | 513 | // if no hadronic energy, use ECAL resolution
|
---|
[ef8a06d] | 514 | if (fHCalTowerEnergy <= fHCalEnergyMin)
|
---|
| 515 | {
|
---|
| 516 | energy = fECalTowerEnergy;
|
---|
| 517 | sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
|
---|
[a1b19ea] | 518 | //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
|
---|
[ef8a06d] | 519 | }
|
---|
| 520 |
|
---|
[fd4b326] | 521 | // if hadronic fraction > 0, use HCAL resolution
|
---|
[ef8a06d] | 522 | else
|
---|
| 523 | {
|
---|
| 524 | energy = fECalTowerEnergy + fHCalTowerEnergy;
|
---|
| 525 | sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
|
---|
[a1b19ea] | 526 | //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
|
---|
[ef8a06d] | 527 | }
|
---|
| 528 |
|
---|
| 529 | energy = LogNormal(energy, sigma);
|
---|
| 530 | //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
|
---|
[fd4b326] | 531 |
|
---|
[ef8a06d] | 532 | if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
|
---|
| 533 |
|
---|
| 534 | // for now keep this the same
|
---|
| 535 | ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
|
---|
| 536 | hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
|
---|
| 537 |
|
---|
| 538 | ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
|
---|
| 539 | hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
|
---|
| 540 |
|
---|
[61dccd3] | 541 | time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
|
---|
| 542 |
|
---|
[ef8a06d] | 543 | ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
|
---|
| 544 | hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
|
---|
| 545 |
|
---|
| 546 | if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
|
---|
| 547 | if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
|
---|
| 548 |
|
---|
| 549 | if(fSmearTowerCenter)
|
---|
| 550 | {
|
---|
| 551 | eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
|
---|
| 552 | phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
|
---|
| 553 | }
|
---|
| 554 | else
|
---|
| 555 | {
|
---|
| 556 | eta = fTowerEta;
|
---|
| 557 | phi = fTowerPhi;
|
---|
| 558 | }
|
---|
| 559 |
|
---|
| 560 | pt = energy / TMath::CosH(eta);
|
---|
| 561 |
|
---|
[5eda6767] | 562 | // check whether barrel or endcap tower
|
---|
[61dccd3] | 563 | if ((fTowerRmax - fTower->Position.Perp()) < 1.e-06 && TMath::Abs(eta) > 0.)
|
---|
[5eda6767] | 564 | r = fTower->Position.Z()/TMath::SinH(eta);
|
---|
| 565 | else
|
---|
| 566 | r = fTower->Position.Pt();
|
---|
[d1678fd] | 567 |
|
---|
[61dccd3] | 568 | fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
|
---|
| 569 |
|
---|
| 570 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
| 571 | fTower->L = fTower->Position.Vect().Mag();
|
---|
| 572 | //cout<<" tower pt, eta, phi, l, tof: "<<fTower->Momentum.E()<<", "<<fTower->Momentum.Eta()<<", "<<fTower->Momentum.Phi()<<", "<<fTower->L<<", "<<fTower->Position.T()/2.99792458E2<<endl;
|
---|
[ef8a06d] | 573 |
|
---|
| 574 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
|
---|
| 575 | fTower->Eem = ecalEnergy;
|
---|
| 576 | fTower->Ehad = hcalEnergy;
|
---|
[61dccd3] | 577 | fTower->Etrk = fTrackEnergy;
|
---|
[ef8a06d] | 578 | fTower->Edges[0] = fTowerEdges[0];
|
---|
| 579 | fTower->Edges[1] = fTowerEdges[1];
|
---|
| 580 | fTower->Edges[2] = fTowerEdges[2];
|
---|
| 581 | fTower->Edges[3] = fTowerEdges[3];
|
---|
| 582 |
|
---|
| 583 | if(energy > 0.0)
|
---|
| 584 | {
|
---|
| 585 | if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
|
---|
| 586 | {
|
---|
| 587 | fPhotonOutputArray->Add(fTower);
|
---|
| 588 | }
|
---|
| 589 | fTowerOutputArray->Add(fTower);
|
---|
| 590 | }
|
---|
[a1b19ea] | 591 |
|
---|
| 592 | // fill energy flow candidates
|
---|
[fd4b326] | 593 |
|
---|
[a1b19ea] | 594 | fTrackSigma = TMath::Sqrt(fTrackSigma);
|
---|
| 595 | neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
|
---|
| 596 | neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
|
---|
| 597 |
|
---|
| 598 | if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
|
---|
| 599 | {
|
---|
| 600 | // create new photon tower
|
---|
| 601 | tower = static_cast<Candidate*>(fTower->Clone());
|
---|
[4afb18d] | 602 | pt = neutralEnergy / TMath::CosH(eta);
|
---|
[a1b19ea] | 603 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
|
---|
| 604 |
|
---|
[fd4b326] | 605 | // if no hadronic energy, use ECAL resolution
|
---|
[f382b7e] | 606 | if (fHCalTowerEnergy <= fHCalEnergyMin)
|
---|
| 607 | {
|
---|
| 608 | tower->Eem = neutralEnergy;
|
---|
| 609 | tower->Ehad = 0.0;
|
---|
| 610 | tower->PID = 22;
|
---|
[4afb18d] | 611 | fEFlowPhotonOutputArray->Add(tower);
|
---|
[f382b7e] | 612 | }
|
---|
[fd4b326] | 613 | // if hadronic fraction > 0, use HCAL resolution
|
---|
[f382b7e] | 614 | else
|
---|
| 615 | {
|
---|
| 616 | tower->Eem = 0;
|
---|
| 617 | tower->Ehad = neutralEnergy;
|
---|
| 618 | tower->PID = 130;
|
---|
[4afb18d] | 619 | fEFlowNeutralHadronOutputArray->Add(tower);
|
---|
[f382b7e] | 620 | }
|
---|
[a1b19ea] | 621 |
|
---|
| 622 | //clone tracks
|
---|
| 623 | fItTowerTrackArray->Reset();
|
---|
| 624 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
|
---|
| 625 | {
|
---|
| 626 | mother = track;
|
---|
| 627 | track = static_cast<Candidate*>(track->Clone());
|
---|
| 628 | track->AddCandidate(mother);
|
---|
| 629 | fEFlowTrackOutputArray->Add(track);
|
---|
| 630 | }
|
---|
| 631 | }
|
---|
| 632 |
|
---|
[fd4b326] | 633 |
|
---|
[a1b19ea] | 634 | // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
|
---|
| 635 | else if(fTrackEnergy > 0.0)
|
---|
| 636 | {
|
---|
| 637 | //cout<<"no significant neutral excess found:"<<endl;
|
---|
| 638 | weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
|
---|
| 639 | weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
|
---|
| 640 |
|
---|
[fd4b326] | 641 | bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
|
---|
[a1b19ea] | 642 | rescaleFactor = bestEnergyEstimate/fTrackEnergy;
|
---|
| 643 |
|
---|
| 644 | //rescale tracks
|
---|
| 645 | fItTowerTrackArray->Reset();
|
---|
| 646 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
|
---|
[fd4b326] | 647 | {
|
---|
[a1b19ea] | 648 | mother = track;
|
---|
[fd4b326] | 649 | track = static_cast<Candidate *>(track->Clone());
|
---|
[a1b19ea] | 650 | track->AddCandidate(mother);
|
---|
[fd4b326] | 651 | track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
|
---|
[a1b19ea] | 652 | fEFlowTrackOutputArray->Add(track);
|
---|
| 653 | }
|
---|
| 654 | }
|
---|
[fd4b326] | 655 |
|
---|
[a1b19ea] | 656 |
|
---|
[ef8a06d] | 657 | }
|
---|
| 658 |
|
---|
| 659 | //------------------------------------------------------------------------------
|
---|
| 660 |
|
---|
| 661 | Double_t DualReadoutCalorimeter::LogNormal(Double_t mean, Double_t sigma)
|
---|
| 662 | {
|
---|
| 663 | Double_t a, b;
|
---|
| 664 |
|
---|
| 665 | if(mean > 0.0)
|
---|
| 666 | {
|
---|
| 667 | b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
|
---|
| 668 | a = TMath::Log(mean) - 0.5*b*b;
|
---|
| 669 |
|
---|
| 670 | return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
|
---|
| 671 | }
|
---|
| 672 | else
|
---|
| 673 | {
|
---|
| 674 | return 0.0;
|
---|
| 675 | }
|
---|
| 676 | }
|
---|