Fork me on GitHub

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


Ignore:
Files:
3 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • external/PUPPI/PuppiAlgo.cc

    r5c402f7 r425cd17  
     1//#include "FWCore/ParameterSet/interface/ParameterSet.h"
     2//#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
    13#include "PuppiAlgo.hh"
     4#include "external/fastjet/internal/base.hh"
    25#include "Math/QuantFuncMathCore.h"
    36#include "Math/SpecFuncMathCore.h"
    … …  
    1619  std::vector<AlgoSubObj> lAlgos = iAlgo.subAlgos;
    1720  fNAlgos = lAlgos.size();
    18   //std::cout << "==> "  << fEtaMin << " - " << fEtaMax << " - " << fPtMin  << " - " << fNeutralPtMin  << " - " << fNeutralPtSlope  << " - " << fRMSEtaSF  << " - " << std::endl;
     21  std::cout << "==> "  << fEtaMin << " - " << fEtaMax << " - " << fPtMin  << " - " << fNeutralPtMin  << " - " << fNeutralPtSlope  << " - " << fRMSEtaSF  << " - " << std::endl;
    1922  for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++)  {
    2023    int    pAlgoId      = lAlgos[i0].metricId;
    … …  
    2528    double pRMSPtMin    = lAlgos[i0].rmsPtMin;
    2629    double pRMSSF       = lAlgos[i0].rmsScaleFactor;
    27     //std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
     30    std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
    2831    fAlgoId        .push_back(pAlgoId);
    2932    fCharged       .push_back(pCharged);

  • external/PUPPI/PuppiAlgo.hh

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

  • external/PUPPI/PuppiContainer.cc

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

  • external/PUPPI/PuppiContainer.hh

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

  • external/PUPPI/puppiCleanContainer.hh

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

  • modules/Calorimeter.cc

    r5c402f7 r425cd17  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60  
     60  Int_t i;
     61
    6162  fECalResolutionFormula = new DelphesFormula;
    6263  fHCalResolutionFormula = new DelphesFormula;
    6364
    64   fECalTowerTrackArray = new TObjArray;
    65   fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
    66 
    67   fHCalTowerTrackArray = new TObjArray;
    68   fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
    69  
     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  }
    7073}
    7174
    … …  
    7477Calorimeter::~Calorimeter()
    7578{
    76  
     79  Int_t i;
     80
    7781  if(fECalResolutionFormula) delete fECalResolutionFormula;
    7882  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    7983
    80   if(fECalTowerTrackArray) delete fECalTowerTrackArray;
    81   if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
    82 
    83   if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
    84   if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
    85  
     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  }
    8692}
    8793
    … …  
    212218  Double_t ecalFraction, hcalFraction;
    213219  Double_t ecalEnergy, hcalEnergy;
     220  Double_t ecalSigma, hcalSigma;
    214221  Int_t pdgCode;
    215222
    … …  
    361368      fHCalTowerEnergy = 0.0;
    362369
    363       fECalTrackEnergy = 0.0;
    364       fHCalTrackEnergy = 0.0;
    365      
    366       fECalTrackSigma = 0.0;
    367       fHCalTrackSigma = 0.0;
    368      
     370      fECalTrackEnergy[0] = 0.0;
     371      fECalTrackEnergy[1] = 0.0;
     372
     373      fHCalTrackEnergy[0] = 0.0;
     374      fHCalTrackEnergy[1] = 0.0;
     375
    369376      fTowerTrackHits = 0;
    370377      fTowerPhotonHits = 0;
    371378
    372       fECalTowerTrackArray->Clear();
    373       fHCalTowerTrackArray->Clear();
    374    
     379      fECalTowerTrackArray[0]->Clear();
     380      fECalTowerTrackArray[1]->Clear();
     381
     382      fHCalTowerTrackArray[0]->Clear();
     383      fHCalTowerTrackArray[1]->Clear();
    375384    }
    376385
    … …  
    397406      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    398407      {
    399         fECalTrackEnergy += ecalEnergy;
    400         fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
    401         fECalTowerTrackArray->Add(track);
     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        }
    402419      }
    403      
    404420      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    405421      {
    406         fHCalTrackEnergy += hcalEnergy;
    407         fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
    408         fHCalTowerTrackArray->Add(track);
     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        }
    409433      }
    410      
    411434      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    412435      {
    … …  
    453476  Double_t energy, pt, eta, phi;
    454477  Double_t ecalEnergy, hcalEnergy;
    455   Double_t ecalNeutralEnergy, hcalNeutralEnergy;
    456  
    457478  Double_t ecalSigma, hcalSigma;
    458   Double_t ecalNeutralSigma, hcalNeutralSigma;
    459 
    460   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    461  
     479
    462480  TLorentzVector momentum;
    463481  TFractionMap::iterator itFractionMap;
    … …  
    466484
    467485  if(!fTower) return;
     486
    468487
    469488  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    535554    fTowerOutputArray->Add(fTower);
    536555  }
    537  
     556
    538557  // fill energy flow candidates
    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)
     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)
    551618  {
    552619    // create new photon tower
    553620    tower = static_cast<Candidate*>(fTower->Clone());
    554     pt =  ecalNeutralEnergy / TMath::CosH(eta);
    555    
    556     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
    557     tower->Eem = ecalNeutralEnergy;
     621
     622    pt = ecalEnergy / TMath::CosH(eta);
     623
     624    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
     625    tower->Eem = ecalEnergy;
    558626    tower->Ehad = 0.0;
    559627    tower->PID = 22;
    560    
     628
    561629    fEFlowPhotonOutputArray->Add(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
     630  }
     631  if(hcalEnergy > 0.0)
     632  {
     633    // create new neutral hadron tower
    604634    tower = static_cast<Candidate*>(fTower->Clone());
    605     pt =  hcalNeutralEnergy / TMath::CosH(eta);
    606    
    607     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
    608     tower->Ehad = hcalNeutralEnergy;
     635
     636    pt = hcalEnergy / TMath::CosH(eta);
     637
     638    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
    609639    tower->Eem = 0.0;
    610    
     640    tower->Ehad = hcalEnergy;
     641
    611642    fEFlowNeutralHadronOutputArray->Add(tower);
    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  
     643  }
    650644}
    651645

  • modules/Calorimeter.h

    r5c402f7 r425cd17  
    5858  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    5959  Double_t fECalTowerEnergy, fHCalTowerEnergy;
    60   Double_t fECalTrackEnergy, fHCalTrackEnergy;
     60  Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
    6161
    6262  Double_t fTimingEnergyMin;
    … …  
    7070  Double_t fECalEnergySignificanceMin;
    7171  Double_t fHCalEnergySignificanceMin;
    72 
    73   Double_t fECalTrackSigma;
    74   Double_t fHCalTrackSigma;
    7572
    7673  Bool_t fSmearTowerCenter;
    … …  
    106103  TObjArray *fEFlowNeutralHadronOutputArray; //!
    107104
    108   TObjArray *fECalTowerTrackArray; //!
    109   TIterator *fItECalTowerTrackArray; //!
     105  TObjArray *fECalTowerTrackArray[2]; //!
     106  TIterator *fItECalTowerTrackArray[2]; //!
    110107
    111   TObjArray *fHCalTowerTrackArray; //!
    112   TIterator *fItHCalTowerTrackArray; //!
     108  TObjArray *fHCalTowerTrackArray[2]; //!
     109  TIterator *fItHCalTowerTrackArray[2]; //!
    113110
    114111  void FinalizeTower();

  • modules/RunPUPPI.cc

    r5c402f7 r425cd17  
    33#include "PUPPI/RecoObj2.hh"
    44#include "PUPPI/AlgoObj.hh"
    5 #include "PUPPI/PuppiContainer.hh"
    6 
    7 #include "fastjet/PseudoJet.hh"
     5//#include "PUPPI/puppiParticle.hh"
     6//#include "PUPPI/puppiAlgoBin.hh"
    87
    98#include "classes/DelphesClasses.h"
    … …  
    1817
    1918using namespace std;
    20 using namespace fastjet;
    2119
    2220//------------------------------------------------------------------------------
    … …  
    240238  fPuppi->initialize(puppiInputVector);
    241239  fPuppi->puppiWeights();
    242   std::vector<PseudoJet> puppiParticles = fPuppi->puppiParticles();
     240  std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
    243241
    244242  // Loop on final particles
    245   for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
     243  for (std::vector<fastjet::PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
    246244    if(it->user_index() <= int(InputParticles.size())){     
    247245      candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());

  • modules/RunPUPPI.h

    r5c402f7 r425cd17  
    33
    44#include "classes/DelphesModule.h"
     5#include "PUPPI/PuppiContainer.hh"
    56#include <vector>
    67
    78class TObjArray;
    89class TIterator;
    9 class PuppiContainer;
     10
    1011
    1112class RunPUPPI: public DelphesModule {

  • modules/SimpleCalorimeter.cc

    r5c402f7 r425cd17  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60  
     60  Int_t i;
     61
    6162  fResolutionFormula = new DelphesFormula;
    62   fTowerTrackArray = new TObjArray;
    63   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    64  
     63
     64  for(i = 0; i < 2; ++i)
     65  {
     66    fTowerTrackArray[i] = new TObjArray;
     67    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
     68  }
    6569}
    6670
    … …  
    6973SimpleCalorimeter::~SimpleCalorimeter()
    7074{
    71  
     75  Int_t i;
     76
    7277  if(fResolutionFormula) delete fResolutionFormula;
    73   if(fTowerTrackArray) delete fTowerTrackArray;
    74   if(fItTowerTrackArray) delete fItTowerTrackArray;
    75  
     78
     79  for(i = 0; i < 2; ++i)
     80  {
     81    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
     82    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
     83  }
    7684}
    7785
    … …  
    190198  Double_t fraction;
    191199  Double_t energy;
     200  Double_t sigma;
    192201  Int_t pdgCode;
    193202
    … …  
    331340      fTowerEnergy = 0.0;
    332341
    333       fTrackEnergy = 0.0;
    334       fTrackSigma = 0.0;
     342      fTrackEnergy[0] = 0.0;
     343      fTrackEnergy[1] = 0.0;
    335344
    336345      fTowerTime = 0.0;
    … …  
    342351      fTowerPhotonHits = 0;
    343352
    344       fTowerTrackArray->Clear();
    345      }
     353      fTowerTrackArray[0]->Clear();
     354      fTowerTrackArray[1]->Clear();
     355    }
    346356
    347357    // check for track hits
    … …  
    356366      energy = momentum.E() * fTrackFractions[number];
    357367
    358       fTrackTime += TMath::Sqrt(energy)*position.T();
    359       fTrackTimeWeight += TMath::Sqrt(energy);
     368      fTrackTime += energy*position.T();
     369      fTrackTimeWeight += energy;
    360370
    361371      if(fTrackFractions[number] > 1.0E-9)
    362372      {
    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      
     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        }
    370384      }
    371        
    372385      else
    373386      {
    … …  
    405418{
    406419  Candidate *tower, *track, *mother;
    407   Double_t energy,neutralEnergy, pt, eta, phi;
    408   Double_t sigma, neutralSigma;
     420  Double_t energy, pt, eta, phi;
     421  Double_t sigma;
    409422  Double_t time;
    410    
    411   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    412423
    413424  TLorentzVector momentum;
    … …  
    425436
    426437  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    427 
    428438
    429439  if(fSmearTowerCenter)
    … …  
    454464  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    455465
    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)
     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)
    469499  {
    470500    // create new photon tower
    471501    tower = static_cast<Candidate*>(fTower->Clone());
    472     pt = neutralEnergy / TMath::CosH(eta);
    473 
    474     tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
    475     tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     502    pt = energy / TMath::CosH(eta);
     503
     504    tower->Eem = (!fIsEcal) ? 0 : energy;
     505    tower->Ehad = (fIsEcal) ? 0 : energy;
     506   
     507    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     508   
    476509    tower->PID = (fIsEcal) ? 22 : 0;
    477  
    478     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     510   
    479511    fEFlowTowerOutputArray->Add(tower);
    480    
    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   }
    491    
    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  }
     512  }
     513}
    515514
    516515//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

    r5c402f7 r425cd17  
    5959  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    6060  Double_t fTowerEnergy;
    61   Double_t fTrackEnergy;
     61  Double_t fTrackEnergy[2];
    6262
    6363  Double_t fTowerTime;
    … …  
    7272
    7373  Double_t fEnergySignificanceMin;
    74 
    75   Double_t fTrackSigma;
    7674
    7775  Bool_t fSmearTowerCenter;
    … …  
    104102  TObjArray *fEFlowTowerOutputArray; //!
    105103
    106   TObjArray *fTowerTrackArray; //!
    107   TIterator *fItTowerTrackArray; //!
     104  TObjArray *fTowerTrackArray[2]; //!
     105  TIterator *fItTowerTrackArray[2]; //!
    108106
    109107  void FinalizeTower();

  • modules/VertexFinderDA4D.cc

    r5c402f7 r425cd17  
    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

    r5c402f7 r425cd17  
    2525
    2626#include "Pythia.h"
    27 #include "Pythia8Plugins/CombineMatchingInput.h"
    2827
    2928#include "TROOT.h"
    … …  
    225224
    226225  Pythia8::Pythia *pythia = 0;
    227  
    228   // for matching
    229   Pythia8::CombineMatchingInput *combined = 0;
    230   Pythia8::UserHooks* matching = 0;
    231226
    232227  if(argc != 4)
    … …  
    275270    // Initialize Pythia
    276271    pythia = new Pythia8::Pythia;
    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  
     272    //Pythia8::Event& event = pythia->event;
     273    //Pythia8::ParticleData& pdt = pythia->particleData;
    286274
    287275    if(pythia == NULL)
Note: See TracChangeset for help on using the changeset viewer.