Changeset 341014c in git for modules/SimpleCalorimeter.cc
- Timestamp:
- Feb 12, 2019, 9:29:17 PM (6 years ago)
- Branches:
- ImprovedOutputFile, Timing, llp, master
- Children:
- 6455202
- Parents:
- 45e58be
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/SimpleCalorimeter.cc
r45e58be r341014c 17 17 */ 18 18 19 20 19 /** \class SimpleCalorimeter 21 20 * … … 33 32 #include "classes/DelphesFormula.h" 34 33 34 #include "ExRootAnalysis/ExRootClassifier.h" 35 #include "ExRootAnalysis/ExRootFilter.h" 35 36 #include "ExRootAnalysis/ExRootResult.h" 36 #include "ExRootAnalysis/ExRootFilter.h" 37 #include "ExRootAnalysis/ExRootClassifier.h" 38 37 38 #include "TDatabasePDG.h" 39 #include "TFormula.h" 40 #include "TLorentzVector.h" 39 41 #include "TMath.h" 42 #include "TObjArray.h" 43 #include "TRandom3.h" 40 44 #include "TString.h" 41 #include "TFormula.h"42 #include "TRandom3.h"43 #include "TObjArray.h"44 #include "TDatabasePDG.h"45 #include "TLorentzVector.h"46 45 47 46 #include <algorithm> 48 #include <stdexcept>49 47 #include <iostream> 50 48 #include <sstream> 49 #include <stdexcept> 51 50 52 51 using namespace std; … … 58 57 fItParticleInputArray(0), fItTrackInputArray(0) 59 58 { 60 59 61 60 fResolutionFormula = new DelphesFormula; 62 61 fTowerTrackArray = new TObjArray; 63 62 fItTowerTrackArray = fTowerTrackArray->MakeIterator(); 64 65 63 } 66 64 … … 69 67 SimpleCalorimeter::~SimpleCalorimeter() 70 68 { 71 69 72 70 if(fResolutionFormula) delete fResolutionFormula; 73 71 if(fTowerTrackArray) delete fTowerTrackArray; 74 72 if(fItTowerTrackArray) delete fItTowerTrackArray; 75 76 73 } 77 74 … … 84 81 Double_t fraction; 85 82 TBinMap::iterator itEtaBin; 86 set< Double_t>::iterator itPhiBin;87 vector< Double_t> *phiBins;83 set<Double_t>::iterator itPhiBin; 84 vector<Double_t> *phiBins; 88 85 89 86 // read eta and phi bins … … 93 90 fEtaBins.clear(); 94 91 fPhiBins.clear(); 95 for(i = 0; i < size /2; ++i)96 { 97 paramEtaBins = param[i *2];92 for(i = 0; i < size / 2; ++i) 93 { 94 paramEtaBins = param[i * 2]; 98 95 sizeEtaBins = paramEtaBins.GetSize(); 99 paramPhiBins = param[i *2 + 1];96 paramPhiBins = param[i * 2 + 1]; 100 97 sizePhiBins = paramPhiBins.GetSize(); 101 98 … … 114 111 { 115 112 fEtaBins.push_back(itEtaBin->first); 116 phiBins = new vector< double>(itEtaBin->second.size());113 phiBins = new vector<double>(itEtaBin->second.size()); 117 114 fPhiBins.push_back(phiBins); 118 115 phiBins->clear(); … … 131 128 fFractionMap[0] = 1.0; 132 129 133 for(i = 0; i < size /2; ++i)134 { 135 paramFractions = param[i *2 + 1];130 for(i = 0; i < size / 2; ++i) 131 { 132 paramFractions = param[i * 2 + 1]; 136 133 fraction = paramFractions[0].GetDouble(); 137 fFractionMap[param[i *2].GetInt()] = fraction;134 fFractionMap[param[i * 2].GetInt()] = fraction; 138 135 } 139 136 … … 170 167 void SimpleCalorimeter::Finish() 171 168 { 172 vector< vector< Double_t >*>::iterator itPhiBin;169 vector<vector<Double_t> *>::iterator itPhiBin; 173 170 if(fItParticleInputArray) delete fItParticleInputArray; 174 171 if(fItTrackInputArray) delete fItTrackInputArray; … … 197 194 TFractionMap::iterator itFractionMap; 198 195 199 vector< Double_t>::iterator itEtaBin;200 vector< Double_t>::iterator itPhiBin;201 vector< Double_t> *phiBins;202 203 vector< Long64_t>::iterator itTowerHits;196 vector<Double_t>::iterator itEtaBin; 197 vector<Double_t>::iterator itPhiBin; 198 vector<Double_t> *phiBins; 199 200 vector<Long64_t>::iterator itTowerHits; 204 201 205 202 DelphesFactory *factory = GetFactory(); … … 211 208 fItParticleInputArray->Reset(); 212 209 number = -1; 213 while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))210 while((particle = static_cast<Candidate *>(fItParticleInputArray->Next()))) 214 211 { 215 212 const TLorentzVector &particlePosition = particle->Position; … … 254 251 fItTrackInputArray->Reset(); 255 252 number = -1; 256 while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))253 while((track = static_cast<Candidate *>(fItTrackInputArray->Next()))) 257 254 { 258 255 const TLorentzVector &trackPosition = track->Position; … … 303 300 towerHit = (*itTowerHits); 304 301 flags = (towerHit >> 24) & 0x00000000000000FFLL; 305 number = (towerHit) &0x0000000000FFFFFFLL;302 number = (towerHit)&0x0000000000FFFFFFLL; 306 303 hitEtaPhi = towerHit >> 32; 307 304 … … 324 321 325 322 // calculate eta and phi of the tower's center 326 fTowerEta = 0.5 *(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);327 fTowerPhi = 0.5 *((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);323 fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]); 324 fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]); 328 325 329 326 fTowerEdges[0] = fEtaBins[etaBin - 1]; … … 346 343 347 344 fTowerTrackArray->Clear(); 348 345 } 349 346 350 347 // check for track hits … … 353 350 ++fTowerTrackHits; 354 351 355 track = static_cast<Candidate *>(fTrackInputArray->At(number));352 track = static_cast<Candidate *>(fTrackInputArray->At(number)); 356 353 momentum = track->Momentum; 357 354 position = track->Position; … … 359 356 energy = momentum.E() * fTrackFractions[number]; 360 357 361 fTrackTime += TMath::Sqrt(energy) *position.T();358 fTrackTime += TMath::Sqrt(energy) * position.T(); 362 359 fTrackTimeWeight += TMath::Sqrt(energy); 363 360 364 361 if(fTrackFractions[number] > 1.0E-9) 365 362 { 366 367 // compute total charged energy 368 fTrackEnergy += energy; 369 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 370 if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy; 371 else energyGuess = momentum.E(); 372 373 fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess); 374 fTowerTrackArray->Add(track); 363 364 // compute total charged energy 365 fTrackEnergy += energy; 366 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 367 if(sigma / momentum.E() < track->TrackResolution) 368 energyGuess = energy; 369 else 370 energyGuess = momentum.E(); 371 372 fTrackSigma += ((track->TrackResolution) * energyGuess) * ((track->TrackResolution) * energyGuess); 373 fTowerTrackArray->Add(track); 375 374 } 376 375 377 376 else 378 377 { … … 386 385 if(flags & 2) ++fTowerPhotonHits; 387 386 388 particle = static_cast<Candidate *>(fParticleInputArray->At(number));387 particle = static_cast<Candidate *>(fParticleInputArray->At(number)); 389 388 momentum = particle->Momentum; 390 389 position = particle->Position; … … 395 394 fTowerEnergy += energy; 396 395 397 fTowerTime += energy *position.T();396 fTowerTime += energy * position.T(); 398 397 fTowerTimeWeight += energy; 399 398 … … 410 409 { 411 410 Candidate *tower, *track, *mother; 412 Double_t energy, neutralEnergy, pt, eta, phi;411 Double_t energy, neutralEnergy, pt, eta, phi; 413 412 Double_t sigma, neutralSigma; 414 413 Double_t time; 415 414 416 415 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor; 417 416 … … 425 424 energy = LogNormal(fTowerEnergy, sigma); 426 425 427 time = (fTowerTimeWeight < 1.0E-09 ) ? 0.0 : fTowerTime/fTowerTimeWeight;426 time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight; 428 427 429 428 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy); 430 429 431 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0; 432 430 if(energy < fEnergyMin || energy < fEnergySignificanceMin * sigma) energy = 0.0; 433 431 434 432 if(fSmearTowerCenter) … … 459 457 if(energy > 0.0) fTowerOutputArray->Add(fTower); 460 458 461 462 459 // e-flow candidates 463 460 464 461 //compute neutral excess 465 462 466 463 fTrackSigma = TMath::Sqrt(fTrackSigma); 467 neutralEnergy = max( (energy - fTrackEnergy), 0.0);468 464 neutralEnergy = max((energy - fTrackEnergy), 0.0); 465 469 466 //compute sigma_trk total 470 neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma *fTrackSigma+ sigma*sigma);471 467 neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma * fTrackSigma + sigma * sigma); 468 472 469 // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack 473 470 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin) 474 471 { 475 472 // create new photon tower 476 tower = static_cast<Candidate *>(fTower->Clone());473 tower = static_cast<Candidate *>(fTower->Clone()); 477 474 pt = neutralEnergy / TMath::CosH(eta); 478 475 … … 480 477 tower->Ehad = (fIsEcal) ? 0 : neutralEnergy; 481 478 tower->PID = (fIsEcal) ? 22 : 0; 482 479 483 480 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy); 484 481 fEFlowTowerOutputArray->Add(tower); 485 482 486 483 fItTowerTrackArray->Reset(); 487 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))484 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next()))) 488 485 { 489 486 mother = track; 490 track = static_cast<Candidate *>(track->Clone());487 track = static_cast<Candidate *>(track->Clone()); 491 488 track->AddCandidate(mother); 492 489 … … 494 491 } 495 492 } 496 493 497 494 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking 498 else if 499 { 500 weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma *fTrackSigma) : 0.0;501 weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;502 503 bestEnergyEstimate = (weightTrack *fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);504 rescaleFactor = bestEnergyEstimate /fTrackEnergy;505 495 else if(fTrackEnergy > 0.0) 496 { 497 weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma * fTrackSigma) : 0.0; 498 weightCalo = (sigma > 0.0) ? 1 / (sigma * sigma) : 0.0; 499 500 bestEnergyEstimate = (weightTrack * fTrackEnergy + weightCalo * energy) / (weightTrack + weightCalo); 501 rescaleFactor = bestEnergyEstimate / fTrackEnergy; 502 506 503 fItTowerTrackArray->Reset(); 507 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))504 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next()))) 508 505 { 509 506 mother = track; 510 track = static_cast<Candidate *>(track->Clone());507 track = static_cast<Candidate *>(track->Clone()); 511 508 track->AddCandidate(mother); 512 509 … … 516 513 } 517 514 } 518 519 } 515 } 520 516 521 517 //------------------------------------------------------------------------------ … … 527 523 if(mean > 0.0) 528 524 { 529 b = TMath::Sqrt(TMath::Log((1.0 + (sigma *sigma)/(mean*mean))));530 a = TMath::Log(mean) - 0.5 *b*b;531 532 return TMath::Exp(a + b *gRandom->Gaus(0.0, 1.0));525 b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean)))); 526 a = TMath::Log(mean) - 0.5 * b * b; 527 528 return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0)); 533 529 } 534 530 else
Note:
See TracChangeset
for help on using the changeset viewer.