- Timestamp:
- Mar 17, 2021, 4:19:54 PM (3 years ago)
- Branches:
- master
- Children:
- 9cc5aeb
- Parents:
- b8a6aa3
- Location:
- modules
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/AngularSmearing.cc
rb8a6aa3 rfd4b326 98 98 { 99 99 Candidate *candidate, *mother; 100 Double_t pt, eta, phi, e ;100 Double_t pt, eta, phi, e, m; 101 101 102 102 fItInputArray->Reset(); 103 103 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 104 104 { 105 const TLorentzVector &candidatePosition = candidate->Position;106 105 const TLorentzVector &candidateMomentum = candidate->Momentum; 107 eta = candidate Position.Eta();108 phi = candidate Position.Phi();106 eta = candidateMomentum.Eta(); 107 phi = candidateMomentum.Phi(); 109 108 pt = candidateMomentum.Pt(); 110 109 e = candidateMomentum.E(); 110 m = candidateMomentum.M(); 111 111 112 112 // apply smearing formula for eta,phi 113 114 113 eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e, candidate)); 115 114 phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e, candidate)); … … 119 118 mother = candidate; 120 119 candidate = static_cast<Candidate *>(candidate->Clone()); 121 eta = candidateMomentum.Eta(); 122 phi = candidateMomentum.Phi(); 123 candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta)); 120 candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m); 124 121 candidate->AddCandidate(mother); 125 122 -
modules/Calorimeter.cc
rb8a6aa3 rfd4b326 559 559 if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin) 560 560 { 561 // create new photon tower 561 // create new photon tower assuming null mass 562 562 tower = static_cast<Candidate *>(fTower->Clone()); 563 563 pt = ecalNeutralEnergy / TMath::CosH(eta); … … 646 646 track = static_cast<Candidate *>(track->Clone()); 647 647 track->AddCandidate(mother); 648 649 648 track->Momentum *= rescaleFactor; 649 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M()); 650 650 651 651 fEFlowTrackOutputArray->Add(track); -
modules/DenseTrackFilter.cc
rb8a6aa3 rfd4b326 237 237 238 238 Candidate *candidate, *track; 239 Double_t pt, eta, phi ;239 Double_t pt, eta, phi, m; 240 240 Int_t numberOfCandidates; 241 241 … … 251 251 eta = candidate->Momentum.Eta(); 252 252 phi = candidate->Momentum.Phi(); 253 m = candidate->Momentum.M(); 254 253 255 eta = gRandom->Gaus(eta, fEtaPhiRes); 254 256 phi = gRandom->Gaus(phi, fEtaPhiRes); 255 candidate->Momentum.SetPtEtaPhi E(pt, eta, phi, pt * TMath::CosH(eta));257 candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m); 256 258 candidate->AddCandidate(track); 257 259 -
modules/DualReadoutCalorimeter.cc
rb8a6aa3 rfd4b326 24 24 // This implementation of dual calorimetry relies on several approximations: 25 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)) 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. 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. 27 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 28 … … 63 63 fItParticleInputArray(0), fItTrackInputArray(0) 64 64 { 65 65 66 66 fECalResolutionFormula = new DelphesFormula; 67 67 fHCalResolutionFormula = new DelphesFormula; … … 75 75 fTowerTrackArray = new TObjArray; 76 76 fItTowerTrackArray = fTowerTrackArray->MakeIterator(); 77 77 78 78 } 79 79 … … 82 82 DualReadoutCalorimeter::~DualReadoutCalorimeter() 83 83 { 84 84 85 85 if(fECalResolutionFormula) delete fECalResolutionFormula; 86 86 if(fHCalResolutionFormula) delete fHCalResolutionFormula; … … 94 94 if(fTowerTrackArray) delete fTowerTrackArray; 95 95 if(fItTowerTrackArray) delete fItTowerTrackArray; 96 96 97 97 } 98 98 … … 379 379 fHCalTrackEnergy = 0.0; 380 380 fTrackEnergy = 0.0; 381 381 382 382 fECalTrackSigma = 0.0; 383 383 fHCalTrackSigma = 0.0; 384 384 fTrackSigma = 0.0; 385 385 386 386 fTowerTrackHits = 0; 387 387 fTowerPhotonHits = 0; … … 390 390 fHCalTowerTrackArray->Clear(); 391 391 fTowerTrackArray->Clear(); 392 392 393 393 } 394 394 … … 414 414 } 415 415 416 417 /* 418 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 419 { 420 fECalTrackEnergy += ecalEnergy; 421 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 422 if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy; 423 else energyGuess = momentum.E(); 424 425 fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 426 fECalTowerTrackArray->Add(track); 427 } 428 429 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9) 430 { 431 fHCalTrackEnergy += hcalEnergy; 432 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 433 if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy; 434 else energyGuess = momentum.E(); 435 436 fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 437 fHCalTowerTrackArray->Add(track); 438 } 439 440 // muons 441 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 442 { 443 fEFlowTrackOutputArray->Add(track); 444 } 445 */ 446 416 447 417 // in Dual Readout we do not care if tracks are ECAL of HCAL 448 418 if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9) 449 { 419 { 450 420 fTrackEnergy += energy; 451 // 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) 421 // 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) 452 422 sigma = 0.0; 453 423 if(fHCalTrackFractions[number] > 0) … … 455 425 else 456 426 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 457 458 if(sigma/momentum.E() < track->TrackResolution) 459 energyGuess = ecalEnergy + hcalEnergy; 427 428 if(sigma/momentum.E() < track->TrackResolution) 429 energyGuess = ecalEnergy + hcalEnergy; 460 430 else 461 431 energyGuess = momentum.E(); 462 432 463 433 fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 464 434 fTowerTrackArray->Add(track); 465 435 466 436 } 467 437 else … … 511 481 Double_t ecalEnergy, hcalEnergy; 512 482 Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy; 513 483 514 484 Double_t ecalSigma, hcalSigma, sigma; 515 485 Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma; 516 486 517 487 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor; 518 488 519 489 TLorentzVector momentum; 520 490 TFractionMap::iterator itFractionMap; … … 536 506 //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl; 537 507 538 // if no hadronic energy, use ECAL resolution 508 // if no hadronic energy, use ECAL resolution 539 509 if (fHCalTowerEnergy <= fHCalEnergyMin) 540 510 { … … 544 514 } 545 515 546 // if hadronic fraction > 0, use HCAL resolution 516 // if hadronic fraction > 0, use HCAL resolution 547 517 else 548 518 { … … 554 524 energy = LogNormal(energy, sigma); 555 525 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl; 556 526 557 527 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0; 558 528 … … 626 596 627 597 // fill energy flow candidates 628 598 629 599 fTrackSigma = TMath::Sqrt(fTrackSigma); 630 600 neutralEnergy = max( (energy - fTrackEnergy) , 0.0); … … 638 608 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin) 639 609 { 640 610 641 611 //cout<<"significant neutral excess found:"<<endl; 642 612 // create new photon tower … … 646 616 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy); 647 617 648 // if no hadronic energy, use ECAL resolution 618 // if no hadronic energy, use ECAL resolution 649 619 if (fHCalTowerEnergy <= fHCalEnergyMin) 650 620 { … … 654 624 } 655 625 656 // if hadronic fraction > 0, use HCAL resolution 626 // if hadronic fraction > 0, use HCAL resolution 657 627 else 658 628 { … … 663 633 664 634 fEFlowPhotonOutputArray->Add(tower); 665 635 666 636 667 637 //clone tracks … … 669 639 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 670 640 { 671 //cout<<"looping over tracks"<<endl;672 641 mother = track; 673 642 track = static_cast<Candidate*>(track->Clone()); … … 677 646 } 678 647 679 648 680 649 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking 681 650 else if(fTrackEnergy > 0.0) … … 685 654 weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0; 686 655 687 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 656 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 688 657 rescaleFactor = bestEnergyEstimate/fTrackEnergy; 689 658 … … 691 660 fItTowerTrackArray->Reset(); 692 661 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 693 { 662 { 694 663 mother = track; 695 track = static_cast<Candidate *>(track->Clone());664 track = static_cast<Candidate *>(track->Clone()); 696 665 track->AddCandidate(mother); 697 698 track->Momentum *= rescaleFactor; 699 666 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M()); 700 667 fEFlowTrackOutputArray->Add(track); 701 668 } 702 669 } 703 670 704 671 705 672 } -
modules/EnergySmearing.cc
rb8a6aa3 rfd4b326 95 95 { 96 96 Candidate *candidate, *mother; 97 Double_t pt, energy, eta, phi ;97 Double_t pt, energy, eta, phi, m; 98 98 99 99 fItInputArray->Reset(); … … 107 107 phi = candidatePosition.Phi(); 108 108 energy = candidateMomentum.E(); 109 m = candidateMomentum.M(); 109 110 110 111 // apply smearing formula … … 117 118 eta = candidateMomentum.Eta(); 118 119 phi = candidateMomentum.Phi(); 119 candidate->Momentum.SetPtEtaPhiE(energy / TMath::CosH(eta), eta, phi, energy); 120 pt = (energy > m) ? TMath::Sqrt(energy*energy - m*m)/TMath::CosH(eta) : 0; 121 candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 120 122 candidate->TrackResolution = fFormula->Eval(pt, eta, phi, energy) / candidateMomentum.E(); 121 123 candidate->AddCandidate(mother); -
modules/MomentumSmearing.cc
rb8a6aa3 rfd4b326 95 95 { 96 96 Candidate *candidate, *mother; 97 Double_t pt, eta, phi, e, res;97 Double_t pt, eta, phi, e, m, res; 98 98 99 99 fItInputArray->Reset(); … … 106 106 pt = candidateMomentum.Pt(); 107 107 e = candidateMomentum.E(); 108 m = candidateMomentum.M(); 108 109 res = fFormula->Eval(pt, eta, phi, e, candidate); 109 110 … … 121 122 eta = candidateMomentum.Eta(); 122 123 phi = candidateMomentum.Phi(); 123 candidate->Momentum.SetPtEtaPhi E(pt, eta, phi, pt * TMath::CosH(eta));124 candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m); 124 125 //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e); 125 126 candidate->TrackResolution = res; -
modules/SimpleCalorimeter.cc
rb8a6aa3 rfd4b326 507 507 track = static_cast<Candidate *>(track->Clone()); 508 508 track->AddCandidate(mother); 509 510 track->Momentum *= rescaleFactor; 511 509 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M()); 512 510 fEFlowTrackOutputArray->Add(track); 513 511 } -
modules/TrackCovariance.cc
rb8a6aa3 rfd4b326 1 /*1 /* 2 2 * Delphes: a framework for fast simulation of a generic collider experiment 3 3 * Copyright (C) 2020 Universite catholique de Louvain (UCLouvain), Belgium -
modules/TrackSmearing.cc
rb8a6aa3 rfd4b326 158 158 TLorentzVector beamSpotPosition; 159 159 Candidate *candidate, *mother; 160 Double_t pt, eta, e, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;160 Double_t pt, eta, e, m, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi; 161 161 Double_t x, y, z, t, px, py, pz, theta; 162 162 Double_t q, r; … … 224 224 eta = momentum.Eta(); 225 225 e = momentum.E(); 226 226 m = momentum.M(); 227 227 228 d0 = trueD0 = candidate->D0; 228 229 dz = trueDZ = candidate->DZ; … … 332 333 candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta)); 333 334 candidate->Momentum.SetPz(p * TMath::Cos(theta)); 334 candidate->Momentum.SetE( candidate->Momentum.Pt() * TMath::CosH(eta));335 candidate->Momentum.SetE(TMath::Sqrt(p*p + m*m)); 335 336 candidate->PT = candidate->Momentum.Pt(); 336 337 -
modules/TreeWriter.cc
rb8a6aa3 rfd4b326 320 320 Candidate *particle = 0; 321 321 Track *entry = 0; 322 Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi ;322 Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m; 323 323 const Double_t c_light = 2.99792458E8; 324 324 … … 387 387 p = momentum.P(); 388 388 phi = momentum.Phi(); 389 m = momentum.M(); 389 390 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10; 390 391 … … 400 401 entry->CtgTheta = ctgTheta; 401 402 entry->C = candidate->C; 403 entry->Mass = m; 402 404 403 405 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); … … 470 472 Candidate *particle = 0; 471 473 ParticleFlowCandidate *entry = 0; 472 Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi ;474 Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m; 473 475 const Double_t c_light = 2.99792458E8; 474 476 … … 541 543 p = momentum.P(); 542 544 phi = momentum.Phi(); 545 m = momentum.M(); 543 546 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10; 544 547 … … 550 553 entry->CtgTheta = ctgTheta; 551 554 entry->C = candidate->C; 555 entry->Mass = m; 552 556 553 557 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
Note:
See TracChangeset
for help on using the changeset viewer.