Fork me on GitHub

Changes in / [659c7b6:f688c89] in git


Ignore:
Files:
11 added
17 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r659c7b6 rf688c89  
    319319        modules/EnergySmearing.h \
    320320        modules/MomentumSmearing.h \
     321        modules/TrackSmearing.h \
    321322        modules/ImpactParameterSmearing.h \
    322323        modules/TimeSmearing.h \
    … …  
    341342        modules/StatusPidFilter.h \
    342343        modules/PdgCodeFilter.h \
     344        modules/BeamSpotFilter.h \
    343345        modules/Cloner.h \
    344346        modules/Weighter.h \
    … …  
    346348        modules/JetFlavorAssociation.h \
    347349        modules/JetFakeParticle.h \
     350        modules/VertexSorter.h \
     351        modules/VertexFinder.h \
    348352        modules/ExampleModule.h
    349353tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    552556        classes/DelphesFactory.h \
    553557        classes/DelphesFormula.h
     558tmp/modules/BeamSpotFilter.$(ObjSuf): \
     559        modules/BeamSpotFilter.$(SrcSuf) \
     560        modules/BeamSpotFilter.h \
     561        classes/DelphesClasses.h \
     562        classes/DelphesFactory.h \
     563        classes/DelphesFormula.h \
     564        external/ExRootAnalysis/ExRootResult.h \
     565        external/ExRootAnalysis/ExRootFilter.h \
     566        external/ExRootAnalysis/ExRootClassifier.h
    554567tmp/modules/Calorimeter.$(ObjSuf): \
    555568        modules/Calorimeter.$(SrcSuf) \
    … …  
    842855        external/ExRootAnalysis/ExRootFilter.h \
    843856        external/ExRootAnalysis/ExRootClassifier.h
     857tmp/modules/TrackSmearing.$(ObjSuf): \
     858        modules/TrackSmearing.$(SrcSuf) \
     859        modules/TrackSmearing.h \
     860        classes/DelphesClasses.h \
     861        classes/DelphesFactory.h \
     862        classes/DelphesFormula.h \
     863        external/ExRootAnalysis/ExRootResult.h \
     864        external/ExRootAnalysis/ExRootFilter.h \
     865        external/ExRootAnalysis/ExRootClassifier.h
    844866tmp/modules/TreeWriter.$(ObjSuf): \
    845867        modules/TreeWriter.$(SrcSuf) \
    … …  
    858880        classes/DelphesFactory.h \
    859881        classes/DelphesFormula.h \
     882        external/ExRootAnalysis/ExRootResult.h \
     883        external/ExRootAnalysis/ExRootFilter.h \
     884        external/ExRootAnalysis/ExRootClassifier.h
     885tmp/modules/VertexFinder.$(ObjSuf): \
     886        modules/VertexFinder.$(SrcSuf) \
     887        modules/VertexFinder.h \
     888        classes/DelphesClasses.h \
     889        classes/DelphesFactory.h \
     890        classes/DelphesFormula.h \
     891        classes/DelphesPileUpReader.h \
     892        external/ExRootAnalysis/ExRootResult.h \
     893        external/ExRootAnalysis/ExRootFilter.h \
     894        external/ExRootAnalysis/ExRootClassifier.h
     895tmp/modules/VertexSorter.$(ObjSuf): \
     896        modules/VertexSorter.$(SrcSuf) \
     897        modules/VertexSorter.h \
     898        classes/DelphesClasses.h \
     899        classes/DelphesFactory.h \
     900        classes/DelphesFormula.h \
     901        classes/DelphesPileUpReader.h \
    860902        external/ExRootAnalysis/ExRootResult.h \
    861903        external/ExRootAnalysis/ExRootFilter.h \
    … …  
    921963        tmp/modules/AngularSmearing.$(ObjSuf) \
    922964        tmp/modules/BTagging.$(ObjSuf) \
     965        tmp/modules/BeamSpotFilter.$(ObjSuf) \
    923966        tmp/modules/Calorimeter.$(ObjSuf) \
    924967        tmp/modules/Cloner.$(ObjSuf) \
    … …  
    952995        tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
    953996        tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
     997        tmp/modules/TrackSmearing.$(ObjSuf) \
    954998        tmp/modules/TreeWriter.$(ObjSuf) \
    955999        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
     1000        tmp/modules/VertexFinder.$(ObjSuf) \
     1001        tmp/modules/VertexSorter.$(ObjSuf) \
    9561002        tmp/modules/Weighter.$(ObjSuf)
    9571003
    … …  
    15581604        tmp/external/tcl/tclVar.$(ObjSuf)
    15591605
     1606modules/TrackSmearing.h: \
     1607        classes/DelphesModule.h
     1608        @touch $@
     1609
    15601610external/fastjet/ClusterSequence.hh: \
    15611611        external/fastjet/PseudoJet.hh \
    … …  
    18151865        @touch $@
    18161866
     1867modules/VertexSorter.h: \
     1868        classes/DelphesModule.h \
     1869        classes/DelphesClasses.h
     1870        @touch $@
     1871
    18171872modules/Delphes.h: \
    18181873        classes/DelphesModule.h
     1874        @touch $@
     1875
     1876modules/VertexFinder.h: \
     1877        classes/DelphesModule.h \
     1878        classes/DelphesClasses.h
    18191879        @touch $@
    18201880
    … …  
    19832043
    19842044modules/FastJetFinder.h: \
     2045        classes/DelphesModule.h
     2046        @touch $@
     2047
     2048modules/BeamSpotFilter.h: \
    19852049        classes/DelphesModule.h
    19862050        @touch $@

  • classes/DelphesClasses.cc

    r659c7b6 rf688c89  
    4141CompBase *Tower::fgCompare = CompE<Tower>::Instance();
    4242CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
     43CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
    4344CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
    4445
    … …  
    121122  Charge(0), Mass(0.0),
    122123  IsPU(0), IsRecoPU(0), IsConstituent(0), IsFromConversion(0),
     124  ClusterIndex(-1), ClusterNDF(0), ClusterSigma(0), SumPT2(0), BTVSumPT2(0), GenDeltaZ(0), GenSumPT2(0),
    123125  Flavor(0), FlavorAlgo(0), FlavorPhys(0),
    124126  BTag(0), BTagAlgo(0), BTagPhys(0),
    … …  
    127129  Momentum(0.0, 0.0, 0.0, 0.0),
    128130  Position(0.0, 0.0, 0.0, 0.0),
     131  PositionError(0.0, 0.0, 0.0, 0.0),
     132  InitialPosition(0.0, 0.0, 0.0, 0.0),
    129133  Area(0.0, 0.0, 0.0, 0.0),
    130   Dxy(0), SDxy(0), Xd(0), Yd(0), Zd(0),
     134  L(0),
     135  D0(0), ErrorD0(0),
     136  DZ(0), ErrorDZ(0),
     137  P(0),  ErrorP(0),
     138  PT(0), ErrorPT(0),
     139  CtgTheta(0), ErrorCtgTheta(0),
     140  Phi(0), ErrorPhi(0), 
     141  Xd(0), Yd(0), Zd(0),
    131142  TrackResolution(0),
    132143  NCharged(0),
    … …  
    245256  object.IsConstituent = IsConstituent;
    246257  object.IsFromConversion = IsFromConversion;
     258  object.ClusterIndex = ClusterIndex;
     259  object.ClusterNDF = ClusterNDF;
     260  object.ClusterSigma = ClusterSigma;
     261  object.SumPT2 = SumPT2;
     262  object.BTVSumPT2 = BTVSumPT2;
     263  object.GenDeltaZ = GenDeltaZ;
     264  object.GenSumPT2 = GenSumPT2;
    247265  object.Flavor = Flavor;
    248266  object.FlavorAlgo = FlavorAlgo;
    … …  
    262280  object.Momentum = Momentum;
    263281  object.Position = Position;
     282  object.InitialPosition = InitialPosition;
     283  object.PositionError = PositionError;
    264284  object.Area = Area;
    265   object.Dxy = Dxy;
    266   object.SDxy = SDxy;
     285  object.L = L;
     286  object.ErrorT = ErrorT;
     287  object.D0 = D0;
     288  object.ErrorD0 = ErrorD0;
     289  object.DZ = DZ;
     290  object.ErrorDZ = ErrorDZ;
     291  object.P = P;
     292  object.ErrorP = ErrorP;
     293  object.PT = PT;
     294  object.ErrorPT = ErrorPT;
     295  object.CtgTheta = CtgTheta ;
     296  object.ErrorCtgTheta = ErrorCtgTheta;
     297  object.Phi = Phi;
     298  object.ErrorPhi = ErrorPhi; 
    267299  object.Xd = Xd;
    268300  object.Yd = Yd;
    … …  
    282314  object.SumPtChargedPU = SumPtChargedPU;
    283315  object.SumPt = SumPt;
    284 
     316  object.ClusterIndex = ClusterIndex;
     317  object.ClusterNDF = ClusterNDF;
     318  object.ClusterSigma = ClusterSigma;
     319  object.SumPT2 = SumPT2;
     320 
    285321  object.FracPt[0] = FracPt[0];
    286322  object.FracPt[1] = FracPt[1];
    … …  
    363399  Momentum.SetXYZT(0.0, 0.0, 0.0, 0.0);
    364400  Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
     401  InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    365402  Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
    366   Dxy = 0.0;
    367   SDxy = 0.0;
     403  L = 0.0;
     404  ErrorT = 0.0;
     405  D0 = 0.0; 
     406  ErrorD0 = 0.0;
     407  DZ = 0.0;
     408  ErrorDZ = 0.0;
     409  P =0.0;
     410  ErrorP =0.0;
     411  PT = 0.0;
     412  ErrorPT = 0.0;
     413  CtgTheta = 0.0;
     414  ErrorCtgTheta = 0.0;
     415  Phi = 0.0;
     416  ErrorPhi = 0.0;
    368417  Xd = 0.0;
    369418  Yd = 0.0;
    … …  
    387436  SumPt = -999;
    388437
     438  ClusterIndex = -1;
     439  ClusterNDF = -99;
     440  ClusterSigma = 0.0;
     441  SumPT2 = 0.0;
     442  BTVSumPT2 = 0.0;
     443  GenDeltaZ = 0.0;
     444  GenSumPT2 = 0.0;
     445 
    389446  FracPt[0] = 0.0;
    390447  FracPt[1] = 0.0;

  • classes/DelphesClasses.h

    r659c7b6 rf688c89  
    147147  Float_t Pz; // particle momentum vector (z component) | hepevt.phep[number][2]
    148148
    149   Float_t PT; // particle transverse momentum
     149  Float_t D0;
     150  Float_t DZ;
     151  Float_t P;
     152  Float_t PT;
     153  Float_t CtgTheta;
     154  Float_t Phi;
    150155  Float_t Eta; // particle pseudorapidity
    151   Float_t Phi; // particle azimuthal angle
    152 
    153156  Float_t Rapidity; // particle rapidity
    154157
    … …  
    168171//---------------------------------------------------------------------------
    169172
    170 class Vertex: public TObject
     173class Vertex: public SortableObject
    171174{
    172175public:
    … …  
    176179  Float_t Z; // vertex position (z component)
    177180
    178   ClassDef(Vertex, 1)
     181  Double_t ErrorX;
     182  Double_t ErrorY;
     183  Double_t ErrorZ;
     184  Double_t ErrorT;
     185
     186  Int_t Index;
     187  Int_t NDF;
     188  Double_t Sigma;
     189  Double_t SumPT2;
     190  Double_t BTVSumPT2;
     191  Double_t GenDeltaZ;
     192  Double_t GenSumPT2;
     193
     194  TRefArray Constituents; // references to constituents
     195
     196  static CompBase *fgCompare; //!
     197  const CompBase *GetCompare() const { return fgCompare; }
     198
     199  ClassDef(Vertex, 3)
    179200};
    180201
    … …  
    397418  Int_t Charge; // track charge
    398419
    399   Float_t PT; // track transverse momentum
    400 
    401420  Float_t Eta; // track pseudorapidity
    402   Float_t Phi; // track azimuthal angle
    403421
    404422  Float_t EtaOuter; // track pseudorapidity at the tracker edge
    … …  
    415433  Float_t TOuter; // track position (z component) at the tracker edge
    416434
    417   Float_t Dxy;     // track signed transverse impact parameter
    418   Float_t SDxy;    // signed error on the track signed transverse impact parameter
     435  Float_t L; // track path length
     436  Float_t ErrorT; // error on the time measurement
     437
     438  Float_t D0;     // track signed transverse impact parameter
     439  Float_t ErrorD0;    // signed error on the track signed transverse impact parameter
     440
     441  Float_t DZ; // track transverse momentum
     442  Float_t ErrorDZ; // track transverse momentum error
     443
     444  Float_t P; // track transverse momentum
     445  Float_t ErrorP; // track transverse momentum error
     446
     447  Float_t PT; // track transverse momentum
     448  Float_t ErrorPT; // track transverse momentum error
     449
     450  Float_t CtgTheta; // track transverse momentum
     451  Float_t ErrorCtgTheta; // track transverse momentum error
     452
     453  Float_t Phi; // track azimuthal angle
     454  Float_t ErrorPhi; // track azimuthal angle
     455
    419456  Float_t Xd;      // X coordinate of point of closest approach to vertex
    420457  Float_t Yd;      // Y coordinate of point of closest approach to vertex
    … …  
    423460  TRef Particle; // reference to generated particle
    424461
     462  Int_t VertexIndex; // reference to vertex
     463
    425464  static CompBase *fgCompare; //!
    426465  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    526565  Float_t DeltaPhi;
    527566
    528   TLorentzVector Momentum, Position, Area;
    529 
    530   Float_t Dxy;
    531   Float_t SDxy;
     567  TLorentzVector Momentum, Position, InitialPosition, PositionError, Area;
     568
     569  Float_t L; // path length
     570  Float_t ErrorT; // path length
     571  Float_t D0;
     572  Float_t ErrorD0;
     573  Float_t DZ;
     574  Float_t ErrorDZ;
     575  Float_t P;
     576  Float_t ErrorP;
     577  Float_t PT;
     578  Float_t ErrorPT;
     579  Float_t CtgTheta;
     580  Float_t ErrorCtgTheta;
     581  Float_t Phi;
     582  Float_t ErrorPhi;
     583
    532584  Float_t Xd;
    533585  Float_t Yd;
    … …  
    535587
    536588  // tracking resolution
    537  
     589
    538590  Float_t TrackResolution;
    539591
    … …  
    562614  Float_t SumPt;
    563615
     616  // vertex variables
     617
     618  Int_t ClusterIndex;
     619  Int_t ClusterNDF;
     620  Double_t ClusterSigma;
     621  Double_t SumPT2;
     622  Double_t BTVSumPT2;
     623  Double_t GenDeltaZ;
     624  Double_t GenSumPT2;
     625
    564626  // N-subjettiness variables
    565627
    … …  
    595657  void SetFactory(DelphesFactory *factory) { fFactory = factory; }
    596658
    597   ClassDef(Candidate, 4)
     659  ClassDef(Candidate, 5)
    598660};
    599661

  • classes/SortableObject.h

    r659c7b6 rf688c89  
    156156      return -1;
    157157    else if(t1->ET < t2->ET)
     158      return 1;
     159    else
     160      return 0;
     161  }
     162};
     163
     164//---------------------------------------------------------------------------
     165
     166template <typename T>
     167class CompSumPT2: public CompBase
     168{
     169  CompSumPT2() {}
     170public:
     171  static CompSumPT2 *Instance()
     172  {
     173    static CompSumPT2 single;
     174    return &single;
     175  }
     176
     177  Int_t Compare(const TObject *obj1, const TObject *obj2) const
     178  {
     179    const T *t1 = static_cast<const T*>(obj1);
     180    const T *t2 = static_cast<const T*>(obj2);
     181    if(t1->SumPT2 > t2->SumPT2)
     182      return -1;
     183    else if(t1->SumPT2 < t2->SumPT2)
    158184      return 1;
    159185    else

  • doc/genMakefile.tcl

    r659c7b6 rf688c89  
    263263executableDeps {converters/*.cpp} {examples/*.cpp}
    264264
    265 executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp}
     265executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp} {readers/DelphesROOT.cpp}
    266266
    267267puts {ifeq ($(HAS_CMSSW),true)}

  • modules/Calorimeter.cc

    r659c7b6 rf688c89  
    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;
    627     tower->PID = 22;
    628 
     559   
    629560    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron tower
     561   
     562    //clone tracks
     563    fItECalTowerTrackArray->Reset();
     564    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     565    {
     566      mother = track;
     567      track = static_cast<Candidate*>(track->Clone());
     568      track->AddCandidate(mother);
     569
     570      fEFlowTrackOutputArray->Add(track);
     571    }
     572 
     573  }
     574 
     575  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     576  else if(fECalTrackEnergy > 0.0)
     577  {
     578    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
     579    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
     580 
     581    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
     582    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     583
     584    //rescale tracks
     585    fItECalTowerTrackArray->Reset();
     586    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     587    { 
     588      mother = track;
     589      track = static_cast<Candidate*>(track->Clone());
     590      track->AddCandidate(mother);
     591
     592      track->Momentum *= rescaleFactor;
     593
     594      fEFlowTrackOutputArray->Add(track);
     595    }
     596  }
     597
     598
     599  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
     600  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
     601  {
     602    // create new photon tower
    634603    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     604    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     605   
     606    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     607    tower->Ehad = hcalNeutralEnergy;
    639608    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     609   
    642610    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     611   
     612    //clone tracks
     613    fItHCalTowerTrackArray->Reset();
     614    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     615    {
     616      mother = track;
     617      track = static_cast<Candidate*>(track->Clone());
     618      track->AddCandidate(mother);
     619
     620      fEFlowTrackOutputArray->Add(track);
     621    }
     622 
     623  }
     624 
     625  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     626  else if(fHCalTrackEnergy > 0.0)
     627  {
     628    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
     629    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
     630 
     631    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     632    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
     633
     634    //rescale tracks
     635    fItECalTowerTrackArray->Reset();
     636    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     637    { 
     638      mother = track;
     639      track = static_cast<Candidate*>(track->Clone());
     640      track->AddCandidate(mother);
     641
     642      track->Momentum *= rescaleFactor;
     643
     644      fEFlowTrackOutputArray->Add(track);
     645    }
     646  }
     647 
     648 
    644649}
    645650

  • modules/Calorimeter.h

    r659c7b6 rf688c89  
    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/ImpactParameterSmearing.cc

    r659c7b6 rf688c89  
    9696{
    9797  Candidate *candidate, *particle, *mother;
    98   Double_t xd, yd, zd, dxy, sx, sy, sz, ddxy;
     98  Double_t xd, yd, zd, d0, sx, sy, sz, dd0;
    9999  Double_t pt, eta, px, py, phi, e;
    100100
    … …  
    103103  {
    104104
    105     // take momentum before smearing (otherwise apply double smearing on dxy)
     105    // take momentum before smearing (otherwise apply double smearing on d0)
    106106    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    107107
    … …  
    131131
    132132    // calculate impact parameter (after-smearing)
    133     dxy = (xd*py - yd*px)/pt;
     133    d0 = (xd*py - yd*px)/pt;
    134134
    135     ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
     135    dd0 = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
    136136
    137137    // fill smeared values in candidate
    … …  
    143143    candidate->Zd = zd;
    144144
    145     candidate->Dxy = dxy;
    146     candidate->SDxy = ddxy;
     145    candidate->D0 = d0;
     146    candidate->ErrorD0 = dd0;
    147147
    148148    candidate->AddCandidate(mother);

  • modules/ModulesLinkDef.h

    r659c7b6 rf688c89  
    3535#include "modules/EnergySmearing.h"
    3636#include "modules/MomentumSmearing.h"
     37#include "modules/TrackSmearing.h"
    3738#include "modules/ImpactParameterSmearing.h"
    3839#include "modules/TimeSmearing.h"
    … …  
    5758#include "modules/StatusPidFilter.h"
    5859#include "modules/PdgCodeFilter.h"
     60#include "modules/BeamSpotFilter.h"
    5961#include "modules/Cloner.h"
    6062#include "modules/Weighter.h"
    … …  
    6264#include "modules/JetFlavorAssociation.h"
    6365#include "modules/JetFakeParticle.h"
     66#include "modules/VertexSorter.h"
     67#include "modules/VertexFinder.h"
     68#include "modules/VertexFinderDA4D.h"
    6469#include "modules/ExampleModule.h"
    6570
    … …  
    7984#pragma link C++ class EnergySmearing+;
    8085#pragma link C++ class MomentumSmearing+;
     86#pragma link C++ class TrackSmearing+;
    8187#pragma link C++ class ImpactParameterSmearing+;
    8288#pragma link C++ class TimeSmearing+;
    … …  
    101107#pragma link C++ class StatusPidFilter+;
    102108#pragma link C++ class PdgCodeFilter+;
     109#pragma link C++ class BeamSpotFilter+;
    103110#pragma link C++ class Cloner+;
    104111#pragma link C++ class Weighter+;
    … …  
    106113#pragma link C++ class JetFlavorAssociation+;
    107114#pragma link C++ class JetFakeParticle+;
     115#pragma link C++ class VertexSorter+;
     116#pragma link C++ class VertexFinder+;
     117#pragma link C++ class VertexFinderDA4D+;
    108118#pragma link C++ class ExampleModule+;
    109119

  • modules/ParticlePropagator.cc

    r659c7b6 rf688c89  
    6767}
    6868
     69
    6970//------------------------------------------------------------------------------
    7071
    … …  
    9192  fItInputArray = fInputArray->MakeIterator();
    9293
     94  // import beamspot
     95  try
     96  {
     97    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
     98  }
     99  catch(runtime_error &e)
     100  {
     101    fBeamSpotInputArray = 0;
     102  } 
    93103  // create output arrays
    94104
    … …  
    111121{
    112122  Candidate *candidate, *mother;
    113   TLorentzVector candidatePosition, candidateMomentum;
     123  TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
    114124  Double_t px, py, pz, pt, pt2, e, q;
    115125  Double_t x, y, z, t, r, phi;
    … …  
    120130  Double_t tmp, discr, discr2;
    121131  Double_t delta, gammam, omega, asinrho;
    122   Double_t rcu, rc2, dxy, xd, yd, zd;
    123 
     132  Double_t rcu, rc2, xd, yd, zd;
     133  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
     134  Double_t bsx, bsy, bsz;
     135         
    124136  const Double_t c_light = 2.99792458E8;
    125 
     137 
     138  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
     139    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
     140  else
     141  {
     142    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
     143    beamSpotPosition = beamSpotCandidate.Position;
     144  }
     145 
    126146  fItInputArray->Reset();
    127147  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    … …  
    132152    y = candidatePosition.Y()*1.0E-3;
    133153    z = candidatePosition.Z()*1.0E-3;
     154   
     155    bsx = beamSpotPosition.X()*1.0E-3;
     156    bsy = beamSpotPosition.Y()*1.0E-3;
     157    bsz = beamSpotPosition.Z()*1.0E-3;
     158   
    134159    q = candidate->Charge;
    135160
    … …  
    181206      y_t = y + py*t;
    182207      z_t = z + pz*t;
     208     
     209      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
    183210
    184211      mother = candidate;
    185212      candidate = static_cast<Candidate*>(candidate->Clone());
    186213
     214      candidate->InitialPosition = candidatePosition;
    187215      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
    188 
     216      candidate->L = l*1.0E3;
     217   
    189218      candidate->Momentum = candidateMomentum;
    190219      candidate->AddCandidate(mother);
    … …  
    238267      yd = (rc2 > 0.0) ? yd / rc2 : -999;
    239268      zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
    240 
    241       // calculate impact paramater
    242       dxy = (xd*py - yd*px)/pt;
    243 
     269     
     270      // use perigee momentum rather than original particle
     271      // momentum, since the orignal particle momentum isn't known
     272     
     273      px = TMath::Sign(1.0,r) * pt * (-y_c / r_c);
     274      py = TMath::Sign(1.0,r) * pt * (x_c / r_c);
     275      etap = candidateMomentum.Eta();
     276      phip = TMath::ATan2(py, px);
     277       
     278      candidateMomentum.SetPtEtaPhiE(pt, etap, phip, candidateMomentum.E());
     279     
     280      // calculate additional track parameters (correct for beamspot position)
     281       
     282      d0        = (  (x - bsx) * py - (y - bsy) * px) / pt;
     283      dz        = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
     284      p         = candidateMomentum.P();
     285      ctgTheta  = 1.0 / TMath::Tan (candidateMomentum.Theta ());
     286         
    244287      // 3. time evaluation t = TMath::Min(t_r, t_z)
    245288      //    t_r : time to exit from the sides
    … …  
    287330      r_t = TMath::Hypot(x_t, y_t);
    288331
     332     
     333      // compute path length for an helix
     334     
     335      alpha = pz*1.0E9 / c_light / gammam;
     336      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
     337         
    289338      if(r_t > 0.0)
    290339      {
     340       
     341        // store these variables before cloning
     342           
     343        candidate->D0 = d0*1.0E3;
     344        candidate->DZ = dz*1.0E3;
     345        candidate->P  = p;
     346        candidate->PT = pt;
     347        candidate->CtgTheta = ctgTheta;
     348        candidate->Phi = phip;
     349     
    291350        mother = candidate;
    292351        candidate = static_cast<Candidate*>(candidate->Clone());
    293352
     353        candidate->InitialPosition = candidatePosition;
    294354        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
    295355
    296356        candidate->Momentum = candidateMomentum;
    297         candidate->Dxy = dxy*1.0E3;
    298         candidate->Xd = xd*1.0E3;
     357       
     358            candidate->L  =  l*1.0E3;
     359       
     360            candidate->Xd = xd*1.0E3;
    299361        candidate->Yd = yd*1.0E3;
    300362        candidate->Zd = zd*1.0E3;

  • modules/ParticlePropagator.h

    r659c7b6 rf688c89  
    3535class TClonesArray;
    3636class TIterator;
     37class TLorentzVector;
    3738
    3839class ParticlePropagator: public DelphesModule
    … …  
    5556
    5657  const TObjArray *fInputArray; //!
     58  const TObjArray *fBeamSpotInputArray; //!
    5759
    5860  TObjArray *fOutputArray; //!

  • modules/PileUpMerger.cc

    r659c7b6 rf688c89  
    115115  TDatabasePDG *pdg = TDatabasePDG::Instance();
    116116  TParticlePDG *pdgParticle;
    117   Int_t pid;
     117  Int_t pid, nch, nvtx = -1;
    118118  Float_t x, y, z, t, vx, vy;
    119   Float_t px, py, pz, e;
    120   Double_t dz, dphi, dt;
     119  Float_t px, py, pz, e, pt;
     120  Double_t dz, dphi, dt, sumpt2, dz0, dt0;
    121121  Int_t numberOfEvents, event, numberOfParticles;
    122122  Long64_t allEntries, entry;
    … …  
    132132  fFunction->GetRandom2(dz, dt);
    133133
     134  dz0 = -1.0e6;
     135  dt0 = -1.0e6;
     136
    134137  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    135138  dz *= 1.0E3; // necessary in order to make z in mm
     139
     140  //cout<<dz<<","<<dt<<endl;
     141
    136142  vx = 0.0;
    137143  vy = 0.0;
     144
    138145  numberOfParticles = fInputArray->GetEntriesFast();
     146  nch = 0;
     147  sumpt2 = 0.0;
     148
     149  factory = GetFactory();
     150  vertex = factory->NewCandidate();
     151
    139152  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    140153  {
    … …  
    143156    z = candidate->Position.Z();
    144157    t = candidate->Position.T();
    145     candidate->Position.SetZ(z + dz);
    146     candidate->Position.SetT(t + dt);
     158    pt = candidate->Momentum.Pt();
     159
     160    // take postion and time from first stable particle
     161    if (dz0 < -999999.0)
     162      dz0 = z;
     163    if (dt0 < -999999.0)
     164      dt0 = t;
     165
     166    // cancel any possible offset in position and time the input file
     167    candidate->Position.SetZ(z - dz0 + dz);
     168    candidate->Position.SetT(t - dt0 + dt);
     169
     170    candidate->IsPU = 0;
     171
    147172    fParticleOutputArray->Add(candidate);
     173
     174    if(TMath::Abs(candidate->Charge) >  1.0E-9)
     175    {
     176      nch++;
     177      sumpt2 += pt*pt;
     178      vertex->AddCandidate(candidate);
     179    }
    148180  }
    149181
    150182  if(numberOfParticles > 0)
    151183  {
    152     vx /= numberOfParticles;
    153     vy /= numberOfParticles;
     184    vx /= sumpt2;
     185    vy /= sumpt2;
    154186  }
    155187
    156   factory = GetFactory();
    157 
    158   vertex = factory->NewCandidate();
     188  nvtx++;
    159189  vertex->Position.SetXYZT(vx, vy, dz, dt);
     190  vertex->ClusterIndex = nvtx;
     191  vertex->ClusterNDF = nch;
     192  vertex->SumPT2 = sumpt2;
     193  vertex->GenSumPT2 = sumpt2;
    160194  fVertexOutputArray->Add(vertex);
    161195
    … …  
    170204      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
    171205      break;
     206    case 2:
     207      numberOfEvents = fMeanPileUp;
     208      break;
    172209    default:
    173210      numberOfEvents = gRandom->Poisson(fMeanPileUp);
    … …  
    176213
    177214  allEntries = fReader->GetEntries();
     215
    178216
    179217  for(event = 0; event < numberOfEvents; ++event)
    … …  
    198236    vx = 0.0;
    199237    vy = 0.0;
     238
    200239    numberOfParticles = 0;
     240    sumpt2 = 0.0;
     241
     242    //factory = GetFactory();
     243    vertex = factory->NewCandidate();
     244
    201245    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    202246    {
    … …  
    215259      candidate->Momentum.SetPxPyPzE(px, py, pz, e);
    216260      candidate->Momentum.RotateZ(dphi);
     261      pt = candidate->Momentum.Pt();
    217262
    218263      x -= fInputBeamSpotX;
    … …  
    224269      vx += candidate->Position.X();
    225270      vy += candidate->Position.Y();
     271
    226272      ++numberOfParticles;
     273      if(TMath::Abs(candidate->Charge) >  1.0E-9)
     274      {
     275        nch++;
     276        sumpt2 += pt*pt;
     277        vertex->AddCandidate(candidate);
     278      }
    227279
    228280      fParticleOutputArray->Add(candidate);
    … …  
    235287    }
    236288
    237     vertex = factory->NewCandidate();
     289    nvtx++;
     290
    238291    vertex->Position.SetXYZT(vx, vy, dz, dt);
     292
     293    vertex->ClusterIndex = nvtx;
     294    vertex->ClusterNDF = nch;
     295    vertex->SumPT2 = sumpt2;
     296    vertex->GenSumPT2 = sumpt2;
     297
    239298    vertex->IsPU = 1;
    240299
    241300    fVertexOutputArray->Add(vertex);
     301
    242302  }
    243303}

  • modules/SimpleCalorimeter.cc

    r659c7b6 rf688c89  
    366366      energy = momentum.E() * fTrackFractions[number];
    367367
    368       fTrackTime += TMath::Sqrt(energy)*position.T();
    369       fTrackTimeWeight += TMath::Sqrt(energy);
     368      fTrackTime += energy*position.T();
     369      fTrackTimeWeight += energy;
    370370
    371371      if(fTrackFractions[number] > 1.0E-9)
    … …  
    403403    fTowerEnergy += energy;
    404404
    405     fTowerTime += TMath::Sqrt(energy)*position.T();
    406     fTowerTimeWeight += TMath::Sqrt(energy);
     405    fTowerTime += energy*position.T();
     406    fTowerTimeWeight += energy;
    407407
    408408    fTower->AddCandidate(particle);

  • modules/TimeSmearing.cc

    r659c7b6 rf688c89  
    9393{
    9494  Candidate *candidate, *mother;
    95   Double_t t;
     95  Double_t ti, tf_smeared, tf;
    9696  const Double_t c_light = 2.99792458E8;
    9797
    … …  
    9999  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    100100  {
    101     const TLorentzVector &candidatePosition = candidate->Position;
    102     t = candidatePosition.T()*1.0E-3/c_light;
     101    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
     102    const TLorentzVector &candidateFinalPosition = candidate->Position;
     103
     104    ti = candidateInitialPosition.T()*1.0E-3/c_light;
     105    tf = candidateFinalPosition.T()*1.0E-3/c_light;
    103106
    104107    // apply smearing formula
    105     t = gRandom->Gaus(t, fTimeResolution);
     108    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
     109    ti = ti + tf_smeared - tf;
    106110
    107111    mother = candidate;
    108112    candidate = static_cast<Candidate*>(candidate->Clone());
    109     candidate->Position.SetT(t*1.0E3*c_light);
     113    candidate->InitialPosition.SetT(ti*1.0E3*c_light);
     114    candidate->Position.SetT(tf*1.0E3*c_light);
     115
     116    candidate->ErrorT = fTimeResolution*1.0E3*c_light;
    110117
    111118    candidate->AddCandidate(mother);

  • modules/TrackCountingBTagging.cc

    r659c7b6 rf688c89  
    9696
    9797  Double_t jpx, jpy;
    98   Double_t dr, tpx, tpy, tpt;
    99   Double_t xd, yd, dxy, ddxy, ip, sip;
     98  Double_t dr, tpt;
     99  Double_t xd, yd, d0, dd0, ip, sip;
    100100
    101101  Int_t sign;
    … …  
    117117    {
    118118      const TLorentzVector &trkMomentum = track->Momentum;
    119 
     119     
    120120      dr = jetMomentum.DeltaR(trkMomentum);
    121 
    122121      tpt = trkMomentum.Pt();
    123       tpx = trkMomentum.Px();
    124       tpy = trkMomentum.Py();
    125 
    126122      xd = track->Xd;
    127123      yd = track->Yd;
    128       dxy = TMath::Hypot(xd, yd);
    129       ddxy = track->SDxy;
     124      d0 = TMath::Hypot(xd, yd);
     125      dd0 = track->ErrorD0;
    130126
    131127      if(tpt < fPtMin) continue;
    132128      if(dr > fDeltaR) continue;
    133       if(dxy > fIPmax) continue;
     129      if(d0 > fIPmax) continue;
    134130
    135131      sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
    136132
    137       ip = sign*dxy;
    138       sip = ip / TMath::Abs(ddxy);
     133      ip = sign*d0;
     134      sip = ip / TMath::Abs(dd0);
    139135
    140136      if(sip > fSigMin) count++;

  • modules/TreeWriter.cc

    r659c7b6 rf688c89  
    215215    entry->Pz = momentum.Pz();
    216216
     217    entry->D0            = candidate->D0;
     218    entry->DZ            = candidate->DZ;
     219    entry->P             = candidate->P;
     220    entry->PT            = candidate->PT;
     221    entry->CtgTheta      = candidate->CtgTheta;
     222    entry->Phi           = candidate->Phi;
     223
    217224    entry->Eta = eta;
    218225    entry->Phi = momentum.Phi();
    … …  
    233240{
    234241  TIter iterator(array);
    235   Candidate *candidate = 0;
     242  Candidate *candidate = 0, *constituent = 0;
    236243  Vertex *entry = 0;
    237244
    238245  const Double_t c_light = 2.99792458E8;
    239246
     247  Double_t x, y, z, t, xError, yError, zError, tError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
     248  UInt_t index, ndf;
     249
     250  CompBase *compare = Candidate::fgCompare;
     251  Candidate::fgCompare = CompSumPT2<Candidate>::Instance();
     252  array->Sort();
     253  Candidate::fgCompare = compare;
     254
    240255  // loop over all vertices
    241256  iterator.Reset();
    242257  while((candidate = static_cast<Candidate*>(iterator.Next())))
    243258  {
    244     const TLorentzVector &position = candidate->Position;
     259
     260    index = candidate->ClusterIndex;
     261    ndf = candidate->ClusterNDF;
     262    sigma = candidate->ClusterSigma;
     263    sumPT2 = candidate->SumPT2;
     264    btvSumPT2 = candidate->BTVSumPT2;
     265    genDeltaZ = candidate->GenDeltaZ;
     266    genSumPT2 = candidate->GenSumPT2;
     267
     268    x = candidate->Position.X();
     269    y = candidate->Position.Y();
     270    z = candidate->Position.Z();
     271    t = candidate->Position.T()*1.0E-3/c_light;
     272
     273    xError = candidate->PositionError.X ();
     274    yError = candidate->PositionError.Y ();
     275    zError = candidate->PositionError.Z ();
     276    tError = candidate->PositionError.T ()*1.0E-3/c_light;
    245277
    246278    entry = static_cast<Vertex*>(branch->NewEntry());
    247279
    248     entry->X = position.X();
    249     entry->Y = position.Y();
    250     entry->Z = position.Z();
    251     entry->T = position.T()*1.0E-3/c_light;
    252   }
    253 }
     280    entry->Index = index;
     281    entry->NDF = ndf;
     282    entry->Sigma = sigma;
     283    entry->SumPT2 = sumPT2;
     284    entry->BTVSumPT2 = btvSumPT2;
     285    entry->GenDeltaZ = genDeltaZ;
     286    entry->GenSumPT2 = genSumPT2;
     287
     288    entry->X = x;
     289    entry->Y = y;
     290    entry->Z = z;
     291    entry->T = t;
     292
     293    entry->ErrorX = xError;
     294    entry->ErrorY = yError;
     295    entry->ErrorZ = zError;
     296    entry->ErrorT = tError;
     297
     298
     299    TIter itConstituents(candidate->GetCandidates());
     300    itConstituents.Reset();
     301    entry->Constituents.Clear();
     302    while((constituent = static_cast<Candidate*>(itConstituents.Next())))
     303    {
     304      entry->Constituents.Add(constituent);
     305    }
     306
     307  }
     308}
     309
    254310
    255311//------------------------------------------------------------------------------
    … …  
    292348    entry->TOuter = position.T()*1.0E-3/c_light;
    293349
    294     entry->Dxy = candidate->Dxy;
    295     entry->SDxy = candidate->SDxy ;
     350    entry->L = candidate->L;
     351
     352    entry->D0            = candidate->D0;
     353    entry->ErrorD0       = candidate->ErrorD0;
     354    entry->DZ            = candidate->DZ;
     355    entry->ErrorDZ       = candidate->ErrorDZ;
     356    entry->P             = candidate->P;
     357    entry->ErrorP        = candidate->ErrorP;
     358    entry->PT            = candidate->PT;
     359    entry->ErrorPT       = candidate->ErrorPT;
     360    entry->CtgTheta      = candidate->CtgTheta;
     361    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
     362    entry->Phi           = candidate->Phi;
     363    entry->ErrorPhi      = candidate->ErrorPhi;
     364
    296365    entry->Xd = candidate->Xd;
    297366    entry->Yd = candidate->Yd;
    … …  
    307376
    308377    entry->Eta = eta;
    309     entry->Phi = momentum.Phi();
    310     entry->PT = pt;
    311378
    312379    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    … …  
    319386
    320387    entry->Particle = particle;
     388
     389    entry->VertexIndex = candidate->ClusterIndex;
     390
    321391  }
    322392}

  • readers/DelphesCMSFWLite.cpp

    r659c7b6 rf688c89  
    240240
    241241    branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
    242     branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
     242    branchRwgt = treeWriter->NewBranch("Weight", Weight::Class());
    243243
    244244    confReader = new ExRootConfReader;
Note: See TracChangeset for help on using the changeset viewer.