Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/TrackSmearing.cc

    rcde9f31 r341014c  
    1313#include "classes/DelphesFormula.h"
    1414
     15#include "ExRootAnalysis/ExRootClassifier.h"
     16#include "ExRootAnalysis/ExRootFilter.h"
    1517#include "ExRootAnalysis/ExRootResult.h"
    16 #include "ExRootAnalysis/ExRootFilter.h"
    17 #include "ExRootAnalysis/ExRootClassifier.h"
    18 
     18
     19#include "TDatabasePDG.h"
     20#include "TFile.h"
     21#include "TFormula.h"
     22#include "TLorentzVector.h"
    1923#include "TMath.h"
     24#include "TObjArray.h"
     25#include "TProfile2D.h"
     26#include "TRandom3.h"
    2027#include "TString.h"
    21 #include "TFormula.h"
    22 #include "TRandom3.h"
    23 #include "TObjArray.h"
    24 #include "TDatabasePDG.h"
    25 #include "TLorentzVector.h"
    26 #include "TFile.h"
    27 #include "TProfile2D.h"
    28 
    29 #include <algorithm>
    30 #include <stdexcept>
     28
     29#include <algorithm>
    3130#include <iostream>
    3231#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)
    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)
     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)
    247246      continue;
    248247
    249     if (fUseDZFormula)
     248    if(fUseDZFormula)
    250249      dzError = fDZFormula->Eval(pt, eta);
    251250    else
    252251    {
    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)
     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)
    263261      continue;
    264262
    265     if (fUsePFormula)
     263    if(fUsePFormula)
    266264      pError = fPFormula->Eval(pt, eta) * p;
    267265    else
    268266    {
    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)
     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)
    279276      continue;
    280277
    281     if (fUseCtgThetaFormula)
     278    if(fUseCtgThetaFormula)
    282279      ctgThetaError = fCtgThetaFormula->Eval(pt, eta);
    283280    else
    284281    {
    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)
     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)
    295291      continue;
    296292
    297     if (fUsePhiFormula)
     293    if(fUsePhiFormula)
    298294      phiError = fPhiFormula->Eval(pt, eta);
    299295    else
    300296    {
    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)
     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)
    311306      continue;
    312307
    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);
     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);
    320315    }
    321316
    322317    if(p < 0.0) continue;
    323     while (phi > TMath::Pi ()) phi -= TMath::TwoPi ();
    324     while (phi <= -TMath::Pi ()) phi += TMath::TwoPi ();
     318    while(phi > TMath::Pi()) phi -= TMath::TwoPi();
     319    while(phi <= -TMath::Pi()) phi += TMath::TwoPi();
    325320
    326321    mother = candidate;
    327     candidate = static_cast<Candidate*>(candidate->Clone());
     322    candidate = static_cast<Candidate *>(candidate->Clone());
    328323    candidate->D0 = d0;
    329324    candidate->DZ = dz;
    … …  
    332327    candidate->Phi = phi;
    333328
    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    
     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
    350345    // -- solve for delta: d0' = ( (x+delta)*py' - (y+delta)*px' )/pt'
    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 ();
     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();
    358353
    359354    candidate->InitialPosition.SetT(t);
    … …  
    366361    q = candidate->Charge;
    367362
    368     r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
     363    r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
    369364    phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    370365
    371366    // 2. helix axis coordinates
    372     x_c = x + r*TMath::Sin(phi_0);
    373     y_c = y - r*TMath::Cos(phi_0);
     367    x_c = x + r * TMath::Sin(phi_0);
     368    y_c = y - r * TMath::Cos(phi_0);
    374369    r_c = TMath::Hypot(x_c, y_c);
    375370
    376371    rcu = TMath::Abs(r);
    377     rc2 = r_c*r_c;
     372    rc2 = r_c * r_c;
    378373
    379374    // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
    380     xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
     375    xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
    381376    xd = (rc2 > 0.0) ? xd / rc2 : -999;
    382     yd = y_c*(-rcu*r_c + rc2);
     377    yd = y_c * (-rcu * r_c + rc2);
    383378    yd = (rc2 > 0.0) ? yd / rc2 : -999;
    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    
     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
    401396    candidate->AddCandidate(mother);
    402397    fOutputArray->Add(candidate);
    … …  
    406401}
    407402
    408 Double_t TrackSmearing::ptError (const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta)
     403Double_t TrackSmearing::ptError(const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta)
    409404{
    410405  Double_t a, b;
    411406  a = (p * p * ctgTheta * ctgTheta * dCtgTheta * dCtgTheta) / ((ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1));
    412407  b = (dP * dP) / (ctgTheta * ctgTheta + 1);
    413   return sqrt (a + b);
    414 }
    415 
    416 //------------------------------------------------------------------------------
     408  return sqrt(a + b);
     409}
     410
     411//------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.