- Timestamp:
- Apr 16, 2014, 5:17:35 PM (11 years ago)
- Location:
- trunk
- Files:
-
- 2 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Makefile
r1368 r1369 316 316 modules/ImpactParameterSmearing.h \ 317 317 modules/TimeSmearing.h \ 318 modules/SimpleCalorimeter.h \ 318 319 modules/Calorimeter.h \ 319 320 modules/Isolation.h \ … … 621 622 modules/EnergyScale.$(SrcSuf) \ 622 623 modules/EnergyScale.h \ 624 classes/DelphesClasses.h \ 625 classes/DelphesFactory.h \ 626 classes/DelphesFormula.h \ 627 external/ExRootAnalysis/ExRootResult.h \ 628 external/ExRootAnalysis/ExRootFilter.h \ 629 external/ExRootAnalysis/ExRootClassifier.h 630 tmp/modules/SimpleCalorimeter.$(ObjSuf): \ 631 modules/SimpleCalorimeter.$(SrcSuf) \ 632 modules/SimpleCalorimeter.h \ 623 633 classes/DelphesClasses.h \ 624 634 classes/DelphesFactory.h \ … … 1067 1077 tmp/modules/ExampleModule.$(ObjSuf) \ 1068 1078 tmp/modules/EnergyScale.$(ObjSuf) \ 1079 tmp/modules/SimpleCalorimeter.$(ObjSuf) \ 1069 1080 tmp/modules/ParticlePropagator.$(ObjSuf) \ 1070 1081 tmp/modules/Weighter.$(ObjSuf) \ … … 1558 1569 @touch $@ 1559 1570 1571 modules/SimpleCalorimeter.h: \ 1572 classes/DelphesModule.h 1573 @touch $@ 1574 1560 1575 external/fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh: \ 1561 1576 external/fastjet/JetDefinition.hh \ -
trunk/modules/Calorimeter.cc
r1364 r1369 41 41 42 42 Calorimeter::Calorimeter() : 43 f ResolutionFormula(0),43 fECalResolutionFormula(0), fHCalResolutionFormula(0), 44 44 fItParticleInputArray(0), fItTrackInputArray(0), 45 45 fTowerTrackArray(0), fItTowerTrackArray(0) 46 46 { 47 fResolutionFormula = new DelphesFormula; 48 47 fECalResolutionFormula = new DelphesFormula; 48 fHCalResolutionFormula = new DelphesFormula; 49 49 50 fTowerTrackArray = new TObjArray; 50 51 fItTowerTrackArray = fTowerTrackArray->MakeIterator(); … … 55 56 Calorimeter::~Calorimeter() 56 57 { 57 if(fResolutionFormula) delete fResolutionFormula; 58 58 if(fECalResolutionFormula) delete fECalResolutionFormula; 59 if(fHCalResolutionFormula) delete fHCalResolutionFormula; 60 59 61 if(fTowerTrackArray) delete fTowerTrackArray; 60 62 if(fItTowerTrackArray) delete fItTowerTrackArray; … … 66 68 { 67 69 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions; 68 Long_t i, j, k, size, sizeEtaBins, sizePhiBins ;69 Double_t fraction;70 Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions; 71 Double_t ecalFraction, hcalFraction; 70 72 TBinMap::iterator itEtaBin; 71 73 set< Double_t >::iterator itPhiBin; … … 114 116 // set default energy fractions values 115 117 fFractionMap.clear(); 116 fFractionMap[0] = 1.0;118 fFractionMap[0] = make_pair(0.0, 1.0); 117 119 118 120 for(i = 0; i < size/2; ++i) 119 121 { 120 122 paramFractions = param[i*2 + 1]; 121 fraction = paramFractions[0].GetDouble(); 122 fFractionMap[param[i*2].GetInt()] = fraction; 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); 123 129 } 124 130 /* … … 130 136 */ 131 137 // read resolution formulas 132 fResolutionFormula->Compile(GetString("ResolutionFormula", "0")); 133 138 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0")); 139 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0")); 140 134 141 // import array with output from other modules 135 142 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles")); … … 141 148 // create output arrays 142 149 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers")); 143 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers")); 144 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 145 157 } 146 158 … … 167 179 Int_t number; 168 180 Long64_t towerHit, towerEtaPhi, hitEtaPhi; 169 Double_t fraction;170 Double_t e nergy;181 Double_t ecalFraction, hcalFraction; 182 Double_t ecalEnergy, hcalEnergy; 171 183 Int_t pdgCode; 172 184 … … 181 193 DelphesFactory *factory = GetFactory(); 182 194 fTowerHits.clear(); 183 fTowerFractions.clear(); 184 fTrackFractions.clear(); 185 195 fTowerECalFractions.clear(); 196 fTowerHCalFractions.clear(); 197 fTrackECalFractions.clear(); 198 fTrackHCalFractions.clear(); 199 186 200 // loop over all particles 187 201 fItParticleInputArray->Reset(); … … 200 214 } 201 215 202 fraction = itFractionMap->second; 203 fTowerFractions.push_back(fraction); 204 205 if(fraction < 1.0E-9) continue; 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; 206 223 207 224 // find eta bin [1, fEtaBins.size - 1] … … 243 260 } 244 261 245 fraction = itFractionMap->second; 246 247 fTrackFractions.push_back(fraction); 248 262 ecalFraction = itFractionMap->second.first; 263 hcalFraction = itFractionMap->second.second; 264 265 fTrackECalFractions.push_back(ecalFraction); 266 fTrackHCalFractions.push_back(hcalFraction); 267 249 268 // find eta bin [1, fEtaBins.size - 1] 250 269 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta()); … … 308 327 fTowerEdges[3] = (*phiBins)[phiBin]; 309 328 310 fTowerEnergy = 0.0; 311 fTrackEnergy = 0.0; 312 313 fTowerTime = 0.0; 314 fTrackTime = 0.0; 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; 315 343 316 fTowerWeightTime = 0.0;317 318 344 fTowerTrackHits = 0; 319 345 fTowerPhotonHits = 0; … … 330 356 momentum = track->Momentum; 331 357 position = track->Position; 332 333 energy = momentum.E() * fTrackFractions[number]; 358 334 359 335 fTrackEnergy += energy; 360 ecalEnergy = momentum.E() * fTrackECalFractions[number]; 361 hcalEnergy = momentum.E() * fTrackHCalFractions[number]; 362 363 fTrackECalEnergy += ecalEnergy; 364 fTrackHCalEnergy += hcalEnergy; 336 365 337 fTrackTime += TMath::Sqrt(energy)*position.T(); 338 fTrackWeightTime += TMath::Sqrt(energy); 339 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 340 372 fTowerTrackArray->Add(track); 341 373 … … 351 383 352 384 // fill current tower 353 energy = momentum.E() * fTowerFractions[number]; 354 355 fTowerEnergy += energy; 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); 356 396 357 fTowerTime += TMath::Sqrt(energy)*position.T(); 358 fTowerWeightTime += TMath::Sqrt(energy); 359 397 360 398 fTower->AddCandidate(particle); 361 399 } … … 369 407 void Calorimeter::FinalizeTower() 370 408 { 371 Candidate *t ower;409 Candidate *track, *tower; 372 410 Double_t energy, pt, eta, phi; 373 Double_t sigma; 374 Double_t time; 411 Double_t ecalEnergy, hcalEnergy; 412 Double_t ecalSigma, hcalSigma; 413 Double_t ecalTime, hcalTime, time; 375 414 376 415 if(!fTower) return; 377 416 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; 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; 385 438 386 439 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]); … … 392 445 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time); 393 446 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 394 447 fTower->Eem = ecalEnergy; 448 fTower->Ehad = hcalEnergy; 449 395 450 fTower->Edges[0] = fTowerEdges[0]; 396 451 fTower->Edges[1] = fTowerEdges[1]; … … 400 455 401 456 // fill calorimeter towers 402 if(energy > 0.0) fTowerOutputArray->Add(fTower); 457 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 candidates 468 469 // save all the tracks as energy flow tracks 470 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; 403 483 404 484 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 410 if(energy > 0.0) 485 // save ECAL and/or HCAL energy excess as an energy flow tower 486 if(ecalEnergy > 0.0) 411 487 { 412 488 // create new photon tower 413 489 tower = static_cast<Candidate*>(fTower->Clone()); 414 pt = energy / TMath::CosH(eta); 415 416 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 417 fEFlowTowerOutputArray->Add(tower); 418 } 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 419 516 420 517 } -
trunk/modules/Calorimeter.h
r1364 r1369 38 38 private: 39 39 40 typedef std::map< Long64_t, Double_t> TFractionMap; //!40 typedef std::map< Long64_t, std::pair< Double_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 nergy;46 Double_t fTrackE nergy;45 Double_t fTowerECalEnergy, fTowerHCalEnergy; 46 Double_t fTrackECalEnergy, fTrackHCalEnergy; 47 47 48 Double_t fTower Time;49 Double_t fTrack Time;48 Double_t fTowerECalTime, fTowerHCalTime; 49 Double_t fTrackECalTime, fTrackHCalTime; 50 50 51 Double_t fTower WeightTime;52 Double_t fTrack WeightTime;51 Double_t fTowerECalWeightTime, fTowerHCalWeightTime; 52 Double_t fTrackECalWeightTime, fTrackHCalWeightTime; 53 53 54 54 Int_t fTowerTrackHits, fTowerPhotonHits; … … 62 62 std::vector < Long64_t > fTowerHits; 63 63 64 std::vector < Double_t > fTowerFractions; 65 66 std::vector < Double_t > fTrackFractions; 67 68 DelphesFormula *fResolutionFormula; //! 69 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 70 73 TIterator *fItParticleInputArray; //! 71 74 TIterator *fItTrackInputArray; //! … … 75 78 76 79 TObjArray *fTowerOutputArray; //! 80 TObjArray *fPhotonOutputArray; //! 77 81 78 TObjArray *fEFlowTowerOutputArray; //! 79 82 TObjArray *fEFlowTrackOutputArray; //! 83 TObjArray *fEFlowPhotonOutputArray; //! 84 TObjArray *fEFlowNeutralHadronOutputArray; //! 85 80 86 TObjArray *fTowerTrackArray; //! 81 87 TIterator *fItTowerTrackArray; //! -
trunk/modules/ModulesLinkDef.h
r1367 r1369 21 21 #include "modules/ImpactParameterSmearing.h" 22 22 #include "modules/TimeSmearing.h" 23 #include "modules/SimpleCalorimeter.h" 23 24 #include "modules/Calorimeter.h" 24 25 #include "modules/Isolation.h" … … 57 58 #pragma link C++ class ImpactParameterSmearing+; 58 59 #pragma link C++ class TimeSmearing+; 60 #pragma link C++ class SimpleCalorimeter+; 59 61 #pragma link C++ class Calorimeter+; 60 62 #pragma link C++ class Isolation+;
Note:
See TracChangeset
for help on using the changeset viewer.