Fork me on GitHub

Changeset b9ae4c3 in git for modules


Ignore:
Timestamp:
Sep 2, 2016, 3:46:14 PM (8 years ago)
Author:
GitHub <noreply@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
fa068d3
Parents:
ec5e04b (diff), 23389ff (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Pavel Demin <pavel-demin@…> (09/02/16 15:46:14)
git-committer:
GitHub <noreply@…> (09/02/16 15:46:14)
Message:

Merge pull request #32 from delphes/dev_01

Dev 01

Location:
modules
Files:
10 added
14 edited

Legend:

Unmodified
Added
Removed

  • modules/Calorimeter.cc

    rec5e04b rb9ae4c3  
    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
    … …  
    219213  Double_t ecalEnergy, hcalEnergy;
    220214  Double_t ecalSigma, hcalSigma;
     215  Double_t energyGuess;
    221216  Int_t pdgCode;
    222217
    … …  
    368363      fHCalTowerEnergy = 0.0;
    369364
    370       fECalTrackEnergy[0] = 0.0;
    371       fECalTrackEnergy[1] = 0.0;
    372 
    373       fHCalTrackEnergy[0] = 0.0;
    374       fHCalTrackEnergy[1] = 0.0;
    375 
     365      fECalTrackEnergy = 0.0;
     366      fHCalTrackEnergy = 0.0;
     367     
     368      fECalTrackSigma = 0.0;
     369      fHCalTrackSigma = 0.0;
     370     
    376371      fTowerTrackHits = 0;
    377372      fTowerPhotonHits = 0;
    378373
    379       fECalTowerTrackArray[0]->Clear();
    380       fECalTowerTrackArray[1]->Clear();
    381 
    382       fHCalTowerTrackArray[0]->Clear();
    383       fHCalTowerTrackArray[1]->Clear();
     374      fECalTowerTrackArray->Clear();
     375      fHCalTowerTrackArray->Clear();
     376   
    384377    }
    385378
    … …  
    406399      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    407400      {
    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         }
     401        fECalTrackEnergy += ecalEnergy;
     402        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());       
     403        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;       
     404        else energyGuess = momentum.E();
     405
     406        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
     407        fECalTowerTrackArray->Add(track);
    419408      }
     409     
    420410      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    421411      {
     412        fHCalTrackEnergy += hcalEnergy;
    422413        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         }
     414        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
     415        else energyGuess = momentum.E();
     416
     417        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
     418        fHCalTowerTrackArray->Add(track);
    433419      }
     420     
    434421      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    435422      {
    … …  
    476463  Double_t energy, pt, eta, phi;
    477464  Double_t ecalEnergy, hcalEnergy;
     465  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
     466 
    478467  Double_t ecalSigma, hcalSigma;
    479 
     468  Double_t ecalNeutralSigma, hcalNeutralSigma;
     469
     470  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
     471 
    480472  TLorentzVector momentum;
    481473  TFractionMap::iterator itFractionMap;
    … …  
    484476
    485477  if(!fTower) return;
    486 
    487478
    488479  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    554545    fTowerOutputArray->Add(fTower);
    555546  }
    556 
     547 
    557548  // 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)
     549  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
     550  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
     551
     552  //compute neutral excesses
     553  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
     554  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
     555 
     556  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
     557  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
     558 
     559   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
     560  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    618561  {
    619562    // create new photon tower
    620563    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;
     564    pt =  ecalNeutralEnergy / TMath::CosH(eta);
     565   
     566    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     567    tower->Eem = ecalNeutralEnergy;
    626568    tower->Ehad = 0.0;
    627569    tower->PID = 22;
    628 
     570   
    629571    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron tower
     572   
     573    //clone tracks
     574    fItECalTowerTrackArray->Reset();
     575    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     576    {
     577      mother = track;
     578      track = static_cast<Candidate*>(track->Clone());
     579      track->AddCandidate(mother);
     580
     581      fEFlowTrackOutputArray->Add(track);
     582    }
     583 
     584  }
     585 
     586  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     587  else if(fECalTrackEnergy > 0.0)
     588  {
     589    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
     590    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
     591 
     592    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
     593    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     594
     595    //rescale tracks
     596    fItECalTowerTrackArray->Reset();
     597    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     598    { 
     599      mother = track;
     600      track = static_cast<Candidate*>(track->Clone());
     601      track->AddCandidate(mother);
     602
     603      track->Momentum *= rescaleFactor;
     604
     605      fEFlowTrackOutputArray->Add(track);
     606    }
     607  }
     608
     609
     610  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
     611  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
     612  {
     613    // create new photon tower
    634614    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     615    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     616   
     617    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     618    tower->Ehad = hcalNeutralEnergy;
    639619    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     620   
    642621    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     622   
     623    //clone tracks
     624    fItHCalTowerTrackArray->Reset();
     625    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     626    {
     627      mother = track;
     628      track = static_cast<Candidate*>(track->Clone());
     629      track->AddCandidate(mother);
     630
     631      fEFlowTrackOutputArray->Add(track);
     632    }
     633 
     634  }
     635 
     636  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     637  else if(fHCalTrackEnergy > 0.0)
     638  {
     639    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
     640    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
     641 
     642    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     643    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
     644
     645    //rescale tracks
     646    fItHCalTowerTrackArray->Reset();
     647    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     648    { 
     649      mother = track;
     650      track = static_cast<Candidate*>(track->Clone());
     651      track->AddCandidate(mother);
     652
     653      track->Momentum *= rescaleFactor;
     654
     655      fEFlowTrackOutputArray->Add(track);
     656    }
     657  }
     658 
     659 
    644660}
    645661

  • modules/Calorimeter.h

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    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"
    … …  
    5859#include "modules/StatusPidFilter.h"
    5960#include "modules/PdgCodeFilter.h"
     61#include "modules/BeamSpotFilter.h"
    6062#include "modules/RecoPuFilter.h"
    6163#include "modules/Cloner.h"
    … …  
    6466#include "modules/JetFlavorAssociation.h"
    6567#include "modules/JetFakeParticle.h"
     68#include "modules/VertexSorter.h"
     69#include "modules/VertexFinder.h"
     70#include "modules/VertexFinderDA4D.h"
    6671#include "modules/ExampleModule.h"
    6772
    … …  
    8186#pragma link C++ class EnergySmearing+;
    8287#pragma link C++ class MomentumSmearing+;
     88#pragma link C++ class TrackSmearing+;
    8389#pragma link C++ class ImpactParameterSmearing+;
    8490#pragma link C++ class TimeSmearing+;
    … …  
    104110#pragma link C++ class StatusPidFilter+;
    105111#pragma link C++ class PdgCodeFilter+;
     112#pragma link C++ class BeamSpotFilter+;
    106113#pragma link C++ class RecoPuFilter+;
    107114#pragma link C++ class Cloner+;
    … …  
    110117#pragma link C++ class JetFlavorAssociation+;
    111118#pragma link C++ class JetFakeParticle+;
     119#pragma link C++ class VertexSorter+;
     120#pragma link C++ class VertexFinder+;
     121#pragma link C++ class VertexFinderDA4D+;
    112122#pragma link C++ class ExampleModule+;
    113123

  • modules/MomentumSmearing.cc

    rec5e04b rb9ae4c3  
    9696{
    9797  Candidate *candidate, *mother;
    98   Double_t pt, eta, phi, e;
     98  Double_t pt, eta, phi, e, res;
    9999
    100100  fItInputArray->Reset();
    … …  
    107107    pt = candidateMomentum.Pt();
    108108    e = candidateMomentum.E();
     109    res = fFormula->Eval(pt, eta, phi, e);
     110 
     111    // apply smearing formula
     112    //pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     113   
     114    res = ( res > 1.0 ) ? 1.0 : res;
    109115
    110     // apply smearing formula
    111     pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     116    pt = LogNormal(pt, res * pt );
    112117   
    113     if(pt <= 0.0) continue;
     118    //if(pt <= 0.0) continue;
    114119
    115120    mother = candidate;
    … …  
    118123    phi = candidateMomentum.Phi();
    119124    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt*TMath::CosH(eta));
    120     candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     125    //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     126    candidate->TrackResolution = res;
    121127    candidate->AddCandidate(mother);
    122128       
    … …  
    124130  }
    125131}
     132//----------------------------------------------------------------
     133
     134Double_t MomentumSmearing::LogNormal(Double_t mean, Double_t sigma)
     135{
     136  Double_t a, b;
     137
     138  if(mean > 0.0)
     139  {
     140    b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
     141    a = TMath::Log(mean) - 0.5*b*b;
     142
     143    return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
     144  }
     145  else
     146  {
     147    return 0.0;
     148  }
     149}
     150
    126151
    127152//------------------------------------------------------------------------------

  • modules/MomentumSmearing.h

    rec5e04b rb9ae4c3  
    4747private:
    4848
     49  Double_t LogNormal(Double_t mean, Double_t sigma);
     50
    4951  DelphesFormula *fFormula; //!
    5052

  • modules/ParticlePropagator.cc

    rec5e04b rb9ae4c3  
    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;
     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;
    123135
    124136  const Double_t c_light = 2.99792458E8;
     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  }
    125145
    126146  fItInputArray->Reset();
    … …  
    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
    … …  
    182207      z_t = z + pz*t;
    183208
     209      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
     210
    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);
     216      candidate->L = l*1.0E3;
    188217
    189218      candidate->Momentum = candidateMomentum;
    … …  
    239268      zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
    240269
    241       // calculate impact paramater
    242       dxy = (xd*py - yd*px)/pt;
     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
    243287
    244288      // 3. time evaluation t = TMath::Min(t_r, t_z)
    … …  
    287331      r_t = TMath::Hypot(x_t, y_t);
    288332
     333
     334      // compute path length for an helix
     335
     336      alpha = pz*1.0E9 / c_light / gammam;
     337      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
     338
    289339      if(r_t > 0.0)
    290340      {
     341
     342        // store these variables before cloning
     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

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fResolutionFormula = new DelphesFormula;
    63 
    64   for(i = 0; i < 2; ++i)
    65   {
    66     fTowerTrackArray[i] = new TObjArray;
    67     fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
    68   }
     62  fTowerTrackArray = new TObjArray;
     63  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
     64 
    6965}
    7066
    … …  
    7369SimpleCalorimeter::~SimpleCalorimeter()
    7470{
    75   Int_t i;
    76 
     71 
    7772  if(fResolutionFormula) delete fResolutionFormula;
    78 
    79   for(i = 0; i < 2; ++i)
    80   {
    81     if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
    82     if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
    83   }
     73  if(fTowerTrackArray) delete fTowerTrackArray;
     74  if(fItTowerTrackArray) delete fItTowerTrackArray;
     75 
    8476}
    8577
    … …  
    199191  Double_t energy;
    200192  Double_t sigma;
     193  Double_t energyGuess;
     194
    201195  Int_t pdgCode;
    202196
    … …  
    340334      fTowerEnergy = 0.0;
    341335
    342       fTrackEnergy[0] = 0.0;
    343       fTrackEnergy[1] = 0.0;
     336      fTrackEnergy = 0.0;
     337      fTrackSigma = 0.0;
    344338
    345339      fTowerTime = 0.0;
    … …  
    351345      fTowerPhotonHits = 0;
    352346
    353       fTowerTrackArray[0]->Clear();
    354       fTowerTrackArray[1]->Clear();
    355     }
     347      fTowerTrackArray->Clear();
     348     }
    356349
    357350    // check for track hits
    … …  
    371364      if(fTrackFractions[number] > 1.0E-9)
    372365      {
    373         sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    374         if(sigma/momentum.E() < track->TrackResolution)
    375         {
    376           fTrackEnergy[0] += energy;
    377           fTowerTrackArray[0]->Add(track);
    378         }
    379         else
    380         {
    381           fTrackEnergy[1] += energy;
    382           fTowerTrackArray[1]->Add(track);
    383         }
     366             
     367       // compute total charged energy   
     368       fTrackEnergy += energy;
     369       sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
     370       if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy;
     371       else energyGuess = momentum.E();
     372
     373       fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess);
     374       fTowerTrackArray->Add(track);
    384375      }
     376       
    385377      else
    386378      {
    … …  
    403395    fTowerEnergy += energy;
    404396
    405     fTowerTime += TMath::Sqrt(energy)*position.T();
    406     fTowerTimeWeight += TMath::Sqrt(energy);
     397    fTowerTime += energy*position.T();
     398    fTowerTimeWeight += energy;
    407399
    408400    fTower->AddCandidate(particle);
    … …  
    418410{
    419411  Candidate *tower, *track, *mother;
    420   Double_t energy, pt, eta, phi;
    421   Double_t sigma;
     412  Double_t energy,neutralEnergy, pt, eta, phi;
     413  Double_t sigma, neutralSigma;
    422414  Double_t time;
     415   
     416  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    423417
    424418  TLorentzVector momentum;
    … …  
    436430
    437431  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
     432
    438433
    439434  if(fSmearTowerCenter)
    … …  
    464459  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    465460
    466   // fill e-flow candidates
    467 
    468   energy -= fTrackEnergy[1];
    469 
    470   fItTowerTrackArray[0]->Reset();
    471   while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
    472   {
    473     mother = track;
    474     track = static_cast<Candidate*>(track->Clone());
    475     track->AddCandidate(mother);
    476 
    477     track->Momentum *= energy/fTrackEnergy[0];
    478 
    479     fEFlowTrackOutputArray->Add(track);
    480   }
    481 
    482   fItTowerTrackArray[1]->Reset();
    483   while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
    484   {
    485     mother = track;
    486     track = static_cast<Candidate*>(track->Clone());
    487     track->AddCandidate(mother);
    488 
    489     fEFlowTrackOutputArray->Add(track);
    490   }
    491 
    492   if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
    493 
    494   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
    495   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    496 
    497   // save energy excess as an energy flow tower
    498   if(energy > 0.0)
     461 
     462  // e-flow candidates
     463
     464  //compute neutral excess
     465 
     466  fTrackSigma = TMath::Sqrt(fTrackSigma);
     467  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
     468 
     469  //compute sigma_trk total
     470  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
     471   
     472  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
     473  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
    499474  {
    500475    // create new photon tower
    501476    tower = static_cast<Candidate*>(fTower->Clone());
    502     pt = energy / TMath::CosH(eta);
    503 
    504     tower->Eem = (!fIsEcal) ? 0 : energy;
    505     tower->Ehad = (fIsEcal) ? 0 : energy;
     477    pt = neutralEnergy / TMath::CosH(eta);
     478
     479    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
     480    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     481    tower->PID = (fIsEcal) ? 22 : 0;
     482 
     483    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     484    fEFlowTowerOutputArray->Add(tower);
    506485   
    507     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     486    fItTowerTrackArray->Reset();
     487    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     488    {
     489      mother = track;
     490      track = static_cast<Candidate*>(track->Clone());
     491      track->AddCandidate(mother);
     492
     493      fEFlowTrackOutputArray->Add(track);
     494    }
     495  }
    508496   
    509     tower->PID = (fIsEcal) ? 22 : 0;
    510    
    511     fEFlowTowerOutputArray->Add(tower);
    512   }
    513 }
     497  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     498  else if (fTrackEnergy > 0.0)
     499  {
     500    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
     501    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
     502 
     503    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
     504    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
     505   
     506    fItTowerTrackArray->Reset();
     507    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     508    {
     509      mother = track;
     510      track = static_cast<Candidate*>(track->Clone());
     511      track->AddCandidate(mother);
     512
     513      track->Momentum *= rescaleFactor;
     514
     515      fEFlowTrackOutputArray->Add(track);
     516    }
     517  }
     518   
     519 }
    514520
    515521//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

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

  • modules/TimeSmearing.cc

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    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

    rec5e04b rb9ae4c3  
    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//------------------------------------------------------------------------------
    … …  
    261317  Candidate *particle = 0;
    262318  Track *entry = 0;
    263   Double_t pt, signz, cosTheta, eta, rapidity;
     319  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
    264320  const Double_t c_light = 2.99792458E8;
    265321
    … …  
    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
     357    entry->ErrorP        = candidate->ErrorP;
     358    entry->ErrorPT       = candidate->ErrorPT;
     359    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
     360    entry->ErrorPhi      = candidate->ErrorPhi;
     361
    296362    entry->Xd = candidate->Xd;
    297363    entry->Yd = candidate->Yd;
    … …  
    301367
    302368    pt = momentum.Pt();
     369    p = momentum.P();
     370    phi = momentum.Phi();
     371    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1/TMath::Tan(momentum.Theta()) : 1e10;
     372
    303373    cosTheta = TMath::Abs(momentum.CosTheta());
    304374    signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
    … …  
    306376    rapidity = (cosTheta == 1.0 ? signz*999.9 : momentum.Rapidity());
    307377
     378    entry->PT  = pt;
    308379    entry->Eta = eta;
    309     entry->Phi = momentum.Phi();
    310     entry->PT = pt;
     380    entry->Phi = phi;
     381    entry->CtgTheta = ctgTheta;
    311382
    312383    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    … …  
    319390
    320391    entry->Particle = particle;
     392
     393    entry->VertexIndex = candidate->ClusterIndex;
     394
    321395  }
    322396}
Note: See TracChangeset for help on using the changeset viewer.