Fork me on GitHub

Changes in / [532016b:9be57d9] in git


Ignore:
Files:
2 added
18 edited

Legend:

Unmodified
Added
Removed
  • classes/DelphesClasses.cc

    r532016b r9be57d9  
    120120  PID(0), Status(0), M1(-1), M2(-1), D1(-1), D2(-1),
    121121  Charge(0), Mass(0.0),
    122   IsPU(0), IsConstituent(0),
     122  IsPU(0), IsRecoPU(0), IsConstituent(0),
    123123  BTag(0), TauTag(0), Eem(0.0), Ehad(0.0),
    124124  DeltaEta(0.0), DeltaPhi(0.0),
     
    133133  MeanSqDeltaR(0),
    134134  PTD(0),
     135  Ntimes(-1),
     136  IsolationVar(-999),
     137  IsolationVarRhoCorr(-999),
     138  SumPtCharged(-999),
     139  SumPtNeutral(-999),
     140  SumPtChargedPU(-999),
     141  SumPt(-999),
    135142  fFactory(0),
    136143  fArray(0)
     
    249256  object.MeanSqDeltaR = MeanSqDeltaR;
    250257  object.PTD = PTD;
     258  object.Ntimes = Ntimes;
     259  object.IsolationVar = IsolationVar;
     260  object.IsolationVarRhoCorr = IsolationVarRhoCorr;
     261  object.SumPtCharged = SumPtCharged;
     262  object.SumPtNeutral = SumPtNeutral;
     263  object.SumPtChargedPU = SumPtChargedPU;
     264  object.SumPt = SumPt;
     265
    251266  object.FracPt[0] = FracPt[0];
    252267  object.FracPt[1] = FracPt[1];
     
    262277  object.fFactory = fFactory;
    263278  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));
    264283
    265284  if(fArray && fArray->GetEntriesFast() > 0)
     
    311330  MeanSqDeltaR = 0.0;
    312331  PTD = 0.0;
     332 
     333  Ntimes = 0;
     334  Ecal_E_t.clear();
     335 
     336  IsolationVar = -999;
     337  IsolationVarRhoCorr = -999;
     338  SumPtCharged = -999;
     339  SumPtNeutral = -999;
     340  SumPtChargedPU = -999;
     341  SumPt = -999;
     342 
    313343  FracPt[0] = 0.0;
    314344  FracPt[1] = 0.0;
  • classes/DelphesClasses.h

    r532016b r9be57d9  
    231231  TRefArray Particles; // references to generated particles
    232232
     233  // Isolation variables
     234 
     235  Float_t IsolationVar;
     236  Float_t IsolationVarRhoCorr;
     237  Float_t SumPtCharged;
     238  Float_t SumPtNeutral;
     239  Float_t SumPtChargedPU;
     240  Float_t SumPt;
     241
    233242  static CompBase *fgCompare; //!
    234243  const CompBase *GetCompare() const { return fgCompare; }
     
    257266  TRef Particle; // reference to generated particle
    258267
     268  // Isolation variables
     269 
     270  Float_t IsolationVar;
     271  Float_t IsolationVarRhoCorr;
     272  Float_t SumPtCharged;
     273  Float_t SumPtNeutral;
     274  Float_t SumPtChargedPU;
     275  Float_t SumPt;
     276
    259277  static CompBase *fgCompare; //!
    260278  const CompBase *GetCompare() const { return fgCompare; }
     
    280298
    281299  TRef Particle; // reference to generated particle
     300
     301   // Isolation variables
     302 
     303  Float_t IsolationVar;
     304  Float_t IsolationVarRhoCorr;
     305  Float_t SumPtCharged;
     306  Float_t SumPtNeutral;
     307  Float_t SumPtChargedPU;
     308  Float_t SumPt;
    282309
    283310  static CompBase *fgCompare; //!
     
    334361
    335362  TLorentzVector P4();
     363  TLorentzVector Area;
    336364
    337365  ClassDef(Jet, 2)
     
    392420  Float_t E; // calorimeter tower energy
    393421
    394   Float_t T; //particle arrival time of flight
    395 
     422  Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
     423  Int_t   Ntimes; // number of hits contributing to time measurement
     424 
    396425  Float_t Eem; // calorimeter tower electromagnetic energy
    397426  Float_t Ehad; // calorimeter tower hadronic energy
     
    452481
    453482  Int_t IsPU;
     483  Int_t IsRecoPU;
     484 
    454485  Int_t IsConstituent;
    455486
     
    481512  Float_t  PTD;
    482513  Float_t  FracPt[5];
     514 
     515  //Timing information
     516 
     517  Int_t    Ntimes;
     518  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t;
     519
     520  // Isolation variables
     521 
     522  Float_t IsolationVar;
     523  Float_t IsolationVarRhoCorr;
     524  Float_t SumPtCharged;
     525  Float_t SumPtNeutral;
     526  Float_t SumPtChargedPU;
     527  Float_t SumPt;
    483528
    484529  // N-subjettiness variables
    485530
    486531  Float_t Tau[5];
     532 
     533 
     534 
    487535
    488536  static CompBase *fgCompare; //!
  • modules/Calorimeter.cc

    r532016b r9be57d9  
    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

    r532016b r9be57d9  
    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

    r532016b r9be57d9  
    2525 *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
    2626 *
    27  *  \author P. Demin - UCL, Louvain-la-Neuve
     27 *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
    2828 *
    2929 */
     
    109109  fUsePTSum = GetBool("UsePTSum", false);
    110110
     111  fVetoLeptons = GetBool("VetoLeptons", true); 
     112 
    111113  fClassifier->fPTMin = GetDouble("PTMin", 0.5);
     114
    112115
    113116  // import input array(s)
     
    153156  Candidate *candidate, *isolation, *object;
    154157  TObjArray *isolationArray;
    155   Double_t sum, ratio;
     158  Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
    156159  Int_t counter;
    157160  Double_t eta = 0.0;
     
    180183
    181184    // loop over all input tracks
    182     sum = 0.0;
     185   
     186    sumNeutral = 0.0;
     187    sumCharged = 0.0;
     188    sumChargedPU = 0.0;
     189    sumAllParticles = 0.0;
     190   
    183191    counter = 0;
    184192    itIsolationArray.Reset();
     193   
    185194    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
    186195    {
     
    188197
    189198      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
    190          candidate->GetUniqueID() != isolation->GetUniqueID())
     199         candidate->GetUniqueID() != isolation->GetUniqueID() &&
     200         ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) )
    191201      {
    192         sum += isolationMomentum.Pt();
     202     
     203        sumAllParticles += isolationMomentum.Pt();
     204        if(isolation->Charge !=0)
     205        {
     206          sumCharged += isolationMomentum.Pt();
     207          if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt();
     208        }
     209 
     210        else sumNeutral += isolationMomentum.Pt();
     211 
    193212        ++counter;
    194213      }
     
    209228    }
    210229
    211     // correct sum for pile-up contamination
    212     sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
    213 
    214     ratio = sum/candidateMomentum.Pt();
    215     if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
    216 
     230       
     231     // correct sum for pile-up contamination
     232    sumDBeta = sumCharged + TMath::Max(sumNeutral-0.5*sumChargedPU,0.0);
     233    sumRhoCorr = sumCharged + TMath::Max(sumNeutral-TMath::Max(rho,0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(),0.0);
     234    ratioDBeta = sumDBeta/candidateMomentum.Pt();
     235    ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
     236   
     237    candidate->IsolationVar = ratioDBeta;
     238    candidate->IsolationVarRhoCorr = ratioRhoCorr;
     239    candidate->SumPtCharged = sumCharged;
     240    candidate->SumPtNeutral = sumNeutral;
     241    candidate->SumPtChargedPU = sumChargedPU;
     242    candidate->SumPt = sumAllParticles;
     243
     244    if((fUsePTSum && sumDBeta > fPTSumMax) || ratioDBeta > fPTRatioMax) continue;
    217245    fOutputArray->Add(candidate);
    218246  }
  • modules/Isolation.h

    r532016b r9be57d9  
    5959  Bool_t fUsePTSum;
    6060
     61  Bool_t fVetoLeptons; 
     62 
    6163  IsolationClassifier *fClassifier; //!
    6264
  • modules/ModulesLinkDef.h

    r532016b r9be57d9  
    5050#include "modules/JetPileUpSubtractor.h"
    5151#include "modules/TrackPileUpSubtractor.h"
     52#include "modules/TaggingParticlesSkimmer.h"
    5253#include "modules/PileUpJetID.h"
    5354#include "modules/ConstituentFilter.h"
     
    8990#pragma link C++ class JetPileUpSubtractor+;
    9091#pragma link C++ class TrackPileUpSubtractor+;
     92#pragma link C++ class TaggingParticlesSkimmer+;
    9193#pragma link C++ class PileUpJetID+;
    9294#pragma link C++ class ConstituentFilter+;
  • modules/PileUpJetID.cc

    r532016b r9be57d9  
    1 /*
    2  *  Delphes: a framework for fast simulation of a generic collider experiment
    3  *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    4  *
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  *
    10  *  This program is distributed in the hope that it will be useful,
    11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
    12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    13  *  GNU General Public License for more details.
    14  *
    15  *  You should have received a copy of the GNU General Public License
    16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
    17  */
    18 
    191
    202/** \class PileUpJetID
    213 *
    22  *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
     4 *  CMS PileUp Jet ID Variables
    235 *
    24  *  \author S. Zenz, December 2013
     6 *  \author S. Zenz
    257 *
    268 */
     
    4123#include "TRandom3.h"
    4224#include "TObjArray.h"
    43 #include "TDatabasePDG.h"
     25//#include "TDatabasePDG.h"
    4426#include "TLorentzVector.h"
    4527
     
    5436
    5537PileUpJetID::PileUpJetID() :
    56   fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
     38  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
    5739{
    5840
     
    7052void PileUpJetID::Init()
    7153{
     54
    7255  fJetPTMin = GetDouble("JetPTMin", 20.0);
    7356  fParameterR = GetDouble("ParameterR", 0.5);
    7457  fUseConstituents = GetInt("UseConstituents", 0);
    7558
     59
     60  /*
     61  Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
     62  Double_t fBetaMinBarrel; // |eta| < 2.5
     63  Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
     64  Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
     65  Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
     66  */
     67
     68  fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
     69  fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
     70  fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
     71  fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
     72  fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
     73  fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
     74  fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
     75
     76
     77
     78  cout << "  set MeanSqDeltaRMaxBarrel " << fMeanSqDeltaRMaxBarrel << endl;
     79  cout << "  set BetaMinBarrel " << fBetaMinBarrel << endl;
     80  cout << "  set MeanSqDeltaRMaxEndcap " << fMeanSqDeltaRMaxEndcap << endl;
     81  cout << "  set BetaMinEndcap " << fBetaMinEndcap << endl;
     82  cout << "  set MeanSqDeltaRMaxForward " << fMeanSqDeltaRMaxForward << endl;
     83
     84
     85
    7686  fAverageEachTower = false; // for timing
    7787
     
    8191  fItJetInputArray = fJetInputArray->MakeIterator();
    8292
    83   fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
     93
     94  //  cout << "BeforE SCZ additions in init" << endl;
     95  //  cout << GetString("TrackInputArray", "ParticlePropagator/tracks") << endl;
     96  //  cout << GetString("EFlowTrackInputArray", "ParticlePropagator/tracks") << endl;
     97
     98  fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
    8499  fItTrackInputArray = fTrackInputArray->MakeIterator();
    85100
    86   fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
     101  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
    87102  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    88103
    89   fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
    90   fItVertexInputArray = fVertexInputArray->MakeIterator();
    91 
    92   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    93104
    94105  // create output array(s)
    95106
    96107  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
     108
     109  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
     110
     111
     112  //  cout << " end of INIT " << endl;
     113
    97114}
    98115
     
    101118void PileUpJetID::Finish()
    102119{
     120  //  cout << "In finish" << endl;
     121
    103122  if(fItJetInputArray) delete fItJetInputArray;
    104123  if(fItTrackInputArray) delete fItTrackInputArray;
    105124  if(fItNeutralInputArray) delete fItNeutralInputArray;
    106   if(fItVertexInputArray) delete fItVertexInputArray;
     125
    107126}
    108127
     
    111130void PileUpJetID::Process()
    112131{
     132  //  cout << "start of process" << endl;
     133
    113134  Candidate *candidate, *constituent;
    114135  TLorentzVector momentum, area;
    115   Int_t i, nc, nn;
    116   Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
    117   Double_t dr, pt, pt_ann[5];
    118   Double_t zvtx = 0.0;
    119 
    120   Candidate *track;
    121 
    122   // find z position of primary vertex
    123 
    124   fItVertexInputArray->Reset();
    125   while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
    126   {
    127     if(!candidate->IsPU)
    128     {
    129       zvtx = candidate->Position.Z();
    130       break;
    131     }
    132   }
     136
     137  //  cout << "BeforE SCZ additions in process" << endl;
     138
     139  // SCZ
     140  Candidate *trk;
    133141
    134142  // loop over all input candidates
     
    139147    area = candidate->Area;
    140148
    141     sumpt = 0.0;
    142     sumptch = 0.0;
    143     sumptchpv = 0.0;
    144     sumptchpu = 0.0;
    145     sumdrsqptsq = 0.0;
    146     sumptsq = 0.0;
    147     nc = 0;
    148     nn = 0;
    149 
    150     for(i = 0; i < 5; ++i)
    151     {
    152       pt_ann[i] = 0.0;
    153     }
    154 
    155     if(fUseConstituents)
    156     {
     149    float sumT0 = 0.;
     150    float sumT1 = 0.;
     151    float sumT10 = 0.;
     152    float sumT20 = 0.;
     153    float sumT30 = 0.;
     154    float sumT40 = 0.;
     155    float sumWeightsForT = 0.;
     156    candidate->Ntimes = 0;
     157
     158    float sumpt = 0.;
     159    float sumptch = 0.;
     160    float sumptchpv = 0.;
     161    float sumptchpu = 0.;
     162    float sumdrsqptsq = 0.;
     163    float sumptsq = 0.;
     164    int nc = 0;
     165    int nn = 0;
     166    float pt_ann[5];
     167
     168    for (int i = 0 ; i < 5 ; i++) {
     169      pt_ann[i] = 0.;
     170    }
     171
     172    if (fUseConstituents) {
    157173      TIter itConstituents(candidate->GetCandidates());
    158       while((constituent = static_cast<Candidate*>(itConstituents.Next())))
    159       {
    160         pt = constituent->Momentum.Pt();
    161         dr = candidate->Momentum.DeltaR(constituent->Momentum);
    162         sumpt += pt;
    163         sumdrsqptsq += dr*dr*pt*pt;
    164         sumptsq += pt*pt;
    165         if(constituent->Charge == 0)
    166         {
    167           // neutrals
    168           ++nn;
    169         }
    170         else
    171         {
    172           // charged
    173           if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
    174           {
    175             sumptchpu += pt;
    176           }
    177           else
    178           {
    179             sumptchpv += pt;
    180           }
    181           sumptch += pt;
    182           ++nc;
    183         }
    184         for(i = 0; i < 5; ++i)
    185         {
    186           if(dr > 0.1*i && dr < 0.1*(i + 1))
    187           {
    188             pt_ann[i] += pt;
    189           }
    190         }
    191       }
    192     }
    193     else
    194     {
     174      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
     175        float pt = constituent->Momentum.Pt();
     176        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
     177        //      cout << " There exists a constituent with dr=" << dr << endl;
     178        sumpt += pt;
     179        sumdrsqptsq += dr*dr*pt*pt;
     180        sumptsq += pt*pt;
     181        if (constituent->Charge == 0) {
     182          nn++;
     183        } else {
     184          if (constituent->IsRecoPU) {
     185            sumptchpu += pt;
     186          } else {
     187            sumptchpv += pt;
     188          }
     189          sumptch += pt;
     190          nc++;
     191        }
     192        for (int i = 0 ; i < 5 ; i++) {
     193          if (dr > 0.1*i && dr < 0.1*(i+1)) {
     194            pt_ann[i] += pt;
     195          }
     196        }
     197        float tow_sumT = 0;
     198        float tow_sumW = 0;
     199        for (int i = 0 ; i < constituent->Ecal_E_t.size() ; i++) {
     200          float w = TMath::Sqrt(constituent->Ecal_E_t[i].first);
     201          if (fAverageEachTower) {
     202            tow_sumT += w*constituent->Ecal_E_t[i].second;
     203            tow_sumW += w;
     204          } else {
     205            sumT0 += w*constituent->Ecal_E_t[i].second;
     206            sumT1 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.001);
     207            sumT10 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.010);
     208            sumT20 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.020);
     209            sumT30 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.030);
     210            sumT40 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.040);
     211            sumWeightsForT += w;
     212            candidate->Ntimes++;
     213          }
     214        }
     215        if (fAverageEachTower && tow_sumW > 0.) {
     216          sumT0 += tow_sumT;
     217          sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
     218          sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
     219          sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
     220          sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
     221          sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
     222          sumWeightsForT += tow_sumW;
     223          candidate->Ntimes++;
     224        }
     225      }
     226    } else {
    195227      // Not using constituents, using dr
    196228      fItTrackInputArray->Reset();
    197       while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
    198       {
    199         if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
    200         {
    201           pt = track->Momentum.Pt();
    202           sumpt += pt;
    203           sumptch += pt;
    204           if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
    205           {
    206             sumptchpu += pt;
    207           }
    208           else
    209           {
    210             sumptchpv += pt;
    211           }
    212           dr = candidate->Momentum.DeltaR(track->Momentum);
    213           sumdrsqptsq += dr*dr*pt*pt;
    214           sumptsq += pt*pt;
    215           nc++;
    216           for(i = 0; i < 5; ++i)
    217           {
    218             if(dr > 0.1*i && dr < 0.1*(i + 1))
    219             {
    220               pt_ann[i] += pt;
    221             }
    222           }
    223         }
    224       }
    225 
     229      while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
     230        if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
     231          float pt = trk->Momentum.Pt();
     232          sumpt += pt;
     233          sumptch += pt;
     234          if (trk->IsRecoPU) {
     235            sumptchpu += pt;
     236          } else {
     237            sumptchpv += pt;
     238          }
     239          float dr = candidate->Momentum.DeltaR(trk->Momentum);
     240          sumdrsqptsq += dr*dr*pt*pt;
     241          sumptsq += pt*pt;
     242          nc++;
     243          for (int i = 0 ; i < 5 ; i++) {
     244            if (dr > 0.1*i && dr < 0.1*(i+1)) {
     245              pt_ann[i] += pt;
     246            }
     247          }
     248        }
     249      }
    226250      fItNeutralInputArray->Reset();
    227       while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
    228       {
    229         if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
    230         {
    231           pt = constituent->Momentum.Pt();
    232           sumpt += pt;
    233           dr = candidate->Momentum.DeltaR(constituent->Momentum);
    234           sumdrsqptsq += dr*dr*pt*pt;
    235           sumptsq += pt*pt;
    236           nn++;
    237           for(i = 0; i < 5; ++i)
    238           {
    239             if(dr > 0.1*i && dr < 0.1*(i + 1))
    240             {
    241               pt_ann[i] += pt;
    242             }
    243           }
    244         }
    245       }
    246     }
    247 
    248     if(sumptch > 0.0)
    249     {
    250       candidate->Beta = sumptchpu/sumptch;
    251       candidate->BetaStar = sumptchpv/sumptch;
    252     }
    253     else
    254     {
    255       candidate->Beta = -999.0;
    256       candidate->BetaStar = -999.0;
    257     }
    258     if(sumptsq > 0.0)
    259     {
     251      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
     252        if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
     253          float pt = constituent->Momentum.Pt();
     254          sumpt += pt;
     255          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
     256          sumdrsqptsq += dr*dr*pt*pt;
     257          sumptsq += pt*pt;
     258          nn++;
     259          for (int i = 0 ; i < 5 ; i++) {
     260            if (dr > 0.1*i && dr < 0.1*(i+1)) {
     261            pt_ann[i] += pt;
     262            }
     263          }
     264        }
     265      }
     266    }
     267
     268    if (sumptch > 0.) {
     269      candidate->Beta = sumptchpv/sumptch;
     270      candidate->BetaStar = sumptchpu/sumptch;
     271    } else {
     272      candidate->Beta = -999.;
     273      candidate->BetaStar = -999.;
     274    }
     275    if (sumptsq > 0.) {
    260276      candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
    261     }
    262     else
    263     {
    264       candidate->MeanSqDeltaR = -999.0;
     277    } else {
     278      candidate->MeanSqDeltaR = -999.;
    265279    }
    266280    candidate->NCharged = nc;
    267281    candidate->NNeutrals = nn;
    268     if(sumpt > 0.0)
    269     {
     282    if (sumpt > 0.) {
    270283      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
    271       for(i = 0; i < 5; ++i)
    272       {
     284      for (int i = 0 ; i < 5 ; i++) {
    273285        candidate->FracPt[i] = pt_ann[i]/sumpt;
    274286      }
    275     }
    276     else
    277     {
    278       candidate->PTD = -999.0;
    279       for(i = 0; i < 5; ++i)
    280       {
    281         candidate->FracPt[i] = -999.0;
     287    } else {
     288      candidate->PTD = -999.;
     289      for (int i = 0 ; i < 5 ; i++) {
     290        candidate->FracPt[i] = -999.;
    282291      }
    283292    }
    284293
    285294    fOutputArray->Add(candidate);
     295
     296    // New stuff
     297    /*
     298    fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
     299    fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
     300    fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
     301    fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
     302    fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
     303    */
     304
     305    bool passId = false;
     306    if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
     307      if (fabs(candidate->Momentum.Eta())<1.5) {
     308        passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
     309      } else if (fabs(candidate->Momentum.Eta())<4.0) {
     310        passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
     311      } else {
     312        passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
     313      }
     314    }
     315
     316    //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
     317    //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
     318
     319    if (passId) {
     320      if (fUseConstituents) {
     321        TIter itConstituents(candidate->GetCandidates());
     322        while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
     323          if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
     324            fNeutralsInPassingJets->Add(constituent);
     325            //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
     326          }
     327        }
     328      } else { // use DeltaR
     329        fItNeutralInputArray->Reset();
     330        while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
     331          if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
     332            fNeutralsInPassingJets->Add(constituent);
     333            //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
     334          }
     335        }
     336      }
     337    }
     338
     339
    286340  }
    287341}
    288342
    289343//------------------------------------------------------------------------------
    290 
  • modules/PileUpJetID.h

    r532016b r9be57d9  
    1 /*
    2  *  Delphes: a framework for fast simulation of a generic collider experiment
    3  *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    4  *
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  *
    10  *  This program is distributed in the hope that it will be useful,
    11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
    12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    13  *  GNU General Public License for more details.
    14  *
    15  *  You should have received a copy of the GNU General Public License
    16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
    17  */
    18 
    191#ifndef PileUpJetID_h
    202#define PileUpJetID_h
     
    224/** \class PileUpJetID
    235 *
    24  *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
     6 *  CMS PileUp Jet ID Variables
    257 *
    26  *  \author S. Zenz, December 2013
     8 *  \author S. Zenz
    279 *
    2810 */
     
    3416
    3517class TObjArray;
     18class DelphesFormula;
    3619
    3720class PileUpJetID: public DelphesModule
     
    5134  Double_t fParameterR;
    5235
     36  Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
     37  Double_t fBetaMinBarrel; // |eta| < 2.5
     38  Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
     39  Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
     40  Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
     41
     42  Double_t fNeutralPTMin;
     43  Double_t fJetPTMinForNeutrals;
     44
     45  /*
     46JAY
     47---
     48
     49|Eta|<1.5
     50
     51meanSqDeltaR betaStar SigEff BgdEff
     520.13 0.92 96% 8%
     530.13 0.95 97% 16%
     540.13 0.97 98% 27%
     55
     56|Eta|>1.5
     57
     58meanSqDeltaR betaStar SigEff BgdEff
     590.14 0.91 95% 15%
     600.14 0.94 97% 19%
     610.14 0.97 98% 29%
     62
     63BRYAN
     64-----
     65
     66Barrel (MeanSqDR, Beta, sig eff, bg eff):
     670.10, 0.08, 90%, 8%
     680.11, 0.12, 90%, 6%
     690.13, 0.16, 89%, 5%
     70
     71Endcap (MeanSqDR, Beta, sig eff, bg eff):
     720.07, 0.06, 89%, 4%
     730.08, 0.08, 92%, 6%
     740.09, 0.08, 95%, 10%
     750.10, 0.08, 97%, 13%
     76
     77SETH GUESSES FOR |eta| > 4.0
     78----------------------------
     79
     80MeanSqDeltaR
     810.07
     820.10
     830.14
     840.2
     85  */
     86
    5387  // If set to true, may have weird results for PFCHS
    5488  // If set to false, uses everything within dR < fParameterR even if in other jets &c.
    5589  // Results should be very similar for PF
    56   Int_t fUseConstituents;
     90  Int_t fUseConstituents; 
    5791
    5892  Bool_t fAverageEachTower;
     
    6296  const TObjArray *fJetInputArray; //!
    6397
    64   const TObjArray *fTrackInputArray; //!
    65   const TObjArray *fNeutralInputArray; //!
     98  const TObjArray *fTrackInputArray; // SCZ
     99  const TObjArray *fNeutralInputArray;
    66100
    67   TIterator *fItTrackInputArray; //!
    68   TIterator *fItNeutralInputArray; //!
     101  TIterator *fItTrackInputArray; // SCZ
     102  TIterator *fItNeutralInputArray; // SCZ
    69103
    70104  TObjArray *fOutputArray; //!
     105  TObjArray *fNeutralsInPassingJets; // SCZ
    71106
    72   TIterator *fItVertexInputArray; //!
    73   const TObjArray *fVertexInputArray; //!
    74107
    75   Double_t fZVertexResolution;
    76 
    77   ClassDef(PileUpJetID, 1)
     108  ClassDef(PileUpJetID, 2)
    78109};
    79110
    80111#endif
    81 
  • modules/PileUpMerger.cc

    r532016b r9be57d9  
    8080  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
    8181
     82  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
     83  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
     84  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
     85  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
     86
    8287  // read vertex smearing formula
    8388
     
    111116  TParticlePDG *pdgParticle;
    112117  Int_t pid;
    113   Float_t x, y, z, t;
     118  Float_t x, y, z, t, vx, vy;
    114119  Float_t px, py, pz, e;
    115120  Double_t dz, dphi, dt;
    116   Int_t numberOfEvents, event;
     121  Int_t numberOfEvents, event, numberOfParticles;
    117122  Long64_t allEntries, entry;
    118   Candidate *candidate, *vertexcandidate;
     123  Candidate *candidate, *vertex;
    119124  DelphesFactory *factory;
    120125
     
    123128  fItInputArray->Reset();
    124129
    125   // --- Deal with Primary vertex first  ------
     130  // --- Deal with primary vertex first  ------
    126131
    127132  fFunction->GetRandom2(dz, dt);
     
    129134  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    130135  dz *= 1.0E3; // necessary in order to make z in mm
    131 
     136  vx = 0.0;
     137  vy = 0.0;
     138  numberOfParticles = fInputArray->GetEntriesFast();
    132139  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    133140  {
     141    vx += candidate->Position.X();
     142    vy += candidate->Position.Y();
    134143    z = candidate->Position.Z();
    135144    t = candidate->Position.T();
     
    139148  }
    140149
     150  if(numberOfParticles > 0)
     151  {
     152    vx /= numberOfParticles;
     153    vy /= numberOfParticles;
     154  }
     155
    141156  factory = GetFactory();
    142157
    143   vertexcandidate = factory->NewCandidate();
    144   vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
    145   fVertexOutputArray->Add(vertexcandidate);
     158  vertex = factory->NewCandidate();
     159  vertex->Position.SetXYZT(vx, vy, dz, dt);
     160  fVertexOutputArray->Add(vertex);
    146161
    147162  // --- Then with pile-up vertices  ------
     
    181196    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    182197
    183     vertexcandidate = factory->NewCandidate();
    184     vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
    185     vertexcandidate->IsPU = 1;
    186 
    187     fVertexOutputArray->Add(vertexcandidate);
    188 
     198    vx = 0.0;
     199    vy = 0.0;
     200    numberOfParticles = 0;
    189201    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    190202    {
     
    204216      candidate->Momentum.RotateZ(dphi);
    205217
     218      x -= fInputBeamSpotX;
     219      y -= fInputBeamSpotY;
    206220      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
    207221      candidate->Position.RotateZ(dphi);
     222      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
     223
     224      vx += candidate->Position.X();
     225      vy += candidate->Position.Y();
     226      ++numberOfParticles;
    208227
    209228      fParticleOutputArray->Add(candidate);
    210229    }
    211   }
    212 }
    213 
    214 //------------------------------------------------------------------------------
    215 
     230
     231    if(numberOfParticles > 0)
     232    {
     233      vx /= numberOfParticles;
     234      vy /= numberOfParticles;
     235    }
     236
     237    vertex = factory->NewCandidate();
     238    vertex->Position.SetXYZT(vx, vy, dz, dt);
     239    vertex->IsPU = 1;
     240
     241    fVertexOutputArray->Add(vertex);
     242  }
     243}
     244
     245//------------------------------------------------------------------------------
  • modules/PileUpMerger.h

    r532016b r9be57d9  
    5353  Double_t fTVertexSpread;
    5454
     55  Double_t fInputBeamSpotX;
     56  Double_t fInputBeamSpotY;
     57  Double_t fOutputBeamSpotX;
     58  Double_t fOutputBeamSpotY;
     59
    5560  DelphesTF2 *fFunction; //!
    5661
  • modules/PileUpMergerPythia8.cc

    r532016b r9be57d9  
    2929#include "classes/DelphesClasses.h"
    3030#include "classes/DelphesFactory.h"
    31 #include "classes/DelphesFormula.h"
     31#include "classes/DelphesTF2.h"
    3232#include "classes/DelphesPileUpReader.h"
    3333
     
    5656
    5757PileUpMergerPythia8::PileUpMergerPythia8() :
    58   fPythia(0), fItInputArray(0)
    59 {
     58  fFunction(0), fPythia(0), fItInputArray(0)
     59{
     60  fFunction = new DelphesTF2;
    6061}
    6162
     
    6465PileUpMergerPythia8::~PileUpMergerPythia8()
    6566{
     67  delete fFunction;
    6668}
    6769
     
    7274  const char *fileName;
    7375
     76  fPileUpDistribution = GetInt("PileUpDistribution", 0);
     77
    7478  fMeanPileUp  = GetDouble("MeanPileUp", 10);
    75   fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
     79
     80  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
     81  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
     82
     83  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
     84  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
     85  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
     86  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
    7687
    7788  fPTMin = GetDouble("PTMin", 0.0);
     89
     90  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
     91  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
    7892
    7993  fileName = GetString("ConfigFile", "MinBias.cmnd");
     
    86100
    87101  // create output arrays
    88   fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
     102  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
     103  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
    89104}
    90105
     
    103118  TParticlePDG *pdgParticle;
    104119  Int_t pid, status;
    105   Float_t x, y, z, t;
     120  Float_t x, y, z, t, vx, vy;
    106121  Float_t px, py, pz, e;
    107   Double_t dz, dphi;
    108   Int_t poisson, event, i;
    109   Candidate *candidate;
     122  Double_t dz, dphi, dt;
     123  Int_t numberOfEvents, event, numberOfParticles, i;
     124  Candidate *candidate, *vertex;
    110125  DelphesFactory *factory;
    111126
     127  const Double_t c_light = 2.99792458E8;
     128
    112129  fItInputArray->Reset();
     130
     131  // --- Deal with primary vertex first  ------
     132
     133  fFunction->GetRandom2(dz, dt);
     134
     135  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
     136  dz *= 1.0E3; // necessary in order to make z in mm
     137  vx = 0.0;
     138  vy = 0.0;
     139  numberOfParticles = fInputArray->GetEntriesFast();
    113140  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    114141  {
    115     fOutputArray->Add(candidate);
     142    vx += candidate->Position.X();
     143    vy += candidate->Position.Y();
     144    z = candidate->Position.Z();
     145    t = candidate->Position.T();
     146    candidate->Position.SetZ(z + dz);
     147    candidate->Position.SetT(t + dt);
     148    fParticleOutputArray->Add(candidate);
     149  }
     150
     151  if(numberOfParticles > 0)
     152  {
     153    vx /= numberOfParticles;
     154    vy /= numberOfParticles;
    116155  }
    117156
    118157  factory = GetFactory();
    119158
    120   poisson = gRandom->Poisson(fMeanPileUp);
    121 
    122   for(event = 0; event < poisson; ++event)
     159  vertex = factory->NewCandidate();
     160  vertex->Position.SetXYZT(vx, vy, dz, dt);
     161  fVertexOutputArray->Add(vertex);
     162
     163  // --- Then with pile-up vertices  ------
     164
     165  switch(fPileUpDistribution)
     166  {
     167    case 0:
     168      numberOfEvents = gRandom->Poisson(fMeanPileUp);
     169      break;
     170    case 1:
     171      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
     172      break;
     173    default:
     174      numberOfEvents = gRandom->Poisson(fMeanPileUp);
     175      break;
     176  }
     177
     178  for(event = 0; event < numberOfEvents; ++event)
    123179  {
    124180    while(!fPythia->next());
    125181
    126     dz = gRandom->Gaus(0.0, fZVertexSpread);
     182   // --- Pile-up vertex smearing
     183
     184    fFunction->GetRandom2(dz, dt);
     185
     186    dt *= c_light*1.0E3; // necessary in order to make t in mm/c
     187    dz *= 1.0E3; // necessary in order to make z in mm
     188
    127189    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    128190
    129     for(i = 1; i < fPythia->event.size(); ++i)
     191    vx = 0.0;
     192    vy = 0.0;
     193    numberOfParticles = fPythia->event.size();
     194    for(i = 1; i < numberOfParticles; ++i)
    130195    {
    131196      Pythia8::Particle &particle = fPythia->event[i];
     
    143208      candidate->PID = pid;
    144209
    145       candidate->Status = status;
     210      candidate->Status = 1;
    146211
    147212      pdgParticle = pdg->GetParticle(pid);
     
    154219      candidate->Momentum.RotateZ(dphi);
    155220
    156       candidate->Position.SetXYZT(x, y, z + dz, t);
     221      x -= fInputBeamSpotX;
     222      y -= fInputBeamSpotY;
     223      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
    157224      candidate->Position.RotateZ(dphi);
    158 
    159       fOutputArray->Add(candidate);
     225      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
     226
     227      vx += candidate->Position.X();
     228      vy += candidate->Position.Y();
     229
     230      fParticleOutputArray->Add(candidate);
    160231    }
    161   }
    162 }
    163 
    164 //------------------------------------------------------------------------------
    165 
     232
     233    if(numberOfParticles > 0)
     234    {
     235      vx /= numberOfParticles;
     236      vy /= numberOfParticles;
     237    }
     238
     239    vertex = factory->NewCandidate();
     240    vertex->Position.SetXYZT(vx, vy, dz, dt);
     241    vertex->IsPU = 1;
     242
     243    fVertexOutputArray->Add(vertex);
     244  }
     245}
     246
     247//------------------------------------------------------------------------------
  • modules/PileUpMergerPythia8.h

    r532016b r9be57d9  
    3131
    3232class TObjArray;
     33class DelphesTF2;
    3334
    3435namespace Pythia8
     
    5051private:
    5152
     53  Int_t fPileUpDistribution;
    5254  Double_t fMeanPileUp;
     55
    5356  Double_t fZVertexSpread;
     57  Double_t fTVertexSpread;
     58
     59  Double_t fInputBeamSpotX;
     60  Double_t fInputBeamSpotY;
     61  Double_t fOutputBeamSpotX;
     62  Double_t fOutputBeamSpotY;
     63
    5464  Double_t fPTMin;
     65
     66  DelphesTF2 *fFunction; //!
    5567
    5668  Pythia8::Pythia *fPythia; //!
     
    6072  const TObjArray *fInputArray; //!
    6173
    62   TObjArray *fOutputArray; //!
     74  TObjArray *fParticleOutputArray; //!
     75  TObjArray *fVertexOutputArray; //!
    6376
    6477  ClassDef(PileUpMergerPythia8, 1)
  • modules/TauTagging.cc

    r532016b r9be57d9  
    3333#include "classes/DelphesFactory.h"
    3434#include "classes/DelphesFormula.h"
    35 
    36 #include "ExRootAnalysis/ExRootResult.h"
    37 #include "ExRootAnalysis/ExRootFilter.h"
    38 #include "ExRootAnalysis/ExRootClassifier.h"
    3935
    4036#include "TMath.h"
     
    5349using namespace std;
    5450
    55 //------------------------------------------------------------------------------
    56 
    57 class TauTaggingPartonClassifier : public ExRootClassifier
    58 {
    59 public:
    60 
    61   TauTaggingPartonClassifier(const TObjArray *array);
    62 
    63   Int_t GetCategory(TObject *object);
    64 
    65   Double_t fEtaMax, fPTMin;
    66 
    67   const TObjArray *fParticleInputArray;
    68 };
    6951
    7052//------------------------------------------------------------------------------
     
    7961{
    8062  Candidate *tau = static_cast<Candidate *>(object);
    81   Candidate *daughter = 0;
     63  Candidate *daughter1 = 0;
     64  Candidate *daughter2 = 0;
     65 
    8266  const TLorentzVector &momentum = tau->Momentum;
    83   Int_t pdgCode, i;
     67  Int_t pdgCode, i, j;
    8468
    8569  pdgCode = TMath::Abs(tau->PID);
     
    9882  for(i = tau->D1; i <= tau->D2; ++i)
    9983  {
    100     daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
    101     pdgCode = TMath::Abs(daughter->PID);
    102     if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15 || pdgCode == 24) return -1;
     84    daughter1 = static_cast<Candidate *>(fParticleInputArray->At(i));
     85    pdgCode = TMath::Abs(daughter1->PID);
     86    if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15) return -1;
     87    else if(pdgCode == 24)
     88    {
     89     if(daughter1->D1 < 0) return -1;
     90     for(j = daughter1->D1; j <= daughter1->D2; ++j)
     91     {
     92       daughter2 = static_cast<Candidate*>(fParticleInputArray->At(j));
     93       pdgCode = TMath::Abs(daughter2->PID);
     94       if(pdgCode == 11 || pdgCode == 13) return -1;
     95     }
     96       
     97    }
    10398  }
    10499
  • modules/TauTagging.h

    r532016b r9be57d9  
    3131
    3232#include "classes/DelphesModule.h"
     33#include "ExRootAnalysis/ExRootResult.h"
     34#include "ExRootAnalysis/ExRootFilter.h"
     35#include "ExRootAnalysis/ExRootClassifier.h"
    3336
    3437#include <map>
     
    7679};
    7780
     81
     82//------------------------------------------------------------------------------
     83
     84class TauTaggingPartonClassifier : public ExRootClassifier
     85{
     86public:
     87
     88  TauTaggingPartonClassifier(const TObjArray *array);
     89
     90  Int_t GetCategory(TObject *object);
     91
     92  Double_t fEtaMax, fPTMin;
     93
     94  const TObjArray *fParticleInputArray;
     95};
     96
     97
    7898#endif
  • modules/TrackPileUpSubtractor.cc

    r532016b r9be57d9  
    7474  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    7575
     76  fPTMin = GetDouble("PTMin", 0.);
    7677  // import arrays with output from other modules
    77 
     78   
    7879  ExRootConfParam param = GetParam("InputArray");
    7980  Long_t i, size;
     
    147148      // assume perfect pile-up subtraction for tracks outside fZVertexResolution
    148149     
    149       if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
    150 
    151       array->Add(candidate);
     150      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1;
     151      else
     152      {
     153         candidate->IsRecoPU = 0;
     154         if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
     155      }
    152156    }
    153157  }
  • modules/TrackPileUpSubtractor.h

    r532016b r9be57d9  
    5050  Double_t fZVertexResolution;
    5151
     52  Double_t fPTMin;
     53
    5254  std::map< TIterator *, TObjArray * > fInputMap; //!
    5355
  • modules/TreeWriter.cc

    r532016b r9be57d9  
    362362
    363363    entry->T = position.T()*1.0E-3/c_light;
     364    entry->Ntimes = candidate->Ntimes;   
    364365
    365366    FillParticles(candidate, &entry->Particles);
     
    403404    entry->T = position.T()*1.0E-3/c_light;
    404405
     406    // Isolation variables
     407       
     408    entry->IsolationVar = candidate->IsolationVar;
     409    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
     410    entry->SumPtCharged = candidate->SumPtCharged ;
     411    entry->SumPtNeutral = candidate->SumPtNeutral ;
     412    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
     413    entry->SumPt = candidate->SumPt ;
     414
    405415    entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
    406416
     
    442452    entry->T = position.T()*1.0E-3/c_light;
    443453
     454    // Isolation variables
     455       
     456    entry->IsolationVar = candidate->IsolationVar;
     457    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
     458    entry->SumPtCharged = candidate->SumPtCharged ;
     459    entry->SumPtNeutral = candidate->SumPtNeutral ;
     460    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
     461    entry->SumPt = candidate->SumPt ;
     462
     463
    444464    entry->Charge = candidate->Charge;
    445465
     
    487507
    488508    entry->T = position.T()*1.0E-3/c_light;
     509   
     510    // Isolation variables
     511       
     512    entry->IsolationVar = candidate->IsolationVar;
     513    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
     514    entry->SumPtCharged = candidate->SumPtCharged ;
     515    entry->SumPtNeutral = candidate->SumPtNeutral ;
     516    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
     517    entry->SumPt = candidate->SumPt ;
    489518
    490519    entry->Charge = candidate->Charge;
     
    531560
    532561    entry->Mass = momentum.M();
     562
     563    entry->Area = candidate->Area;
    533564
    534565    entry->DeltaEta = candidate->DeltaEta;
Note: See TracChangeset for help on using the changeset viewer.