Fork me on GitHub

Changeset fd4b326 in git


Ignore:
Timestamp:
Mar 17, 2021, 4:19:54 PM (4 years ago)
Author:
michele <michele.selvaggi@…>
Branches:
master
Children:
9cc5aeb
Parents:
b8a6aa3
Message:

tracks have now montecarlo mass

Files:
11 edited

Legend:

Unmodified
Added
Removed

  • classes/DelphesClasses.h

    rb8a6aa3 rfd4b326  
    385385
    386386  Float_t NeutralEnergyFraction;  // charged energy fraction
    387   Float_t ChargedEnergyFraction;  // neutral energy fraction 
     387  Float_t ChargedEnergyFraction;  // neutral energy fraction
    388388
    389389  Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
    … …  
    439439  Float_t CtgTheta; // track cotangent of theta
    440440  Float_t C; // track curvature inverse
     441  Float_t Mass; // particle mass
    441442
    442443  Float_t EtaOuter; // track pseudorapidity at the tracker edge
    … …  
    472473
    473474  // track covariance off-diagonal terms
    474   Float_t ErrorD0Phi;     
    475   Float_t ErrorD0C;       
    476   Float_t ErrorD0DZ;     
    477   Float_t ErrorD0CtgTheta;   
    478   Float_t ErrorPhiC;     
    479   Float_t ErrorPhiDZ;     
    480   Float_t ErrorPhiCtgTheta ; 
    481   Float_t ErrorCDZ;       
    482   Float_t ErrorCCtgTheta;   
    483   Float_t ErrorDZCtgTheta;   
     475  Float_t ErrorD0Phi;
     476  Float_t ErrorD0C;
     477  Float_t ErrorD0DZ;
     478  Float_t ErrorD0CtgTheta;
     479  Float_t ErrorPhiC;
     480  Float_t ErrorPhiDZ;
     481  Float_t ErrorPhiCtgTheta ;
     482  Float_t ErrorCDZ;
     483  Float_t ErrorCCtgTheta;
     484  Float_t ErrorDZCtgTheta;
    484485
    485486  TRef Particle; // reference to generated particle
    … …  
    542543  Float_t CtgTheta; // track cotangent of theta
    543544  Float_t C; // track curvature inverse
     545  Float_t Mass; // particle mass
    544546
    545547  Float_t EtaOuter; // track pseudorapidity at the tracker edge
    … …  
    575577
    576578  // track covariance off-diagonal terms
    577   Float_t ErrorD0Phi;     
    578   Float_t ErrorD0C;       
    579   Float_t ErrorD0DZ;     
    580   Float_t ErrorD0CtgTheta;   
    581   Float_t ErrorPhiC;     
    582   Float_t ErrorPhiDZ;     
    583   Float_t ErrorPhiCtgTheta ; 
    584   Float_t ErrorCDZ;       
    585   Float_t ErrorCCtgTheta;   
    586   Float_t ErrorDZCtgTheta;   
     579  Float_t ErrorD0Phi;
     580  Float_t ErrorD0C;
     581  Float_t ErrorD0DZ;
     582  Float_t ErrorD0CtgTheta;
     583  Float_t ErrorPhiC;
     584  Float_t ErrorPhiDZ;
     585  Float_t ErrorPhiCtgTheta ;
     586  Float_t ErrorCDZ;
     587  Float_t ErrorCCtgTheta;
     588  Float_t ErrorDZCtgTheta;
    587589
    588590  Int_t VertexIndex; // reference to vertex
    … …  
    710712  Float_t FracPt[5];
    711713  Float_t NeutralEnergyFraction;  // charged energy fraction
    712   Float_t ChargedEnergyFraction;  // neutral energy fraction 
     714  Float_t ChargedEnergyFraction;  // neutral energy fraction
    713715
    714716
    … …  
    767769  // event characteristics variables
    768770  Double_t ParticleDensity; // particle multiplicity density in the proximity of the particle
    769  
     771
    770772  static CompBase *fgCompare; //!
    771773  const CompBase *GetCompare() const { return fgCompare; }

  • modules/AngularSmearing.cc

    rb8a6aa3 rfd4b326  
    9898{
    9999  Candidate *candidate, *mother;
    100   Double_t pt, eta, phi, e;
     100  Double_t pt, eta, phi, e, m;
    101101
    102102  fItInputArray->Reset();
    103103  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    104104  {
    105     const TLorentzVector &candidatePosition = candidate->Position;
    106105    const TLorentzVector &candidateMomentum = candidate->Momentum;
    107     eta = candidatePosition.Eta();
    108     phi = candidatePosition.Phi();
     106    eta = candidateMomentum.Eta();
     107    phi = candidateMomentum.Phi();
    109108    pt = candidateMomentum.Pt();
    110109    e = candidateMomentum.E();
     110    m = candidateMomentum.M();
    111111
    112112    // apply smearing formula for eta,phi
    113 
    114113    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e, candidate));
    115114    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e, candidate));
    … …  
    119118    mother = candidate;
    120119    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);
    124121    candidate->AddCandidate(mother);
    125122

  • modules/Calorimeter.cc

    rb8a6aa3 rfd4b326  
    559559  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    560560  {
    561     // create new photon tower
     561    // create new photon tower assuming null mass
    562562    tower = static_cast<Candidate *>(fTower->Clone());
    563563    pt = ecalNeutralEnergy / TMath::CosH(eta);
    … …  
    646646      track = static_cast<Candidate *>(track->Clone());
    647647      track->AddCandidate(mother);
    648 
    649648      track->Momentum *= rescaleFactor;
     649      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
    650650
    651651      fEFlowTrackOutputArray->Add(track);

  • modules/DenseTrackFilter.cc

    rb8a6aa3 rfd4b326  
    237237
    238238  Candidate *candidate, *track;
    239   Double_t pt, eta, phi;
     239  Double_t pt, eta, phi, m;
    240240  Int_t numberOfCandidates;
    241241
    … …  
    251251  eta = candidate->Momentum.Eta();
    252252  phi = candidate->Momentum.Phi();
     253  m = candidate->Momentum.M();
     254
    253255  eta = gRandom->Gaus(eta, fEtaPhiRes);
    254256  phi = gRandom->Gaus(phi, fEtaPhiRes);
    255   candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
     257  candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
    256258  candidate->AddCandidate(track);
    257259

  • modules/DualReadoutCalorimeter.cc

    rb8a6aa3 rfd4b326  
    2424  // This implementation of dual calorimetry relies on several approximations:
    2525  // - 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.
    2727  //   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.
    2828
    … …  
    6363  fItParticleInputArray(0), fItTrackInputArray(0)
    6464{
    65  
     65
    6666  fECalResolutionFormula = new DelphesFormula;
    6767  fHCalResolutionFormula = new DelphesFormula;
    … …  
    7575  fTowerTrackArray = new TObjArray;
    7676  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    77  
     77
    7878}
    7979
    … …  
    8282DualReadoutCalorimeter::~DualReadoutCalorimeter()
    8383{
    84  
     84
    8585  if(fECalResolutionFormula) delete fECalResolutionFormula;
    8686  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    … …  
    9494  if(fTowerTrackArray) delete fTowerTrackArray;
    9595  if(fItTowerTrackArray) delete fItTowerTrackArray;
    96  
     96
    9797}
    9898
    … …  
    379379      fHCalTrackEnergy = 0.0;
    380380      fTrackEnergy = 0.0;
    381      
     381
    382382      fECalTrackSigma = 0.0;
    383383      fHCalTrackSigma = 0.0;
    384384      fTrackSigma = 0.0;
    385      
     385
    386386      fTowerTrackHits = 0;
    387387      fTowerPhotonHits = 0;
    … …  
    390390      fHCalTowerTrackArray->Clear();
    391391      fTowerTrackArray->Clear();
    392    
     392
    393393    }
    394394
    … …  
    414414      }
    415415
    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
    447417      // in Dual Readout we do not care if tracks are ECAL of HCAL
    448418      if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
    449       { 
     419      {
    450420        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)
    452422        sigma = 0.0;
    453423        if(fHCalTrackFractions[number] > 0)
    … …  
    455425        else
    456426          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;
    460430        else
    461431          energyGuess = momentum.E();
    462              
     432
    463433        fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
    464434        fTowerTrackArray->Add(track);
    465      
     435
    466436      }
    467437      else
    … …  
    511481  Double_t ecalEnergy, hcalEnergy;
    512482  Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
    513  
     483
    514484  Double_t ecalSigma, hcalSigma, sigma;
    515485  Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
    516486
    517487  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    518  
     488
    519489  TLorentzVector momentum;
    520490  TFractionMap::iterator itFractionMap;
    … …  
    536506  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl;
    537507
    538   // if no hadronic energy, use ECAL resolution 
     508  // if no hadronic energy, use ECAL resolution
    539509  if (fHCalTowerEnergy <= fHCalEnergyMin)
    540510  {
    … …  
    544514  }
    545515
    546   // if hadronic fraction > 0, use HCAL resolution 
     516  // if hadronic fraction > 0, use HCAL resolution
    547517  else
    548518  {
    … …  
    554524  energy = LogNormal(energy, sigma);
    555525  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
    556  
     526
    557527  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    558528
    … …  
    626596
    627597  // fill energy flow candidates
    628  
     598
    629599  fTrackSigma = TMath::Sqrt(fTrackSigma);
    630600  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
    … …  
    638608  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
    639609  {
    640    
     610
    641611    //cout<<"significant neutral excess found:"<<endl;
    642612    // create new photon tower
    … …  
    646616    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
    647617
    648     // if no hadronic energy, use ECAL resolution 
     618    // if no hadronic energy, use ECAL resolution
    649619    if (fHCalTowerEnergy <= fHCalEnergyMin)
    650620    {
    … …  
    654624    }
    655625
    656     // if hadronic fraction > 0, use HCAL resolution 
     626    // if hadronic fraction > 0, use HCAL resolution
    657627    else
    658628    {
    … …  
    663633
    664634    fEFlowPhotonOutputArray->Add(tower);
    665    
     635
    666636
    667637    //clone tracks
    … …  
    669639    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    670640    {
    671       //cout<<"looping over tracks"<<endl;
    672641      mother = track;
    673642      track = static_cast<Candidate*>(track->Clone());
    … …  
    677646  }
    678647
    679  
     648
    680649  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
    681650  else if(fTrackEnergy > 0.0)
    … …  
    685654    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
    686655
    687     bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
     656    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
    688657    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
    689658
    … …  
    691660    fItTowerTrackArray->Reset();
    692661    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    693     { 
     662    {
    694663      mother = track;
    695       track = static_cast<Candidate*>(track->Clone());
     664      track = static_cast<Candidate *>(track->Clone());
    696665      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());
    700667      fEFlowTrackOutputArray->Add(track);
    701668    }
    702669  }
    703  
     670
    704671
    705672}

  • modules/EnergySmearing.cc

    rb8a6aa3 rfd4b326  
    9595{
    9696  Candidate *candidate, *mother;
    97   Double_t pt, energy, eta, phi;
     97  Double_t pt, energy, eta, phi, m;
    9898
    9999  fItInputArray->Reset();
    … …  
    107107    phi = candidatePosition.Phi();
    108108    energy = candidateMomentum.E();
     109    m = candidateMomentum.M();
    109110
    110111    // apply smearing formula
    … …  
    117118    eta = candidateMomentum.Eta();
    118119    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);
    120122    candidate->TrackResolution = fFormula->Eval(pt, eta, phi, energy) / candidateMomentum.E();
    121123    candidate->AddCandidate(mother);

  • modules/MomentumSmearing.cc

    rb8a6aa3 rfd4b326  
    9595{
    9696  Candidate *candidate, *mother;
    97   Double_t pt, eta, phi, e, res;
     97  Double_t pt, eta, phi, e, m, res;
    9898
    9999  fItInputArray->Reset();
    … …  
    106106    pt = candidateMomentum.Pt();
    107107    e = candidateMomentum.E();
     108    m = candidateMomentum.M();
    108109    res = fFormula->Eval(pt, eta, phi, e, candidate);
    109110
    … …  
    121122    eta = candidateMomentum.Eta();
    122123    phi = candidateMomentum.Phi();
    123     candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
     124    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
    124125    //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
    125126    candidate->TrackResolution = res;

  • modules/SimpleCalorimeter.cc

    rb8a6aa3 rfd4b326  
    507507      track = static_cast<Candidate *>(track->Clone());
    508508      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());
    512510      fEFlowTrackOutputArray->Add(track);
    513511    }

  • modules/TrackCovariance.cc

    rb8a6aa3 rfd4b326  
    1 /*
     1  /*
    22 *  Delphes: a framework for fast simulation of a generic collider experiment
    33 *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), Belgium

  • modules/TrackSmearing.cc

    rb8a6aa3 rfd4b326  
    158158  TLorentzVector beamSpotPosition;
    159159  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;
    161161  Double_t x, y, z, t, px, py, pz, theta;
    162162  Double_t q, r;
    … …  
    224224    eta = momentum.Eta();
    225225    e = momentum.E();
    226    
     226    m = momentum.M();
     227
    227228    d0 = trueD0 = candidate->D0;
    228229    dz = trueDZ = candidate->DZ;
    … …  
    332333    candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta));
    333334    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));
    335336    candidate->PT = candidate->Momentum.Pt();
    336337

  • modules/TreeWriter.cc

    rb8a6aa3 rfd4b326  
    320320  Candidate *particle = 0;
    321321  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;
    323323  const Double_t c_light = 2.99792458E8;
    324324
    … …  
    387387    p = momentum.P();
    388388    phi = momentum.Phi();
     389    m = momentum.M();
    389390    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
    390391
    … …  
    400401    entry->CtgTheta = ctgTheta;
    401402    entry->C = candidate->C;
     403    entry->Mass = m;
    402404
    403405    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    … …  
    470472  Candidate *particle = 0;
    471473  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;
    473475  const Double_t c_light = 2.99792458E8;
    474476
    … …  
    541543    p = momentum.P();
    542544    phi = momentum.Phi();
     545    m = momentum.M();
    543546    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
    544547
    … …  
    550553    entry->CtgTheta = ctgTheta;
    551554    entry->C = candidate->C;
     555    entry->Mass = m;
    552556
    553557    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
Note: See TracChangeset for help on using the changeset viewer.