Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/Calorimeter.cc

    r298734e r341014c  
    1717 */
    1818
    19 
    2019/** \class Calorimeter
    2120 *
    … …  
    3332#include "classes/DelphesFormula.h"
    3433
     34#include "ExRootAnalysis/ExRootClassifier.h"
     35#include "ExRootAnalysis/ExRootFilter.h"
    3536#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"
    3941#include "TMath.h"
     42#include "TObjArray.h"
     43#include "TRandom3.h"
    4044#include "TString.h"
    41 #include "TFormula.h"
    42 #include "TRandom3.h"
    43 #include "TObjArray.h"
    44 #include "TDatabasePDG.h"
    45 #include "TLorentzVector.h"
    4645
    4746#include <algorithm>
    48 #include <stdexcept>
    4947#include <iostream>
    5048#include <sstream>
     49#include <stdexcept>
    5150
    5251using namespace std;
    … …  
    5857  fItParticleInputArray(0), fItTrackInputArray(0)
    5958{
    60  
     59
    6160  fECalResolutionFormula = new DelphesFormula;
    6261  fHCalResolutionFormula = new DelphesFormula;
    … …  
    6766  fHCalTowerTrackArray = new TObjArray;
    6867  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
    69  
    7068}
    7169
    … …  
    7472Calorimeter::~Calorimeter()
    7573{
    76  
     74
    7775  if(fECalResolutionFormula) delete fECalResolutionFormula;
    7876  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    … …  
    8381  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
    8482  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
    85  
    8683}
    8784
    … …  
    9491  Double_t ecalFraction, hcalFraction;
    9592  TBinMap::iterator itEtaBin;
    96   set< Double_t >::iterator itPhiBin;
    97   vector< Double_t > *phiBins;
     93  set<Double_t>::iterator itPhiBin;
     94  vector<Double_t> *phiBins;
    9895
    9996  // read eta and phi bins
    … …  
    103100  fEtaBins.clear();
    104101  fPhiBins.clear();
    105   for(i = 0; i < size/2; ++i)
    106   {
    107     paramEtaBins = param[i*2];
     102  for(i = 0; i < size / 2; ++i)
     103  {
     104    paramEtaBins = param[i * 2];
    108105    sizeEtaBins = paramEtaBins.GetSize();
    109     paramPhiBins = param[i*2 + 1];
     106    paramPhiBins = param[i * 2 + 1];
    110107    sizePhiBins = paramPhiBins.GetSize();
    111108
    … …  
    124121  {
    125122    fEtaBins.push_back(itEtaBin->first);
    126     phiBins = new vector< double >(itEtaBin->second.size());
     123    phiBins = new vector<double>(itEtaBin->second.size());
    127124    fPhiBins.push_back(phiBins);
    128125    phiBins->clear();
    … …  
    141138  fFractionMap[0] = make_pair(0.0, 1.0);
    142139
    143   for(i = 0; i < size/2; ++i)
    144   {
    145     paramFractions = param[i*2 + 1];
     140  for(i = 0; i < size / 2; ++i)
     141  {
     142    paramFractions = param[i * 2 + 1];
    146143
    147144    ecalFraction = paramFractions[0].GetDouble();
    148145    hcalFraction = paramFractions[1].GetDouble();
    149146
    150     fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
     147    fFractionMap[param[i * 2].GetInt()] = make_pair(ecalFraction, hcalFraction);
    151148  }
    152149
    153150  // read min E value for timing measurement in ECAL
    154   fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
     151  fTimingEnergyMin = GetDouble("TimingEnergyMin", 4.);
    155152  // For timing
    156153  // So far this flag needs to be false
    … …  
    192189void Calorimeter::Finish()
    193190{
    194   vector< vector< Double_t >* >::iterator itPhiBin;
     191  vector<vector<Double_t> *>::iterator itPhiBin;
    195192  if(fItParticleInputArray) delete fItParticleInputArray;
    196193  if(fItTrackInputArray) delete fItTrackInputArray;
    … …  
    218215  TFractionMap::iterator itFractionMap;
    219216
    220   vector< Double_t >::iterator itEtaBin;
    221   vector< Double_t >::iterator itPhiBin;
    222   vector< Double_t > *phiBins;
    223 
    224   vector< Long64_t >::iterator itTowerHits;
     217  vector<Double_t>::iterator itEtaBin;
     218  vector<Double_t>::iterator itPhiBin;
     219  vector<Double_t> *phiBins;
     220
     221  vector<Long64_t>::iterator itTowerHits;
    225222
    226223  DelphesFactory *factory = GetFactory();
    … …  
    234231  fItParticleInputArray->Reset();
    235232  number = -1;
    236   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
     233  while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
    237234  {
    238235    const TLorentzVector &particlePosition = particle->Position;
    … …  
    280277  fItTrackInputArray->Reset();
    281278  number = -1;
    282   while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
     279  while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
    283280  {
    284281    const TLorentzVector &trackPosition = track->Position;
    … …  
    331328    towerHit = (*itTowerHits);
    332329    flags = (towerHit >> 24) & 0x00000000000000FFLL;
    333     number = (towerHit) & 0x0000000000FFFFFFLL;
     330    number = (towerHit)&0x0000000000FFFFFFLL;
    334331    hitEtaPhi = towerHit >> 32;
    335332
    … …  
    352349
    353350      // calculate eta and phi of the tower's center
    354       fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
    355       fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
     351      fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
     352      fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
    356353
    357354      fTowerEdges[0] = fEtaBins[etaBin - 1];
    … …  
    365362      fECalTrackEnergy = 0.0;
    366363      fHCalTrackEnergy = 0.0;
    367      
     364
    368365      fECalTrackSigma = 0.0;
    369366      fHCalTrackSigma = 0.0;
    370      
     367
    371368      fTowerTrackHits = 0;
    372369      fTowerPhotonHits = 0;
    … …  
    374371      fECalTowerTrackArray->Clear();
    375372      fHCalTowerTrackArray->Clear();
    376    
    377373    }
    378374
    … …  
    382378      ++fTowerTrackHits;
    383379
    384       track = static_cast<Candidate*>(fTrackInputArray->At(number));
     380      track = static_cast<Candidate *>(fTrackInputArray->At(number));
    385381      momentum = track->Momentum;
    386382      position = track->Position;
    … …  
    400396      {
    401397        fECalTrackEnergy += ecalEnergy;
    402         ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());       
    403         if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;       
    404         else energyGuess = momentum.E();
    405 
    406         fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
     398        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
     399        if(ecalSigma / momentum.E() < track->TrackResolution)
     400          energyGuess = ecalEnergy;
     401        else
     402          energyGuess = momentum.E();
     403
     404        fECalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
    407405        fECalTowerTrackArray->Add(track);
    408406      }
    409      
     407
    410408      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    411409      {
    412410        fHCalTrackEnergy += hcalEnergy;
    413411        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    414         if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
    415         else energyGuess = momentum.E();
    416 
    417         fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
     412        if(hcalSigma / momentum.E() < track->TrackResolution)
     413          energyGuess = hcalEnergy;
     414        else
     415          energyGuess = momentum.E();
     416
     417        fHCalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
    418418        fHCalTowerTrackArray->Add(track);
    419419      }
    420      
     420
    421421      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    422422      {
    … …  
    430430    if(flags & 2) ++fTowerPhotonHits;
    431431
    432     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     432    particle = static_cast<Candidate *>(fParticleInputArray->At(number));
    433433    momentum = particle->Momentum;
    434434    position = particle->Position;
    … …  
    443443    if(ecalEnergy > fTimingEnergyMin && fTower)
    444444    {
    445       if (abs(particle->PID) != 11 || !fElectronsFromTrack)
     445      if(abs(particle->PID) != 11 || !fElectronsFromTrack)
    446446      {
    447447        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
    … …  
    464464  Double_t ecalEnergy, hcalEnergy;
    465465  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
    466  
     466
    467467  Double_t ecalSigma, hcalSigma;
    468468  Double_t ecalNeutralSigma, hcalNeutralSigma;
    469469
    470470  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    471  
     471
    472472  TLorentzVector momentum;
    473473  TFractionMap::iterator itFractionMap;
    … …  
    486486  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    487487
    488   if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
    489   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
     488  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin * ecalSigma) ecalEnergy = 0.0;
     489  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin * hcalSigma) hcalEnergy = 0.0;
    490490
    491491  energy = ecalEnergy + hcalEnergy;
    … …  
    519519  if(sumWeight > 0.0)
    520520  {
    521     fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
     521    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime / sumWeight);
    522522  }
    523523  else
    … …  
    525525    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
    526526  }
    527 
    528527
    529528  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    … …  
    545544    fTowerOutputArray->Add(fTower);
    546545  }
    547  
     546
    548547  // fill energy flow candidates
    549548  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
    … …  
    551550
    552551  //compute neutral excesses
    553   ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
    554   hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
    555  
    556   ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
    557   hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
    558  
    559    // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
     552  ecalNeutralEnergy = max((ecalEnergy - fECalTrackEnergy), 0.0);
     553  hcalNeutralEnergy = max((hcalEnergy - fHCalTrackEnergy), 0.0);
     554
     555  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma * fECalTrackSigma + ecalSigma * ecalSigma);
     556  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma * fHCalTrackSigma + hcalSigma * hcalSigma);
     557
     558  // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
    560559  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    561560  {
    562561    // create new photon tower
    563     tower = static_cast<Candidate*>(fTower->Clone());
    564     pt =  ecalNeutralEnergy / TMath::CosH(eta);
    565    
     562    tower = static_cast<Candidate *>(fTower->Clone());
     563    pt = ecalNeutralEnergy / TMath::CosH(eta);
     564
    566565    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
    567566    tower->Eem = ecalNeutralEnergy;
    568567    tower->Ehad = 0.0;
    569568    tower->PID = 22;
    570    
     569
    571570    fEFlowPhotonOutputArray->Add(tower);
    572    
     571
    573572    //clone tracks
    574573    fItECalTowerTrackArray->Reset();
    575     while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     574    while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
    576575    {
    577576      mother = track;
    578       track = static_cast<Candidate*>(track->Clone());
     577      track = static_cast<Candidate *>(track->Clone());
    579578      track->AddCandidate(mother);
    580579
    581580      fEFlowTrackOutputArray->Add(track);
    582581    }
    583  
    584   }
    585  
     582  }
     583
    586584  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    587585  else if(fECalTrackEnergy > 0.0)
    588586  {
    589     weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
    590     weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
    591  
    592     bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
    593     rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     587    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma * fECalTrackSigma) : 0.0;
     588    weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma * ecalSigma) : 0.0;
     589
     590    bestEnergyEstimate = (weightTrack * fECalTrackEnergy + weightCalo * ecalEnergy) / (weightTrack + weightCalo);
     591    rescaleFactor = bestEnergyEstimate / fECalTrackEnergy;
    594592
    595593    //rescale tracks
    596594    fItECalTowerTrackArray->Reset();
    597     while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
    598     { 
     595    while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
     596    {
    599597      mother = track;
    600       track = static_cast<Candidate*>(track->Clone());
     598      track = static_cast<Candidate *>(track->Clone());
    601599      track->AddCandidate(mother);
    602600
    … …  
    606604    }
    607605  }
    608 
    609606
    610607  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
    … …  
    612609  {
    613610    // create new photon tower
    614     tower = static_cast<Candidate*>(fTower->Clone());
    615     pt =  hcalNeutralEnergy / TMath::CosH(eta);
    616    
     611    tower = static_cast<Candidate *>(fTower->Clone());
     612    pt = hcalNeutralEnergy / TMath::CosH(eta);
     613
    617614    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
    618615    tower->Ehad = hcalNeutralEnergy;
    619616    tower->Eem = 0.0;
    620    
     617
    621618    fEFlowNeutralHadronOutputArray->Add(tower);
    622    
     619
    623620    //clone tracks
    624621    fItHCalTowerTrackArray->Reset();
    625     while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     622    while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
    626623    {
    627624      mother = track;
    628       track = static_cast<Candidate*>(track->Clone());
     625      track = static_cast<Candidate *>(track->Clone());
    629626      track->AddCandidate(mother);
    630627
    631628      fEFlowTrackOutputArray->Add(track);
    632629    }
    633  
    634   }
    635  
     630  }
     631
    636632  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    637633  else if(fHCalTrackEnergy > 0.0)
    638634  {
    639     weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
    640     weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
    641  
    642     bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     635    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma * fHCalTrackSigma) : 0.0;
     636    weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma * hcalSigma) : 0.0;
     637
     638    bestEnergyEstimate = (weightTrack * fHCalTrackEnergy + weightCalo * hcalEnergy) / (weightTrack + weightCalo);
    643639    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
    644640
    645641    //rescale tracks
    646642    fItHCalTowerTrackArray->Reset();
    647     while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
    648     { 
     643    while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
     644    {
    649645      mother = track;
    650       track = static_cast<Candidate*>(track->Clone());
     646      track = static_cast<Candidate *>(track->Clone());
    651647      track->AddCandidate(mother);
    652648
    … …  
    656652    }
    657653  }
    658  
    659  
    660654}
    661655
    … …  
    668662  if(mean > 0.0)
    669663  {
    670     b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
    671     a = TMath::Log(mean) - 0.5*b*b;
    672 
    673     return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
     664    b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
     665    a = TMath::Log(mean) - 0.5 * b * b;
     666
     667    return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0));
    674668  }
    675669  else
Note: See TracChangeset for help on using the changeset viewer.