Fork me on GitHub

Changeset 3db5282 in git


Ignore:
Timestamp:
Jan 27, 2015, 5:32:59 PM (10 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
2862770
Parents:
2d494a6
Message:

included timing information in Calorimeter

Files:
6 edited

Legend:

Unmodified
Added
Removed
  • classes/DelphesClasses.cc

    r2d494a6 r3db5282  
    133133  MeanSqDeltaR(0),
    134134  PTD(0),
     135  Ntimes(-1),
    135136  IsolationVar(-999),
    136137  IsolationVarRhoCorr(-999),
     
    255256  object.MeanSqDeltaR = MeanSqDeltaR;
    256257  object.PTD = PTD;
    257  
     258  object.Ntimes = Ntimes;
    258259  object.IsolationVar = IsolationVar;
    259260  object.IsolationVarRhoCorr = IsolationVarRhoCorr;
     
    276277  object.fFactory = fFactory;
    277278  object.fArray = 0;
     279
     280
     281  // Copy cluster timing info
     282  copy(Ecal_E_t.begin(),Ecal_E_t.end(),back_inserter(object.Ecal_E_t));
    278283
    279284  if(fArray && fArray->GetEntriesFast() > 0)
     
    326331  PTD = 0.0;
    327332 
     333  Ntimes = 0;
     334  Ecal_E_t.clear();
     335 
    328336  IsolationVar = -999;
    329337  IsolationVarRhoCorr = -999;
  • classes/DelphesClasses.h

    r2d494a6 r3db5282  
    419419  Float_t E; // calorimeter tower energy
    420420
    421   Float_t T; //particle arrival time of flight
    422 
     421  Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
     422  Int_t   Ntimes; // number of hits contributing to time measurement
     423 
    423424  Float_t Eem; // calorimeter tower electromagnetic energy
    424425  Float_t Ehad; // calorimeter tower hadronic energy
     
    510511  Float_t  PTD;
    511512  Float_t  FracPt[5];
     513 
     514  //Timing information
     515 
     516  Int_t    Ntimes;
     517  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t;
    512518
    513519  // Isolation variables
     
    523529
    524530  Float_t Tau[5];
     531 
     532 
     533 
    525534
    526535  static CompBase *fgCompare; //!
  • modules/Calorimeter.cc

    r2d494a6 r3db5282  
    150150*/
    151151
     152  // read min E value for timing measurement in ECAL
     153  fTimingEMin = GetDouble("TimingEMin",4.);
     154  // For timing
     155  // So far this flag needs to be false
     156  // Curved extrapolation not supported
     157  fElectronsFromTrack = false;
     158
     159 
    152160  // read min E value for towers to be saved
    153161  fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
     
    356364      fTrackHCalEnergy = 0.0;
    357365
    358       fTowerECalTime = 0.0;
    359       fTowerHCalTime = 0.0;
    360 
    361       fTrackECalTime = 0.0;
    362       fTrackHCalTime = 0.0;
    363 
    364       fTowerECalTimeWeight = 0.0;
    365       fTowerHCalTimeWeight = 0.0;
    366 
    367366      fTowerTrackHits = 0;
    368367      fTowerPhotonHits = 0;
     
    380379      position = track->Position;
    381380
    382 
    383381      ecalEnergy = momentum.E() * fTrackECalFractions[number];
    384382      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
     
    387385      fTrackHCalEnergy += hcalEnergy;
    388386
    389       fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
    390       fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
    391 
    392       fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
    393       fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
     387      bool dbg_scz = false;
     388      if (dbg_scz) {
     389        cout << "   Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T()
     390             << endl;
     391        Candidate *prt = static_cast<Candidate*>(track->GetCandidates()->Last());
     392        const TLorentzVector &ini = prt->Position;
     393
     394        cout << "                and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T();
     395
     396      }
     397     
     398      if (ecalEnergy > fTimingEMin && fTower) {
     399        if (fElectronsFromTrack) {
     400          //      cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
     401          //           << " PID=" << track->PID << endl;
     402          fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T()));
     403        } else {
     404          //      cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
     405          //           << " PID=" << track->PID << endl;
     406        }
     407      }
     408
    394409
    395410      fTowerTrackArray->Add(track);
     
    397412      continue;
    398413    }
    399 
     414 
    400415    // check for photon and electron hits in current tower
    401416    if(flags & 2) ++fTowerPhotonHits;
     
    412427    fTowerHCalEnergy += hcalEnergy;
    413428
    414     fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
    415     fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
    416 
    417     fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
    418     fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
    419 
     429    if (ecalEnergy > fTimingEMin && fTower) {
     430      if (abs(particle->PID) != 11 || !fElectronsFromTrack) {
     431        //      cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
     432        //   << " PID=" << particle->PID << endl;
     433        fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T()));
     434      } else {
     435       
     436        // N.B. Only charged particle set to leave ecal energy is the electrons
     437        //      cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
     438        //           << " PID=" << particle->PID << endl;
     439       
     440      }
     441    }
    420442
    421443    fTower->AddCandidate(particle);
     
    434456  Double_t ecalEnergy, hcalEnergy;
    435457  Double_t ecalSigma, hcalSigma;
    436   Double_t ecalTime, hcalTime, time;
    437 
     458 
    438459  if(!fTower) return;
    439460
     
    444465  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
    445466
    446   ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
    447   hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
    448 
    449467  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    450468  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
     
    454472
    455473  energy = ecalEnergy + hcalEnergy;
    456   time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
    457 
     474   
    458475  if(fSmearTowerCenter)
    459476  {
     
    469486  pt = energy / TMath::CosH(eta);
    470487
    471   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
     488  // Time calculation for tower
     489  fTower->Ntimes = 0;
     490  Float_t tow_sumT = 0;
     491  Float_t tow_sumW = 0;
     492 
     493  for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++)
     494  {
     495    Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first);
     496    tow_sumT += w*fTower->Ecal_E_t[i].second;
     497    tow_sumW += w;
     498    fTower->Ntimes++;
     499  }
     500 
     501  if (tow_sumW > 0.) {
     502    fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW);
     503  } else {
     504    fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.);
     505  }
     506
     507
    472508  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    473509  fTower->Eem = ecalEnergy;
  • modules/Calorimeter.h

    r2d494a6 r3db5282  
    6060  Double_t fTrackECalEnergy, fTrackHCalEnergy;
    6161
    62   Double_t fTowerECalTime, fTowerHCalTime;
    63   Double_t fTrackECalTime, fTrackHCalTime;
    64 
    65   Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
    66   Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
    67 
     62  Double_t fTimingEMin;
     63  Bool_t fElectronsFromTrack; // for timing
     64 
    6865  Int_t fTowerTrackHits, fTowerPhotonHits;
    6966
  • modules/Isolation.cc

    r2d494a6 r3db5282  
    184184    // loop over all input tracks
    185185   
     186    sumNeutral = 0.0;
    186187    sumCharged = 0.0;
    187188    sumChargedPU = 0.0;
  • modules/TreeWriter.cc

    r2d494a6 r3db5282  
    362362
    363363    entry->T = position.T()*1.0E-3/c_light;
     364    entry->Ntimes = candidate->Ntimes;   
    364365
    365366    FillParticles(candidate, &entry->Particles);
Note: See TracChangeset for help on using the changeset viewer.