Fork me on GitHub

Changes in / [425cd17:5c402f7] in git


Ignore:
Files:
3 added
13 edited

Legend:

Unmodified
Added
Removed

  • external/PUPPI/PuppiAlgo.cc

    r425cd17 r5c402f7  
    1 //#include "FWCore/ParameterSet/interface/ParameterSet.h"
    2 //#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
    31#include "PuppiAlgo.hh"
    4 #include "external/fastjet/internal/base.hh"
    52#include "Math/QuantFuncMathCore.h"
    63#include "Math/SpecFuncMathCore.h"
    … …  
    1916  std::vector<AlgoSubObj> lAlgos = iAlgo.subAlgos;
    2017  fNAlgos = lAlgos.size();
    21   std::cout << "==> "  << fEtaMin << " - " << fEtaMax << " - " << fPtMin  << " - " << fNeutralPtMin  << " - " << fNeutralPtSlope  << " - " << fRMSEtaSF  << " - " << std::endl;
     18  //std::cout << "==> "  << fEtaMin << " - " << fEtaMax << " - " << fPtMin  << " - " << fNeutralPtMin  << " - " << fNeutralPtSlope  << " - " << fRMSEtaSF  << " - " << std::endl;
    2219  for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++)  {
    2320    int    pAlgoId      = lAlgos[i0].metricId;
    … …  
    2825    double pRMSPtMin    = lAlgos[i0].rmsPtMin;
    2926    double pRMSSF       = lAlgos[i0].rmsScaleFactor;
    30     std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
     27    //std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
    3128    fAlgoId        .push_back(pAlgoId);
    3229    fCharged       .push_back(pCharged);

  • external/PUPPI/PuppiAlgo.hh

    r425cd17 r5c402f7  
    1 #include "external/fastjet/internal/base.hh"
    2 #include "external/fastjet/PseudoJet.hh"
     1#include "fastjet/PseudoJet.hh"
    32#include "AlgoObj.hh"
    43#include <vector>

  • external/PUPPI/PuppiContainer.cc

    r425cd17 r5c402f7  
    11#include "PuppiContainer.hh"
    2 #include "external/fastjet/internal/base.hh"
    3 #include "external/fastjet/Selector.hh"
    4 //#include "external/fastjet/FunctionOfPseudoJet.hh"
     2#include "fastjet/Selector.hh"
    53#include "Math/ProbFunc.h"
    64#include "TMath.h"

  • external/PUPPI/PuppiContainer.hh

    r425cd17 r5c402f7  
    11#include "PuppiAlgo.hh"
    22#include "RecoObj2.hh"
    3 //#include "external/fastjet/internal/base.hh"
    4 #include "external/fastjet/PseudoJet.hh"
     3#include "fastjet/PseudoJet.hh"
    54#include <vector>
    65

  • external/PUPPI/puppiCleanContainer.hh

    r425cd17 r5c402f7  
    77#include "PUPPI/puppiAlgoBin.hh"
    88
    9 #include "fastjet/internal/base.hh"
    109#include "fastjet/PseudoJet.hh"
    1110#include <algorithm>

  • modules/Calorimeter.cc

    r425cd17 r5c402f7  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fECalResolutionFormula = new DelphesFormula;
    6362  fHCalResolutionFormula = new DelphesFormula;
    6463
    65   for(i = 0; i < 2; ++i)
    66   {
    67     fECalTowerTrackArray[i] = new TObjArray;
    68     fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
    69 
    70     fHCalTowerTrackArray[i] = new TObjArray;
    71     fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
    72   }
     64  fECalTowerTrackArray = new TObjArray;
     65  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
     66
     67  fHCalTowerTrackArray = new TObjArray;
     68  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
     69 
    7370}
    7471
    … …  
    7774Calorimeter::~Calorimeter()
    7875{
    79   Int_t i;
    80 
     76 
    8177  if(fECalResolutionFormula) delete fECalResolutionFormula;
    8278  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    8379
    84   for(i = 0; i < 2; ++i)
    85   {
    86     if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
    87     if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
    88 
    89     if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
    90     if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
    91   }
     80  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
     81  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
     82
     83  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
     84  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
     85 
    9286}
    9387
    … …  
    218212  Double_t ecalFraction, hcalFraction;
    219213  Double_t ecalEnergy, hcalEnergy;
    220   Double_t ecalSigma, hcalSigma;
    221214  Int_t pdgCode;
    222215
    … …  
    368361      fHCalTowerEnergy = 0.0;
    369362
    370       fECalTrackEnergy[0] = 0.0;
    371       fECalTrackEnergy[1] = 0.0;
    372 
    373       fHCalTrackEnergy[0] = 0.0;
    374       fHCalTrackEnergy[1] = 0.0;
    375 
     363      fECalTrackEnergy = 0.0;
     364      fHCalTrackEnergy = 0.0;
     365     
     366      fECalTrackSigma = 0.0;
     367      fHCalTrackSigma = 0.0;
     368     
    376369      fTowerTrackHits = 0;
    377370      fTowerPhotonHits = 0;
    378371
    379       fECalTowerTrackArray[0]->Clear();
    380       fECalTowerTrackArray[1]->Clear();
    381 
    382       fHCalTowerTrackArray[0]->Clear();
    383       fHCalTowerTrackArray[1]->Clear();
     372      fECalTowerTrackArray->Clear();
     373      fHCalTowerTrackArray->Clear();
     374   
    384375    }
    385376
    … …  
    406397      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    407398      {
    408         ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    409         if(ecalSigma/momentum.E() < track->TrackResolution)
    410         {
    411           fECalTrackEnergy[0] += ecalEnergy;
    412           fECalTowerTrackArray[0]->Add(track);
    413         }
    414         else
    415         {
    416           fECalTrackEnergy[1] += ecalEnergy;
    417           fECalTowerTrackArray[1]->Add(track);
    418         }
     399        fECalTrackEnergy += ecalEnergy;
     400        fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     401        fECalTowerTrackArray->Add(track);
    419402      }
     403     
    420404      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    421405      {
    422         hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    423         if(hcalSigma/momentum.E() < track->TrackResolution)
    424         {
    425           fHCalTrackEnergy[0] += hcalEnergy;
    426           fHCalTowerTrackArray[0]->Add(track);
    427         }
    428         else
    429         {
    430           fHCalTrackEnergy[1] += hcalEnergy;
    431           fHCalTowerTrackArray[1]->Add(track);
    432         }
     406        fHCalTrackEnergy += hcalEnergy;
     407        fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     408        fHCalTowerTrackArray->Add(track);
    433409      }
     410     
    434411      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    435412      {
    … …  
    476453  Double_t energy, pt, eta, phi;
    477454  Double_t ecalEnergy, hcalEnergy;
     455  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
     456 
    478457  Double_t ecalSigma, hcalSigma;
    479 
     458  Double_t ecalNeutralSigma, hcalNeutralSigma;
     459
     460  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
     461 
    480462  TLorentzVector momentum;
    481463  TFractionMap::iterator itFractionMap;
    … …  
    484466
    485467  if(!fTower) return;
    486 
    487468
    488469  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    554535    fTowerOutputArray->Add(fTower);
    555536  }
    556 
     537 
    557538  // fill energy flow candidates
    558 
    559   ecalEnergy -= fECalTrackEnergy[1];
    560   hcalEnergy -= fHCalTrackEnergy[1];
    561 
    562   fItECalTowerTrackArray[0]->Reset();
    563   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
    564   {
    565     mother = track;
    566     track = static_cast<Candidate*>(track->Clone());
    567     track->AddCandidate(mother);
    568 
    569     track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
    570 
    571     fEFlowTrackOutputArray->Add(track);
    572   }
    573 
    574   fItECalTowerTrackArray[1]->Reset();
    575   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
    576   {
    577     mother = track;
    578     track = static_cast<Candidate*>(track->Clone());
    579     track->AddCandidate(mother);
    580 
    581     fEFlowTrackOutputArray->Add(track);
    582   }
    583 
    584   fItHCalTowerTrackArray[0]->Reset();
    585   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
    586   {
    587     mother = track;
    588     track = static_cast<Candidate*>(track->Clone());
    589     track->AddCandidate(mother);
    590 
    591     track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
    592 
    593     fEFlowTrackOutputArray->Add(track);
    594   }
    595 
    596   fItHCalTowerTrackArray[1]->Reset();
    597   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
    598   {
    599     mother = track;
    600     track = static_cast<Candidate*>(track->Clone());
    601     track->AddCandidate(mother);
    602 
    603     fEFlowTrackOutputArray->Add(track);
    604   }
    605 
    606   if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
    607   if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0;
    608 
    609   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    610   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    611 
    612   if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
    613   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
    614 
    615   energy = ecalEnergy + hcalEnergy;
    616 
    617   if(ecalEnergy > 0.0)
     539  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
     540  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
     541
     542  //compute neutral excesses
     543  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
     544  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
     545 
     546  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
     547  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
     548 
     549   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
     550  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    618551  {
    619552    // create new photon tower
    620553    tower = static_cast<Candidate*>(fTower->Clone());
    621 
    622     pt = ecalEnergy / TMath::CosH(eta);
    623 
    624     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
    625     tower->Eem = ecalEnergy;
     554    pt =  ecalNeutralEnergy / TMath::CosH(eta);
     555   
     556    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     557    tower->Eem = ecalNeutralEnergy;
    626558    tower->Ehad = 0.0;
    627559    tower->PID = 22;
    628 
     560   
    629561    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron tower
     562   
     563    //clone tracks
     564    fItECalTowerTrackArray->Reset();
     565    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     566    {
     567      mother = track;
     568      track = static_cast<Candidate*>(track->Clone());
     569      track->AddCandidate(mother);
     570
     571      fEFlowTrackOutputArray->Add(track);
     572    }
     573 
     574  }
     575 
     576  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     577  else if(fECalTrackEnergy > 0.0)
     578  {
     579    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
     580    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
     581 
     582    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
     583    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     584
     585    //rescale tracks
     586    fItECalTowerTrackArray->Reset();
     587    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     588    { 
     589      mother = track;
     590      track = static_cast<Candidate*>(track->Clone());
     591      track->AddCandidate(mother);
     592
     593      track->Momentum *= rescaleFactor;
     594
     595      fEFlowTrackOutputArray->Add(track);
     596    }
     597  }
     598
     599
     600  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
     601  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
     602  {
     603    // create new photon tower
    634604    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     605    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     606   
     607    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     608    tower->Ehad = hcalNeutralEnergy;
    639609    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     610   
    642611    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     612   
     613    //clone tracks
     614    fItHCalTowerTrackArray->Reset();
     615    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     616    {
     617      mother = track;
     618      track = static_cast<Candidate*>(track->Clone());
     619      track->AddCandidate(mother);
     620
     621      fEFlowTrackOutputArray->Add(track);
     622    }
     623 
     624  }
     625 
     626  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     627  else if(fHCalTrackEnergy > 0.0)
     628  {
     629    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
     630    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
     631 
     632    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     633    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
     634
     635    //rescale tracks
     636    fItHCalTowerTrackArray->Reset();
     637    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     638    { 
     639      mother = track;
     640      track = static_cast<Candidate*>(track->Clone());
     641      track->AddCandidate(mother);
     642
     643      track->Momentum *= rescaleFactor;
     644
     645      fEFlowTrackOutputArray->Add(track);
     646    }
     647  }
     648 
     649 
    644650}
    645651

  • modules/Calorimeter.h

    r425cd17 r5c402f7  
    5858  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    5959  Double_t fECalTowerEnergy, fHCalTowerEnergy;
    60   Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
     60  Double_t fECalTrackEnergy, fHCalTrackEnergy;
    6161
    6262  Double_t fTimingEnergyMin;
    … …  
    7070  Double_t fECalEnergySignificanceMin;
    7171  Double_t fHCalEnergySignificanceMin;
     72
     73  Double_t fECalTrackSigma;
     74  Double_t fHCalTrackSigma;
    7275
    7376  Bool_t fSmearTowerCenter;
    … …  
    103106  TObjArray *fEFlowNeutralHadronOutputArray; //!
    104107
    105   TObjArray *fECalTowerTrackArray[2]; //!
    106   TIterator *fItECalTowerTrackArray[2]; //!
     108  TObjArray *fECalTowerTrackArray; //!
     109  TIterator *fItECalTowerTrackArray; //!
    107110
    108   TObjArray *fHCalTowerTrackArray[2]; //!
    109   TIterator *fItHCalTowerTrackArray[2]; //!
     111  TObjArray *fHCalTowerTrackArray; //!
     112  TIterator *fItHCalTowerTrackArray; //!
    110113
    111114  void FinalizeTower();

  • modules/RunPUPPI.cc

    r425cd17 r5c402f7  
    33#include "PUPPI/RecoObj2.hh"
    44#include "PUPPI/AlgoObj.hh"
    5 //#include "PUPPI/puppiParticle.hh"
    6 //#include "PUPPI/puppiAlgoBin.hh"
     5#include "PUPPI/PuppiContainer.hh"
     6
     7#include "fastjet/PseudoJet.hh"
    78
    89#include "classes/DelphesClasses.h"
    … …  
    1718
    1819using namespace std;
     20using namespace fastjet;
    1921
    2022//------------------------------------------------------------------------------
    … …  
    238240  fPuppi->initialize(puppiInputVector);
    239241  fPuppi->puppiWeights();
    240   std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
     242  std::vector<PseudoJet> puppiParticles = fPuppi->puppiParticles();
    241243
    242244  // Loop on final particles
    243   for (std::vector<fastjet::PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
     245  for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
    244246    if(it->user_index() <= int(InputParticles.size())){     
    245247      candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());

  • modules/RunPUPPI.h

    r425cd17 r5c402f7  
    33
    44#include "classes/DelphesModule.h"
    5 #include "PUPPI/PuppiContainer.hh"
    65#include <vector>
    76
    87class TObjArray;
    98class TIterator;
    10 
     9class PuppiContainer;
    1110
    1211class RunPUPPI: public DelphesModule {

  • modules/SimpleCalorimeter.cc

    r425cd17 r5c402f7  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fResolutionFormula = new DelphesFormula;
    63 
    64   for(i = 0; i < 2; ++i)
    65   {
    66     fTowerTrackArray[i] = new TObjArray;
    67     fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
    68   }
     62  fTowerTrackArray = new TObjArray;
     63  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
     64 
    6965}
    7066
    … …  
    7369SimpleCalorimeter::~SimpleCalorimeter()
    7470{
    75   Int_t i;
    76 
     71 
    7772  if(fResolutionFormula) delete fResolutionFormula;
    78 
    79   for(i = 0; i < 2; ++i)
    80   {
    81     if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
    82     if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
    83   }
     73  if(fTowerTrackArray) delete fTowerTrackArray;
     74  if(fItTowerTrackArray) delete fItTowerTrackArray;
     75 
    8476}
    8577
    … …  
    198190  Double_t fraction;
    199191  Double_t energy;
    200   Double_t sigma;
    201192  Int_t pdgCode;
    202193
    … …  
    340331      fTowerEnergy = 0.0;
    341332
    342       fTrackEnergy[0] = 0.0;
    343       fTrackEnergy[1] = 0.0;
     333      fTrackEnergy = 0.0;
     334      fTrackSigma = 0.0;
    344335
    345336      fTowerTime = 0.0;
    … …  
    351342      fTowerPhotonHits = 0;
    352343
    353       fTowerTrackArray[0]->Clear();
    354       fTowerTrackArray[1]->Clear();
    355     }
     344      fTowerTrackArray->Clear();
     345     }
    356346
    357347    // check for track hits
    … …  
    366356      energy = momentum.E() * fTrackFractions[number];
    367357
    368       fTrackTime += energy*position.T();
    369       fTrackTimeWeight += energy;
     358      fTrackTime += TMath::Sqrt(energy)*position.T();
     359      fTrackTimeWeight += TMath::Sqrt(energy);
    370360
    371361      if(fTrackFractions[number] > 1.0E-9)
    372362      {
    373         sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    374         if(sigma/momentum.E() < track->TrackResolution)
    375         {
    376           fTrackEnergy[0] += energy;
    377           fTowerTrackArray[0]->Add(track);
    378         }
    379         else
    380         {
    381           fTrackEnergy[1] += energy;
    382           fTowerTrackArray[1]->Add(track);
    383         }
     363             
     364       // compute total charged energy   
     365       fTrackEnergy += energy;
     366       fTrackSigma += ((track->TrackResolution)*momentum.E())*((track->TrackResolution)*momentum.E());
     367       
     368       fTowerTrackArray->Add(track);
     369     
    384370      }
     371       
    385372      else
    386373      {
    … …  
    418405{
    419406  Candidate *tower, *track, *mother;
    420   Double_t energy, pt, eta, phi;
    421   Double_t sigma;
     407  Double_t energy,neutralEnergy, pt, eta, phi;
     408  Double_t sigma, neutralSigma;
    422409  Double_t time;
     410   
     411  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    423412
    424413  TLorentzVector momentum;
    … …  
    436425
    437426  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
     427
    438428
    439429  if(fSmearTowerCenter)
    … …  
    464454  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    465455
    466   // fill e-flow candidates
    467 
    468   energy -= fTrackEnergy[1];
    469 
    470   fItTowerTrackArray[0]->Reset();
    471   while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
    472   {
    473     mother = track;
    474     track = static_cast<Candidate*>(track->Clone());
    475     track->AddCandidate(mother);
    476 
    477     track->Momentum *= energy/fTrackEnergy[0];
    478 
    479     fEFlowTrackOutputArray->Add(track);
    480   }
    481 
    482   fItTowerTrackArray[1]->Reset();
    483   while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
    484   {
    485     mother = track;
    486     track = static_cast<Candidate*>(track->Clone());
    487     track->AddCandidate(mother);
    488 
    489     fEFlowTrackOutputArray->Add(track);
    490   }
    491 
    492   if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
    493 
    494   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
    495   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    496 
    497   // save energy excess as an energy flow tower
    498   if(energy > 0.0)
     456 
     457  // e-flow candidates
     458
     459  //compute neutral excess
     460 
     461  fTrackSigma = TMath::Sqrt(fTrackSigma);
     462  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
     463 
     464  //compute sigma_trk total
     465  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
     466   
     467  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
     468  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
    499469  {
    500470    // create new photon tower
    501471    tower = static_cast<Candidate*>(fTower->Clone());
    502     pt = energy / TMath::CosH(eta);
    503 
    504     tower->Eem = (!fIsEcal) ? 0 : energy;
    505     tower->Ehad = (fIsEcal) ? 0 : energy;
     472    pt = neutralEnergy / TMath::CosH(eta);
     473
     474    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
     475    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     476    tower->PID = (fIsEcal) ? 22 : 0;
     477 
     478    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     479    fEFlowTowerOutputArray->Add(tower);
    506480   
    507     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     481    fItTowerTrackArray->Reset();
     482    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     483    {
     484      mother = track;
     485      track = static_cast<Candidate*>(track->Clone());
     486      track->AddCandidate(mother);
     487
     488      fEFlowTrackOutputArray->Add(track);
     489    }
     490  }
    508491   
    509     tower->PID = (fIsEcal) ? 22 : 0;
    510    
    511     fEFlowTowerOutputArray->Add(tower);
    512   }
    513 }
     492  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     493  else if (fTrackEnergy > 0.0)
     494  {
     495    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
     496    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
     497 
     498    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
     499    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
     500   
     501    fItTowerTrackArray->Reset();
     502    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     503    {
     504      mother = track;
     505      track = static_cast<Candidate*>(track->Clone());
     506      track->AddCandidate(mother);
     507
     508      track->Momentum *= rescaleFactor;
     509
     510      fEFlowTrackOutputArray->Add(track);
     511    }
     512  }
     513   
     514 }
    514515
    515516//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

    r425cd17 r5c402f7  
    5959  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    6060  Double_t fTowerEnergy;
    61   Double_t fTrackEnergy[2];
     61  Double_t fTrackEnergy;
    6262
    6363  Double_t fTowerTime;
    … …  
    7272
    7373  Double_t fEnergySignificanceMin;
     74
     75  Double_t fTrackSigma;
    7476
    7577  Bool_t fSmearTowerCenter;
    … …  
    102104  TObjArray *fEFlowTowerOutputArray; //!
    103105
    104   TObjArray *fTowerTrackArray[2]; //!
    105   TIterator *fItTowerTrackArray[2]; //!
     106  TObjArray *fTowerTrackArray; //!
     107  TIterator *fItTowerTrackArray; //!
    106108
    107109  void FinalizeTower();

  • modules/VertexFinderDA4D.cc

    r425cd17 r5c402f7  
    560560  cout.precision(4);
    561561  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    562     cout  <<  setw(8) << fixed << k->z;
     562    //cout  <<  setw(8) << fixed << k->z;
    563563  }
    564564  cout << endl << "                                                               t= ";
    565565  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    566     cout  <<  setw(8) << fixed << k->t;
    567   }
    568   cout << endl << "T=" << setw(15) << 1./beta <<"                                             Tc= ";
     566    //cout  <<  setw(8) << fixed << k->t;
     567  }
     568  //cout << endl << "T=" << setw(15) << 1./beta <<"                                             Tc= ";
    569569  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    570     cout  <<  setw(8) << fixed << k->Tc ;
     570    //cout  <<  setw(8) << fixed << k->Tc ;
    571571  }
    572572
    … …  
    574574  double sumpk=0;
    575575  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    576     cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
     576    //cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
    577577    sumpk+=k->pk;
    578578  }
    … …  
    588588      double tz= tks[i].z;
    589589      double tt= tks[i].t;
    590       cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
    591            << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
     590      //cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
     591      //     << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
    592592
    593593      double sump=0.;
    … …  
    597597          double p=k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
    598598          if( p > 0.0001){
    599             cout <<  setw (8) <<  setprecision(3) << p;
     599            //cout <<  setw (8) <<  setprecision(3) << p;
    600600          }else{
    601601            cout << "    .   ";

  • readers/DelphesPythia8.cpp

    r425cd17 r5c402f7  
    2525
    2626#include "Pythia.h"
     27#include "Pythia8Plugins/CombineMatchingInput.h"
    2728
    2829#include "TROOT.h"
    … …  
    224225
    225226  Pythia8::Pythia *pythia = 0;
     227 
     228  // for matching
     229  Pythia8::CombineMatchingInput *combined = 0;
     230  Pythia8::UserHooks* matching = 0;
    226231
    227232  if(argc != 4)
    … …  
    270275    // Initialize Pythia
    271276    pythia = new Pythia8::Pythia;
    272     //Pythia8::Event& event = pythia->event;
    273     //Pythia8::ParticleData& pdt = pythia->particleData;
     277 
     278    // jet matching
     279    matching = combined->getHook(*pythia);
     280    if (!matching)
     281    {
     282      throw runtime_error("can't do matching");
     283    }
     284    pythia->setUserHooksPtr(matching);
     285 
    274286
    275287    if(pythia == NULL)
Note: See TracChangeset for help on using the changeset viewer.