Changeset 82575a3 in git for modules/Calorimeter.cc
- Timestamp:
- Apr 16, 2014, 5:17:35 PM (11 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- ab16005
- Parents:
- 9687203
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/Calorimeter.cc
r9687203 r82575a3 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 }
Note:
See TracChangeset
for help on using the changeset viewer.