Changeset 1364 in svn
- Timestamp:
- Apr 16, 2014, 3:29:31 PM (11 years ago)
- Location:
- trunk/modules
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/modules/Calorimeter.cc
r1356 r1364 41 41 42 42 Calorimeter::Calorimeter() : 43 f ECalResolutionFormula(0), fHCalResolutionFormula(0),43 fResolutionFormula(0), 44 44 fItParticleInputArray(0), fItTrackInputArray(0), 45 45 fTowerTrackArray(0), fItTowerTrackArray(0) 46 46 { 47 fECalResolutionFormula = new DelphesFormula; 48 fHCalResolutionFormula = new DelphesFormula; 49 47 fResolutionFormula = new DelphesFormula; 48 50 49 fTowerTrackArray = new TObjArray; 51 50 fItTowerTrackArray = fTowerTrackArray->MakeIterator(); … … 56 55 Calorimeter::~Calorimeter() 57 56 { 58 if(fECalResolutionFormula) delete fECalResolutionFormula; 59 if(fHCalResolutionFormula) delete fHCalResolutionFormula; 60 57 if(fResolutionFormula) delete fResolutionFormula; 58 61 59 if(fTowerTrackArray) delete fTowerTrackArray; 62 60 if(fItTowerTrackArray) delete fItTowerTrackArray; … … 68 66 { 69 67 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions; 70 Long_t i, j, k, size, sizeEtaBins, sizePhiBins , sizeFractions;71 Double_t ecalFraction, hcalFraction;68 Long_t i, j, k, size, sizeEtaBins, sizePhiBins; 69 Double_t fraction; 72 70 TBinMap::iterator itEtaBin; 73 71 set< Double_t >::iterator itPhiBin; … … 116 114 // set default energy fractions values 117 115 fFractionMap.clear(); 118 fFractionMap[0] = make_pair(0.0, 1.0);116 fFractionMap[0] = 1.0; 119 117 120 118 for(i = 0; i < size/2; ++i) 121 119 { 122 120 paramFractions = param[i*2 + 1]; 123 sizeFractions = paramFractions.GetSize(); 124 125 ecalFraction = paramFractions[0].GetDouble(); 126 hcalFraction = paramFractions[1].GetDouble(); 127 128 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction); 121 fraction = paramFractions[0].GetDouble(); 122 fFractionMap[param[i*2].GetInt()] = fraction; 129 123 } 130 124 /* … … 136 130 */ 137 131 // read resolution formulas 138 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0")); 139 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0")); 140 132 fResolutionFormula->Compile(GetString("ResolutionFormula", "0")); 133 141 134 // import array with output from other modules 142 135 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles")); … … 148 141 // create output arrays 149 142 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers")); 150 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons")); 151 152 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks")); 153 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons")); 154 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons")); 155 156 143 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers")); 144 157 145 } 158 146 … … 179 167 Int_t number; 180 168 Long64_t towerHit, towerEtaPhi, hitEtaPhi; 181 Double_t ecalFraction, hcalFraction;182 Double_t e calEnergy, hcalEnergy;169 Double_t fraction; 170 Double_t energy; 183 171 Int_t pdgCode; 184 172 … … 193 181 DelphesFactory *factory = GetFactory(); 194 182 fTowerHits.clear(); 195 fTowerECalFractions.clear(); 196 fTowerHCalFractions.clear(); 197 fTrackECalFractions.clear(); 198 fTrackHCalFractions.clear(); 199 183 fTowerFractions.clear(); 184 fTrackFractions.clear(); 185 200 186 // loop over all particles 201 187 fItParticleInputArray->Reset(); … … 214 200 } 215 201 216 ecalFraction = itFractionMap->second.first; 217 hcalFraction = itFractionMap->second.second; 218 219 fTowerECalFractions.push_back(ecalFraction); 220 fTowerHCalFractions.push_back(hcalFraction); 221 222 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue; 202 fraction = itFractionMap->second; 203 fTowerFractions.push_back(fraction); 204 205 if(fraction < 1.0E-9) continue; 223 206 224 207 // find eta bin [1, fEtaBins.size - 1] … … 260 243 } 261 244 262 ecalFraction = itFractionMap->second.first; 263 hcalFraction = itFractionMap->second.second; 264 265 fTrackECalFractions.push_back(ecalFraction); 266 fTrackHCalFractions.push_back(hcalFraction); 267 245 fraction = itFractionMap->second; 246 247 fTrackFractions.push_back(fraction); 248 268 249 // find eta bin [1, fEtaBins.size - 1] 269 250 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta()); … … 327 308 fTowerEdges[3] = (*phiBins)[phiBin]; 328 309 329 fTowerECalEnergy = 0.0; 330 fTowerHCalEnergy = 0.0; 331 332 fTrackECalEnergy = 0.0; 333 fTrackHCalEnergy = 0.0; 334 335 fTowerECalTime = 0.0; 336 fTowerHCalTime = 0.0; 337 338 fTrackECalTime = 0.0; 339 fTrackHCalTime = 0.0; 340 341 fTowerECalWeightTime = 0.0; 342 fTowerHCalWeightTime = 0.0; 310 fTowerEnergy = 0.0; 311 fTrackEnergy = 0.0; 312 313 fTowerTime = 0.0; 314 fTrackTime = 0.0; 343 315 316 fTowerWeightTime = 0.0; 317 344 318 fTowerTrackHits = 0; 345 319 fTowerPhotonHits = 0; … … 356 330 momentum = track->Momentum; 357 331 position = track->Position; 358 332 333 energy = momentum.E() * fTrackFractions[number]; 359 334 360 ecalEnergy = momentum.E() * fTrackECalFractions[number]; 361 hcalEnergy = momentum.E() * fTrackHCalFractions[number]; 362 363 fTrackECalEnergy += ecalEnergy; 364 fTrackHCalEnergy += hcalEnergy; 335 fTrackEnergy += energy; 365 336 366 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 367 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 368 369 fTrackECalWeightTime += TMath::Sqrt(ecalEnergy); 370 fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy); 371 337 fTrackTime += TMath::Sqrt(energy)*position.T(); 338 fTrackWeightTime += TMath::Sqrt(energy); 339 372 340 fTowerTrackArray->Add(track); 373 341 … … 383 351 384 352 // fill current tower 385 ecalEnergy = momentum.E() * fTowerECalFractions[number]; 386 hcalEnergy = momentum.E() * fTowerHCalFractions[number]; 387 388 fTowerECalEnergy += ecalEnergy; 389 fTowerHCalEnergy += hcalEnergy; 390 391 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 392 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 393 394 fTowerECalWeightTime += TMath::Sqrt(ecalEnergy); 395 fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy); 353 energy = momentum.E() * fTowerFractions[number]; 354 355 fTowerEnergy += energy; 396 356 397 357 fTowerTime += TMath::Sqrt(energy)*position.T(); 358 fTowerWeightTime += TMath::Sqrt(energy); 359 398 360 fTower->AddCandidate(particle); 399 361 } … … 407 369 void Calorimeter::FinalizeTower() 408 370 { 409 Candidate *t rack, *tower;371 Candidate *tower; 410 372 Double_t energy, pt, eta, phi; 411 Double_t ecalEnergy, hcalEnergy; 412 Double_t ecalSigma, hcalSigma; 413 Double_t ecalTime, hcalTime, time; 373 Double_t sigma; 374 Double_t time; 414 375 415 376 if(!fTower) return; 416 377 417 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy); 418 419 // ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma); 420 // if(ecalEnergy < 0.0) ecalEnergy = 0.0; 421 422 ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma); 423 ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime; 424 425 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy); 426 427 // hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma); 428 // if(hcalEnergy < 0.0) hcalEnergy = 0.0; 429 430 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma); 431 hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime; 432 433 energy = ecalEnergy + hcalEnergy; 434 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 435 436 // eta = fTowerEta; 437 // phi = fTowerPhi; 378 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy); 379 380 // energy = gRandom->Gaus(fTowerEnergy, sigma); 381 // if(energy < 0.0) energy = 0.0; 382 383 energy = LogNormal(fTowerEnergy, sigma); 384 time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime; 438 385 439 386 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]); … … 445 392 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time); 446 393 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 447 fTower->Eem = ecalEnergy; 448 fTower->Ehad = hcalEnergy; 449 394 450 395 fTower->Edges[0] = fTowerEdges[0]; 451 396 fTower->Edges[1] = fTowerEdges[1]; … … 455 400 456 401 // fill calorimeter towers 402 if(energy > 0.0) fTowerOutputArray->Add(fTower); 403 404 405 // fill energy flow candidates 406 energy -= fTrackEnergy; 407 if(energy < 0.0) energy = 0.0; 408 409 // save energy excess as an energy flow tower 457 410 if(energy > 0.0) 458 {459 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)460 {461 fPhotonOutputArray->Add(fTower);462 }463 464 fTowerOutputArray->Add(fTower);465 }466 467 // fill energy flow candidates468 469 // save all the tracks as energy flow tracks470 fItTowerTrackArray->Reset();471 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))472 {473 fEFlowTrackOutputArray->Add(track);474 }475 476 ecalEnergy -= fTrackECalEnergy;477 if(ecalEnergy < 0.0) ecalEnergy = 0.0;478 479 hcalEnergy -= fTrackHCalEnergy;480 if(hcalEnergy < 0.0) hcalEnergy = 0.0;481 482 energy = ecalEnergy + hcalEnergy;483 484 485 // save ECAL and/or HCAL energy excess as an energy flow tower486 if(ecalEnergy > 0.0)487 411 { 488 412 // create new photon tower 489 413 tower = static_cast<Candidate*>(fTower->Clone()); 490 491 pt = ecalEnergy / TMath::CosH(eta); 492 493 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy); 494 tower->Eem = ecalEnergy; 495 tower->Ehad = 0; 496 497 fEFlowPhotonOutputArray->Add(tower); 498 } 499 500 if(hcalEnergy > 0.0) 501 { 502 // create new neutral hadron tower 503 tower = static_cast<Candidate*>(fTower->Clone()); 504 505 pt = hcalEnergy / TMath::CosH(eta); 506 507 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy); 508 tower->Eem = 0; 509 tower->Ehad = hcalEnergy; 510 511 fEFlowNeutralHadronOutputArray->Add(tower); 512 } 513 514 515 414 pt = energy / TMath::CosH(eta); 415 416 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 417 fEFlowTowerOutputArray->Add(tower); 418 } 516 419 517 420 } -
trunk/modules/Calorimeter.h
r1356 r1364 38 38 private: 39 39 40 typedef std::map< Long64_t, std::pair< Double_t, Double_t >> TFractionMap; //!40 typedef std::map< Long64_t, Double_t > TFractionMap; //! 41 41 typedef std::map< Double_t, std::set< Double_t > > TBinMap; //! 42 42 43 43 Candidate *fTower; 44 44 Double_t fTowerEta, fTowerPhi, fTowerEdges[4]; 45 Double_t fTowerE CalEnergy, fTowerHCalEnergy;46 Double_t fTrackE CalEnergy, fTrackHCalEnergy;45 Double_t fTowerEnergy; 46 Double_t fTrackEnergy; 47 47 48 Double_t fTower ECalTime, fTowerHCalTime;49 Double_t fTrack ECalTime, fTrackHCalTime;48 Double_t fTowerTime; 49 Double_t fTrackTime; 50 50 51 Double_t fTower ECalWeightTime, fTowerHCalWeightTime;52 Double_t fTrack ECalWeightTime, fTrackHCalWeightTime;51 Double_t fTowerWeightTime; 52 Double_t fTrackWeightTime; 53 53 54 54 Int_t fTowerTrackHits, fTowerPhotonHits; … … 62 62 std::vector < Long64_t > fTowerHits; 63 63 64 std::vector < Double_t > fTowerECalFractions; 65 std::vector < Double_t > fTowerHCalFractions; 66 67 std::vector < Double_t > fTrackECalFractions; 68 std::vector < Double_t > fTrackHCalFractions; 69 70 DelphesFormula *fECalResolutionFormula; //! 71 DelphesFormula *fHCalResolutionFormula; //! 72 64 std::vector < Double_t > fTowerFractions; 65 66 std::vector < Double_t > fTrackFractions; 67 68 DelphesFormula *fResolutionFormula; //! 69 73 70 TIterator *fItParticleInputArray; //! 74 71 TIterator *fItTrackInputArray; //! … … 78 75 79 76 TObjArray *fTowerOutputArray; //! 80 TObjArray *fPhotonOutputArray; //!81 77 82 TObjArray *fEFlowTrackOutputArray; //! 83 TObjArray *fEFlowPhotonOutputArray; //! 84 TObjArray *fEFlowNeutralHadronOutputArray; //! 85 78 TObjArray *fEFlowTowerOutputArray; //! 79 86 80 TObjArray *fTowerTrackArray; //! 87 81 TIterator *fItTowerTrackArray; //!
Note:
See TracChangeset
for help on using the changeset viewer.