Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/TrackSmearing.cc

    r341014c rcde9f31  
    1313#include "classes/DelphesFormula.h"
    1414
     15#include "ExRootAnalysis/ExRootResult.h"
     16#include "ExRootAnalysis/ExRootFilter.h"
    1517#include "ExRootAnalysis/ExRootClassifier.h"
    16 #include "ExRootAnalysis/ExRootFilter.h"
    17 #include "ExRootAnalysis/ExRootResult.h"
    18 
     18
     19#include "TMath.h"
     20#include "TString.h"
     21#include "TFormula.h"
     22#include "TRandom3.h"
     23#include "TObjArray.h"
    1924#include "TDatabasePDG.h"
     25#include "TLorentzVector.h"
    2026#include "TFile.h"
    21 #include "TFormula.h"
    22 #include "TLorentzVector.h"
    23 #include "TMath.h"
    24 #include "TObjArray.h"
    2527#include "TProfile2D.h"
    26 #include "TRandom3.h"
    27 #include "TString.h"
    28 
    29 #include <algorithm>
     28
     29#include <algorithm>
     30#include <stdexcept>
    3031#include <iostream>
    3132#include <sstream>
    32 #include <stdexcept>
    3333
    3434using namespace std;
    … …  
    6666
    6767  // !!! IF WE WANT TO KEEP ROOT INPUT !!!
    68   if(string(GetString("D0ResolutionFormula", "0.0")) != "0.0")
    69   {
    70     fD0Formula->Compile(GetString("D0ResolutionFormula", "0.0"));
    71     fUseD0Formula = true;
    72   }
    73   else
    74   {
    75     fD0ResolutionFile = GetString("D0ResolutionFile", "errors.root");
    76     fD0ResolutionHist = GetString("D0ResolutionHist", "d0");
    77     fUseD0Formula = false;
    78   }
    79   if(string(GetString("DZResolutionFormula", "0.0")) != "0.0")
    80   {
    81     fDZFormula->Compile(GetString("DZResolutionFormula", "0.0"));
    82     fUseDZFormula = true;
    83   }
    84   else
    85   {
    86     fDZResolutionFile = GetString("DZResolutionFile", "errors.root");
    87     fDZResolutionHist = GetString("DZResolutionHist", "dz");
    88     fUseDZFormula = false;
    89   }
    90   if(string(GetString("PResolutionFormula", "0.0")) != "0.0")
    91   {
    92     fPFormula->Compile(GetString("PResolutionFormula", "0.0"));
    93     fUsePFormula = true;
    94   }
    95   else
    96   {
    97     fPResolutionFile = GetString("PResolutionFile", "errors.root");
    98     fPResolutionHist = GetString("PResolutionHist", "p");
    99     fUsePFormula = false;
    100   }
    101   if(string(GetString("CtgThetaResolutionFormula", "0.0")) != "0.0")
    102   {
    103     fCtgThetaFormula->Compile(GetString("CtgThetaResolutionFormula", "0.0"));
    104     fUseCtgThetaFormula = true;
    105   }
    106   else
    107   {
    108     fCtgThetaResolutionFile = GetString("CtgThetaResolutionFile", "errors.root");
    109     fCtgThetaResolutionHist = GetString("CtgThetaResolutionHist", "ctgTheta");
    110     fUseCtgThetaFormula = false;
    111   }
    112   if(string(GetString("PhiResolutionFormula", "0.0")) != "0.0")
    113   {
    114     fPhiFormula->Compile(GetString("PhiResolutionFormula", "0.0"));
    115     fUsePhiFormula = true;
    116   }
    117   else
    118   {
    119     fPhiResolutionFile = GetString("PhiResolutionFile", "errors.root");
    120     fPhiResolutionHist = GetString("PhiResolutionHist", "phi");
    121     fUsePhiFormula = false;
    122   }
     68  if (string (GetString("D0ResolutionFormula", "0.0")) != "0.0")
     69    {
     70      fD0Formula->Compile(GetString("D0ResolutionFormula", "0.0"));
     71      fUseD0Formula = true;
     72    }
     73  else
     74    {
     75      fD0ResolutionFile = GetString("D0ResolutionFile", "errors.root");
     76      fD0ResolutionHist = GetString("D0ResolutionHist", "d0");
     77      fUseD0Formula = false;
     78    }
     79  if (string (GetString("DZResolutionFormula", "0.0")) != "0.0")
     80    {
     81      fDZFormula->Compile(GetString("DZResolutionFormula", "0.0"));
     82      fUseDZFormula = true;
     83    }
     84  else
     85    {
     86      fDZResolutionFile = GetString("DZResolutionFile", "errors.root");
     87      fDZResolutionHist = GetString("DZResolutionHist", "dz");
     88      fUseDZFormula = false;
     89    }
     90  if (string (GetString("PResolutionFormula", "0.0")) != "0.0")
     91    {
     92      fPFormula->Compile(GetString("PResolutionFormula", "0.0"));
     93      fUsePFormula = true;
     94    }
     95  else
     96    {
     97      fPResolutionFile = GetString("PResolutionFile", "errors.root");
     98      fPResolutionHist = GetString("PResolutionHist", "p");
     99      fUsePFormula = false;
     100    }
     101  if (string (GetString("CtgThetaResolutionFormula", "0.0")) != "0.0")
     102    {
     103      fCtgThetaFormula->Compile(GetString("CtgThetaResolutionFormula", "0.0"));
     104      fUseCtgThetaFormula = true;
     105    }
     106  else
     107    {
     108      fCtgThetaResolutionFile = GetString("CtgThetaResolutionFile", "errors.root");
     109      fCtgThetaResolutionHist = GetString("CtgThetaResolutionHist", "ctgTheta");
     110      fUseCtgThetaFormula = false;
     111    }
     112  if (string (GetString("PhiResolutionFormula", "0.0")) != "0.0")
     113    {
     114      fPhiFormula->Compile(GetString("PhiResolutionFormula", "0.0"));
     115      fUsePhiFormula = true;
     116    }
     117  else
     118    {
     119      fPhiResolutionFile = GetString("PhiResolutionFile", "errors.root");
     120      fPhiResolutionHist = GetString("PhiResolutionHist", "phi");
     121      fUsePhiFormula = false;
     122    }
    123123
    124124  fApplyToPileUp = GetBool("ApplyToPileUp", true);
    … …  
    128128  fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles"));
    129129  fItInputArray = fInputArray->MakeIterator();
    130 
     130 
    131131  // import beamspot
    132132  try
    … …  
    137137  {
    138138    fBeamSpotInputArray = 0;
    139   }
    140 
     139  } 
     140 
    141141  // create output array
    142142
    … …  
    170170             *phiErrorHist = NULL;
    171171
    172   if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
     172  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
    173173    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    174174  else
    175175  {
    176     Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0));
     176    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At (0));
    177177    beamSpotPosition = beamSpotCandidate.Position;
    178178  }
    179179
    180   if(!fUseD0Formula)
    181   {
    182     TFile *fin = TFile::Open(fD0ResolutionFile.c_str());
    183     d0ErrorHist = (TProfile2D *)fin->Get(fD0ResolutionHist.c_str());
    184     d0ErrorHist->SetDirectory(0);
    185     fin->Close();
    186   }
    187   if(!fUseDZFormula)
    188   {
    189     TFile *fin = TFile::Open(fDZResolutionFile.c_str());
    190     dzErrorHist = (TProfile2D *)fin->Get(fDZResolutionHist.c_str());
    191     dzErrorHist->SetDirectory(0);
    192     fin->Close();
    193   }
    194   if(!fUsePFormula)
    195   {
    196     TFile *fin = TFile::Open(fPResolutionFile.c_str());
    197     pErrorHist = (TProfile2D *)fin->Get(fPResolutionHist.c_str());
    198     pErrorHist->SetDirectory(0);
    199     fin->Close();
    200   }
    201   if(!fUseCtgThetaFormula)
    202   {
    203     TFile *fin = TFile::Open(fCtgThetaResolutionFile.c_str());
    204     ctgThetaErrorHist = (TProfile2D *)fin->Get(fCtgThetaResolutionHist.c_str());
    205     ctgThetaErrorHist->SetDirectory(0);
    206     fin->Close();
    207   }
    208   if(!fUsePhiFormula)
    209   {
    210     TFile *fin = TFile::Open(fPhiResolutionFile.c_str());
    211     phiErrorHist = (TProfile2D *)fin->Get(fPhiResolutionHist.c_str());
    212     phiErrorHist->SetDirectory(0);
    213     fin->Close();
     180  if (!fUseD0Formula)
     181  {
     182     TFile *fin = TFile::Open (fD0ResolutionFile.c_str ());
     183     d0ErrorHist = (TProfile2D *) fin->Get (fD0ResolutionHist.c_str ());
     184     d0ErrorHist->SetDirectory (0);
     185     fin->Close ();
     186  }
     187  if (!fUseDZFormula)
     188  {
     189     TFile *fin = TFile::Open (fDZResolutionFile.c_str ());
     190     dzErrorHist = (TProfile2D *) fin->Get (fDZResolutionHist.c_str ());
     191     dzErrorHist->SetDirectory (0);
     192     fin->Close ();
     193  }
     194  if (!fUsePFormula)
     195  {
     196     TFile *fin = TFile::Open (fPResolutionFile.c_str ());
     197     pErrorHist = (TProfile2D *) fin->Get (fPResolutionHist.c_str ());
     198     pErrorHist->SetDirectory (0);
     199     fin->Close ();
     200  }
     201  if (!fUseCtgThetaFormula)
     202  {
     203     TFile *fin = TFile::Open (fCtgThetaResolutionFile.c_str ());
     204     ctgThetaErrorHist = (TProfile2D *) fin->Get (fCtgThetaResolutionHist.c_str ());
     205     ctgThetaErrorHist->SetDirectory (0);
     206     fin->Close ();
     207  }
     208  if (!fUsePhiFormula)
     209  {
     210     TFile *fin = TFile::Open (fPhiResolutionFile.c_str ());
     211     phiErrorHist = (TProfile2D *) fin->Get (fPhiResolutionHist.c_str ());
     212     phiErrorHist->SetDirectory (0);
     213     fin->Close ();
    214214  }
    215215
    216216  fItInputArray->Reset();
    217   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    218   {
    219 
     217  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     218  {
     219   
    220220    const TLorentzVector &momentum = candidate->Momentum;
    221221    const TLorentzVector &position = candidate->InitialPosition;
    222 
     222   
    223223    pt = momentum.Pt();
    224224    eta = momentum.Eta();
    … …  
    226226    d0 = trueD0 = candidate->D0;
    227227    dz = trueDZ = candidate->DZ;
    228 
     228     
    229229    p = trueP = candidate->P;
    230230    ctgTheta = trueCtgTheta = candidate->CtgTheta;
    231231    phi = truePhi = candidate->Phi;
    232232
    233     if(fUseD0Formula)
     233    if (fUseD0Formula)
    234234      d0Error = fD0Formula->Eval(pt, eta);
    235235    else
    236236    {
    237       Int_t xbin, ybin;
    238 
    239       xbin = pt < d0ErrorHist->GetXaxis()->GetXmax() ? d0ErrorHist->GetXaxis()->FindBin(pt) : d0ErrorHist->GetXaxis()->GetBinCenter(d0ErrorHist->GetXaxis()->GetNbins());
    240       ybin = d0ErrorHist->GetYaxis()->FindBin(TMath::Abs(eta));
    241       d0Error = d0ErrorHist->GetBinContent(xbin, ybin);
    242       if(!d0Error)
    243         d0Error = -1.0;
    244     }
    245     if(d0Error < 0.0)
     237       Int_t xbin, ybin;
     238
     239       xbin = pt < d0ErrorHist->GetXaxis ()->GetXmax () ? d0ErrorHist->GetXaxis ()->FindBin (pt)
     240            : d0ErrorHist->GetXaxis ()->GetBinCenter (d0ErrorHist->GetXaxis ()->GetNbins ());
     241       ybin = d0ErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
     242       d0Error = d0ErrorHist->GetBinContent (xbin, ybin);
     243       if (!d0Error)
     244          d0Error = -1.0;
     245    }
     246    if (d0Error < 0.0)
    246247      continue;
    247248
    248     if(fUseDZFormula)
     249    if (fUseDZFormula)
    249250      dzError = fDZFormula->Eval(pt, eta);
    250251    else
    251252    {
    252       Int_t xbin, ybin;
    253 
    254       xbin = pt < dzErrorHist->GetXaxis()->GetXmax() ? dzErrorHist->GetXaxis()->FindBin(pt) : dzErrorHist->GetXaxis()->GetBinCenter(dzErrorHist->GetXaxis()->GetNbins());
    255       ybin = dzErrorHist->GetYaxis()->FindBin(TMath::Abs(eta));
    256       dzError = dzErrorHist->GetBinContent(xbin, ybin);
    257       if(!dzError)
    258         dzError = -1.0;
    259     }
    260     if(dzError < 0.0)
     253       Int_t xbin, ybin;
     254
     255       xbin = pt < dzErrorHist->GetXaxis ()->GetXmax () ? dzErrorHist->GetXaxis ()->FindBin (pt)
     256            : dzErrorHist->GetXaxis ()->GetBinCenter (dzErrorHist->GetXaxis ()->GetNbins ());
     257       ybin = dzErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
     258       dzError = dzErrorHist->GetBinContent (xbin, ybin);
     259       if (!dzError)
     260          dzError = -1.0;
     261    }
     262    if (dzError < 0.0)
    261263      continue;
    262264
    263     if(fUsePFormula)
     265    if (fUsePFormula)
    264266      pError = fPFormula->Eval(pt, eta) * p;
    265267    else
    266268    {
    267       Int_t xbin, ybin;
    268 
    269       xbin = pt < pErrorHist->GetXaxis()->GetXmax() ? pErrorHist->GetXaxis()->FindBin(pt) : pErrorHist->GetXaxis()->GetBinCenter(pErrorHist->GetXaxis()->GetNbins());
    270       ybin = pErrorHist->GetYaxis()->FindBin(TMath::Abs(eta));
    271       pError = pErrorHist->GetBinContent(xbin, ybin) * p;
    272       if(!pError)
    273         pError = -1.0;
    274     }
    275     if(pError < 0.0)
     269       Int_t xbin, ybin;
     270
     271       xbin = pt < pErrorHist->GetXaxis ()->GetXmax () ? pErrorHist->GetXaxis ()->FindBin (pt)
     272            : pErrorHist->GetXaxis ()->GetBinCenter (pErrorHist->GetXaxis ()->GetNbins ());
     273       ybin = pErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
     274       pError = pErrorHist->GetBinContent (xbin, ybin) * p;
     275       if (!pError)
     276          pError = -1.0;
     277    }
     278    if (pError < 0.0)
    276279      continue;
    277280
    278     if(fUseCtgThetaFormula)
     281    if (fUseCtgThetaFormula)
    279282      ctgThetaError = fCtgThetaFormula->Eval(pt, eta);
    280283    else
    281284    {
    282       Int_t xbin, ybin;
    283 
    284       xbin = pt < ctgThetaErrorHist->GetXaxis()->GetXmax() ? ctgThetaErrorHist->GetXaxis()->FindBin(pt) : ctgThetaErrorHist->GetXaxis()->GetBinCenter(ctgThetaErrorHist->GetXaxis()->GetNbins());
    285       ybin = ctgThetaErrorHist->GetYaxis()->FindBin(TMath::Abs(eta));
    286       ctgThetaError = ctgThetaErrorHist->GetBinContent(xbin, ybin);
    287       if(!ctgThetaError)
    288         ctgThetaError = -1.0;
    289     }
    290     if(ctgThetaError < 0.0)
     285       Int_t xbin, ybin;
     286
     287       xbin = pt < ctgThetaErrorHist->GetXaxis ()->GetXmax () ? ctgThetaErrorHist->GetXaxis ()->FindBin (pt)
     288            : ctgThetaErrorHist->GetXaxis ()->GetBinCenter (ctgThetaErrorHist->GetXaxis ()->GetNbins ());
     289       ybin = ctgThetaErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
     290       ctgThetaError = ctgThetaErrorHist->GetBinContent (xbin, ybin);
     291       if (!ctgThetaError)
     292         ctgThetaError = -1.0;
     293    }
     294    if (ctgThetaError < 0.0)
    291295      continue;
    292296
    293     if(fUsePhiFormula)
     297    if (fUsePhiFormula)
    294298      phiError = fPhiFormula->Eval(pt, eta);
    295299    else
    296300    {
    297       Int_t xbin, ybin;
    298 
    299       xbin = pt < phiErrorHist->GetXaxis()->GetXmax() ? phiErrorHist->GetXaxis()->FindBin(pt) : phiErrorHist->GetXaxis()->GetBinCenter(phiErrorHist->GetXaxis()->GetNbins());
    300       ybin = phiErrorHist->GetYaxis()->FindBin(TMath::Abs(eta));
    301       phiError = phiErrorHist->GetBinContent(xbin, ybin);
    302       if(!phiError)
    303         phiError = -1.0;
    304     }
    305     if(phiError < 0.0)
     301       Int_t xbin, ybin;
     302
     303       xbin = pt < phiErrorHist->GetXaxis ()->GetXmax () ? phiErrorHist->GetXaxis ()->FindBin (pt)
     304            : phiErrorHist->GetXaxis ()->GetBinCenter (phiErrorHist->GetXaxis ()->GetNbins ());
     305       ybin = phiErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
     306       phiError = phiErrorHist->GetBinContent (xbin, ybin);
     307       if (!phiError)
     308          phiError = -1.0;
     309    }
     310    if (phiError < 0.0)
    306311      continue;
    307312
    308     if(fApplyToPileUp || !candidate->IsPU)
    309     {
    310       d0 = gRandom->Gaus(d0, d0Error);
    311       dz = gRandom->Gaus(dz, dzError);
    312       p = gRandom->Gaus(p, pError);
    313       ctgTheta = gRandom->Gaus(ctgTheta, ctgThetaError);
    314       phi = gRandom->Gaus(phi, phiError);
     313    if (fApplyToPileUp || !candidate->IsPU)
     314    {
     315       d0 = gRandom->Gaus(d0, d0Error);
     316       dz = gRandom->Gaus(dz, dzError);
     317       p = gRandom->Gaus(p, pError);
     318       ctgTheta = gRandom->Gaus(ctgTheta, ctgThetaError);
     319       phi = gRandom->Gaus(phi, phiError);
    315320    }
    316321
    317322    if(p < 0.0) continue;
    318     while(phi > TMath::Pi()) phi -= TMath::TwoPi();
    319     while(phi <= -TMath::Pi()) phi += TMath::TwoPi();
     323    while (phi > TMath::Pi ()) phi -= TMath::TwoPi ();
     324    while (phi <= -TMath::Pi ()) phi += TMath::TwoPi ();
    320325
    321326    mother = candidate;
    322     candidate = static_cast<Candidate *>(candidate->Clone());
     327    candidate = static_cast<Candidate*>(candidate->Clone());
    323328    candidate->D0 = d0;
    324329    candidate->DZ = dz;
    … …  
    327332    candidate->Phi = phi;
    328333
    329     theta = TMath::ACos(ctgTheta / TMath::Sqrt(1.0 + ctgTheta * ctgTheta));
    330     candidate->Momentum.SetPx(p * TMath::Cos(phi) * TMath::Sin(theta));
    331     candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta));
    332     candidate->Momentum.SetPz(p * TMath::Cos(theta));
    333     candidate->Momentum.SetE(candidate->Momentum.Pt() * TMath::CosH(eta));
    334     candidate->PT = candidate->Momentum.Pt();
    335 
    336     x = position.X();
    337     y = position.Y();
    338     z = position.Z();
    339     t = position.T();
    340     px = candidate->Momentum.Px();
    341     py = candidate->Momentum.Py();
    342     pz = candidate->Momentum.Pz();
    343     pt = candidate->Momentum.Pt();
    344 
     334    theta = TMath::ACos(ctgTheta / TMath::Sqrt (1.0 + ctgTheta * ctgTheta));
     335    candidate->Momentum.SetPx (p * TMath::Cos (phi) * TMath::Sin (theta));
     336    candidate->Momentum.SetPy (p * TMath::Sin (phi) * TMath::Sin (theta));
     337    candidate->Momentum.SetPz (p * TMath::Cos (theta));
     338    candidate->Momentum.SetE (candidate->Momentum.Pt () * TMath::CosH (eta));
     339    candidate->PT = candidate->Momentum.Pt ();
     340
     341    x = position.X ();
     342    y = position.Y ();
     343    z = position.Z ();
     344    t = position.T ();
     345    px = candidate->Momentum.Px ();
     346    py = candidate->Momentum.Py ();
     347    pz = candidate->Momentum.Pz ();
     348    pt = candidate->Momentum.Pt ();
     349   
    345350    // -- solve for delta: d0' = ( (x+delta)*py' - (y+delta)*px' )/pt'
    346 
    347     candidate->InitialPosition.SetX(x + ((px * y - py * x + d0 * pt) / (py - px)));
    348     candidate->InitialPosition.SetY(y + ((px * y - py * x + d0 * pt) / (py - px)));
    349     x = candidate->InitialPosition.X();
    350     y = candidate->InitialPosition.Y();
    351     candidate->InitialPosition.SetZ(z + ((pz * (px * (x - beamSpotPosition.X()) + py * (y - beamSpotPosition.Y())) + pt * pt * (dz - z)) / (pt * pt)));
    352     z = candidate->InitialPosition.Z();
     351   
     352    candidate->InitialPosition.SetX (x + ((px * y - py * x + d0 * pt) / (py - px)));
     353    candidate->InitialPosition.SetY (y + ((px * y - py * x + d0 * pt) / (py - px)));
     354    x = candidate->InitialPosition.X ();
     355    y = candidate->InitialPosition.Y ();
     356    candidate->InitialPosition.SetZ (z + ((pz * (px * (x - beamSpotPosition.X ()) + py * (y - beamSpotPosition.Y ())) + pt * pt * (dz - z)) / (pt * pt)));
     357    z = candidate->InitialPosition.Z ();
    353358
    354359    candidate->InitialPosition.SetT(t);
    … …  
    361366    q = candidate->Charge;
    362367
    363     r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
     368    r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
    364369    phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    365370
    366371    // 2. helix axis coordinates
    367     x_c = x + r * TMath::Sin(phi_0);
    368     y_c = y - r * TMath::Cos(phi_0);
     372    x_c = x + r*TMath::Sin(phi_0);
     373    y_c = y - r*TMath::Cos(phi_0);
    369374    r_c = TMath::Hypot(x_c, y_c);
    370375
    371376    rcu = TMath::Abs(r);
    372     rc2 = r_c * r_c;
     377    rc2 = r_c*r_c;
    373378
    374379    // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
    375     xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
     380    xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
    376381    xd = (rc2 > 0.0) ? xd / rc2 : -999;
    377     yd = y_c * (-rcu * r_c + rc2);
     382    yd = y_c*(-rcu*r_c + rc2);
    378383    yd = (rc2 > 0.0) ? yd / rc2 : -999;
    379     zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
    380 
    381     candidate->Xd = xd * 1.0E3;
    382     candidate->Yd = yd * 1.0E3;
    383     candidate->Zd = zd * 1.0E3;
    384 
    385     if(fApplyToPileUp || !candidate->IsPU)
    386     {
    387       candidate->ErrorD0 = d0Error;
    388       candidate->ErrorDZ = dzError;
    389       candidate->ErrorP = pError;
    390       candidate->ErrorCtgTheta = ctgThetaError;
    391       candidate->ErrorPhi = phiError;
    392       candidate->ErrorPT = ptError(p, ctgTheta, pError, ctgThetaError);
    393       candidate->TrackResolution = pError / p;
    394     }
    395 
     384    zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
     385
     386    candidate->Xd = xd*1.0E3;
     387    candidate->Yd = yd*1.0E3;
     388    candidate->Zd = zd*1.0E3;
     389
     390    if (fApplyToPileUp || !candidate->IsPU)
     391    {
     392       candidate->ErrorD0 = d0Error;
     393       candidate->ErrorDZ = dzError;
     394       candidate->ErrorP = pError;
     395       candidate->ErrorCtgTheta = ctgThetaError;
     396       candidate->ErrorPhi = phiError;
     397       candidate->ErrorPT = ptError (p, ctgTheta, pError, ctgThetaError);
     398       candidate->TrackResolution = pError/p;
     399    }
     400   
    396401    candidate->AddCandidate(mother);
    397402    fOutputArray->Add(candidate);
    … …  
    401406}
    402407
    403 Double_t TrackSmearing::ptError(const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta)
     408Double_t TrackSmearing::ptError (const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta)
    404409{
    405410  Double_t a, b;
    406411  a = (p * p * ctgTheta * ctgTheta * dCtgTheta * dCtgTheta) / ((ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1));
    407412  b = (dP * dP) / (ctgTheta * ctgTheta + 1);
    408   return sqrt(a + b);
    409 }
    410 
    411 //------------------------------------------------------------------------------
     413  return sqrt (a + b);
     414}
     415
     416//------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.