PWG3/hfe/AliHFEdca.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Class for impact parameter (DCA) of charged particles
  17 // Study resolution and pull: prepare for beauty study
  18 //
  19 // Authors:
  20 //   Hongyan Yang <hongyan@physi.uni-heidelberg.de>
  21 //   Carlo Bombonati <carlo.bombonati@cern.ch>
  22 //
  23
  24 #include "TMath.h"
  25 #include "TH1F.h"
  26 #include "TList.h"
  27 #include <TParticle.h>
  28 #include "AliMCParticle.h"
  29 #include "AliESDtrack.h"
  30 #include "AliESDEvent.h"
  31 #include "AliMCEvent.h"
  32
  33 #include "AliHFEdca.h"
  34
  35 ClassImp(AliHFEdca)
  36
  37 //________________________________________________________________________
  38 const Char_t* AliHFEdca::fgkParticles[12] = {
  39   // particles name
  40   "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
  41   "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
  42   "allNegative", "allPositive"
  43 };
  44 //________________________________________________________________________
  45 const Int_t AliHFEdca::fgkColorPart[12] = {
  46   // colors assigned to particles
  47   kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
  48   kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
  49   kBlack, kGray+1
  50 };
  51
  52
  53
  54 //________________________________________________________________________
  55 const Float_t AliHFEdca::fgkPtIntv[44] = {
  56   // define pT bins
  57   0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
  58   1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8,
  59   3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 5.0, 5.5, 6.0, 6.5,
  60   7.0, 8.0, 9.0, 10., 11., 12., 13., 14., 15., 16.,
  61   17., 18., 19., 20.};
  62
  63 //________________________________________________________________________
  64 const Char_t *AliHFEdca::fgkDcaVar[2] = {
  65   "deltaDcaXY",  "deltaDcaZ"};
  66
  67 //________________________________________________________________________
  68 const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
  69   ";residual #Delta(d_{xy}) [#mum];couts", ";residual #Delta(d_{z}) [#mum];counts"};
  70
  71
  72 //________________________________________________________________________
  73 const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
  74   "pullDcaXY", "pullDcaZ"
  75 };
  76
  77 //________________________________________________________________________
  78 const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
  79   ";(dca_{xy}^{ESD}-dca_{xy}^{MC})/(#sigma_{track}#oplus#sigma_{Vxy});counts",
  80   ";(dca_{z}^{ESD}-dca_{z}^{MC})/(#sigma_{track}#oplus#sigma_{Vz}); counts"
  81 };
  82
  83 //________________________________________________________________________
  84 AliHFEdca::AliHFEdca():
  85    fResidualList(0x0)
  86   , fPullList(0x0)
  87
  88 {
  89   // default constructor
  90
  91 }
  92
  93 //________________________________________________________________________
  94 AliHFEdca::AliHFEdca(const AliHFEdca &dca):
  95   TObject(dca)
  96   , fResidualList(0x0)
  97   , fPullList(0x0)
  98
  99 {
 100   // copy constructor
 101
 102 }
 103 //_______________________________________________________________________________________________
 104 AliHFEdca&
 105 AliHFEdca::operator=(const AliHFEdca &)
 106 {
 107   //
 108   // Assignment operator
 109   //
 110
 111   Printf("Not yet implemented.");
 112   return *this;
 113 }
 114
 115 //________________________________________________________________________
 116 AliHFEdca::~AliHFEdca()
 117 {
 118   // default destructor
 119
 120   for(Int_t j=0; j<kNParticles; j++){
 121     for(Int_t i=0; i<kNPtBins; i++){
 122       if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
 123       if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
 124       if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
 125       if(fHistDcaXYPull[j][i]) delete fHistDcaZPull[j][i];
 126     }
 127   }
 128
 129   if(fResidualList) delete fResidualList;
 130   if(fPullList) delete fPullList;
 131
 132   Printf("analysis done\n");
 133 }
 134
 135 //________________________________________________________________________
 136 void AliHFEdca::InitAnalysis(){
 137
 138   Printf("initialize analysis\n");
 139
 140 }
 141
 142
 143 //________________________________________________________________________
 144 void AliHFEdca::PostAnalysis() const
 145 {
 146   // do fit
 147   // moved to dcaPostAnalysis.C
 148
 149 }
 150 //________________________________________________________________________
 151 void AliHFEdca::CreateHistogramsResidual(TList *residualList){
 152   // define histogram
 153   // 1. residual
 154
 155   // for residuals
 156   fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
 157   fHistDcaZRes[kNParticles][kNPtBins]=0x0;
 158
 159   const Int_t nBins = 1000;
 160   const Float_t maxXYBin = 1000.;
 161   const Float_t maxZBin = 1000.;
 162
 163
 164   for(Int_t k=0; k<kNDcaVar; k++){
 165     TString histTitle((const char*)fgkDcaVarTitle[k]);
 166
 167     for(Int_t j=0; j<kNParticles; j++){
 168       for(Int_t i=0; i<kNPtBins; i++){
 169
 170         TString histName((const char*)fgkParticles[j]);
 171
 172         histName += Form("_%s_pT-%.1f-%.1f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 173
 174         if(k==0){
 175           fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 176           fHistDcaXYRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
 177         }
 178         if(k==1){
 179           fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 180           fHistDcaZRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
 181         }
 182       } // 43 pt bins
 183     } //12 nparticles
 184   } // 2 dca var
 185
 186   //  TList *fResidualList = 0;
 187   residualList->SetOwner();
 188   residualList->SetName("residual");
 189   for(Int_t iPart=0; iPart<kNParticles; iPart++){
 190     for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 191       residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
 192       residualList->Add(fHistDcaZRes[iPart][iPtBin]);
 193     } // loop over pt bins
 194   }  // loop over particles (pos, neg)
 195
 196
 197
 198
 199 }
 200
 201
 202 //________________________________________________________________________
 203 void AliHFEdca::CreateHistogramsPull(TList *pullList){
 204   // define histogram
 205   // 2. pull
 206
 207   const Int_t nBins = 1000;
 208   const Float_t maxXYBin = 20.;
 209   const Float_t maxZBin = 20.;
 210
 211
 212   // for pull -----------------------------------------------------------------------
 213   fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
 214   fHistDcaZPull[kNParticles][kNPtBins]=0x0;
 215
 216
 217   for(Int_t k=0; k<kNDcaVar; k++){
 218     TString histTitle((const char*)fgkPullDcaVarTitle[k]);
 219
 220     for(Int_t j=0; j<kNParticles; j++){
 221       for(Int_t i=0; i<kNPtBins; i++){
 222
 223         TString histName((const char*)fgkParticles[j]);
 224
 225         histName += Form("_%s_pT-%.1f-%.1f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 226
 227         if(k==0){
 228           fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
 229           fHistDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 230         }
 231         if(k==1){
 232           fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
 233           fHistDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 234         }
 235       } // 43 pt bins
 236     } //6 nparticles
 237   } // 2 dca var
 238
 239   //  TList *fPullList = 0;
 240   pullList->SetOwner();
 241   pullList->SetName("pull");
 242   for(Int_t iPart=0; iPart<kNParticles; iPart++){
 243     for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 244       pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
 245       pullList->Add(fHistDcaZPull[iPart][iPtBin]);
 246     } // loop over pt bins
 247   }  // loop over particles (pos, neg)
 248
 249
 250
 251 }
 252
 253
 254 //_______________________________________________________________________________________________
 255 void AliHFEdca::FillHistograms(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
 256 {
 257 // filling historgams track by track
 258
 259 // obtaining reconstructed dca ------------------------------------------------------------------
 260   Float_t esdpx = track->Px();
 261   Float_t esdpy = track->Py();
 262   Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
 263   Float_t b[2];  // dca in cm
 264   Float_t bCov[3];  // covariance matrix
 265   track->GetImpactParameters(b,bCov);
 266
 267 // obtaining errors of dca ------------------------------------------------------------------
 268   const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
 269   Float_t magneticField = 5;  // initialized as 5kG
 270   magneticField = esdEvent->GetMagneticField();  // in kG
 271
 272   Double_t dz[2];   // error of dca in cm
 273   Double_t covardz[3];
 274   track->PropagateToDCA(primVtx,magneticField, 1000., dz, covardz);
 275
 276   // calculate mcDca ------------------------------------------------------------------
 277
 278   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
 279   TParticle *part = mctrack->Particle();
 280   Int_t pdg = part->GetPdgCode();
 281   Int_t charge = 1;
 282   if(pdg==kPDGelectron || pdg==kPDGmuon
 283      || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
 284
 285   Float_t vx = part->Vx();  // in cm
 286   Float_t vy = part->Vy();  // in cm
 287   Float_t vz = part->Vz();   // in cm
 288
 289   Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
 290
 291   Float_t mcpx = part->Px();
 292   Float_t mcpy = part->Py();
 293   Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
 294
 295   const Float_t conv[2] = {1.783/1.6, 2.99792458};
 296   Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla
 297
 298   Float_t nx = esdpx/mcpt;
 299   Float_t ny = esdpy/mcpt;
 300
 301   Float_t radius;
 302   radius = TMath::Abs(radiusMc);
 303   Float_t mcDcaXY = (radius - TMath::Sqrt(vxy*vxy/100./100. + radius*radius + 2*radius*charge*(vx*ny-vy*nx)/100.)) ;  // in meters
 304
 305 //   printf("magnetic Field = %.3f \t", magneticField);
 306 //   printf("pt=esd  %.3f/mc  %.3f GeV/c, radius=esd  %.3f/mc  %.3f meter \n", esdpt, mcpt, radiusEsd, radiusMc);
 307 //   printf("mcDcaXY=%.5f micron, esdDcaXY=%.5f micron \t ",  mcDcaXY*1.e6, b[0]*1e4);
 308 //   printf("mcDcaZ=%.5f micron, esdDcaZ=%.5f micron\n\n", vz*1e6, b[1]*1e4);
 309
 310   Double_t mcDca[2] = {mcDcaXY*100, vz};  // in cm
 311
 312   Double_t residual[2] = {0, 0};
 313   Double_t pull[2] = {0, 0};
 314   Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
 315   for(Int_t i=0; i<2; i++){
 316     residual[i] = dz[i] - mcDca[i]; // in centimeters
 317     if(error[i]!=0)pull[i] = residual[i]/error[i];   // unitless
 318     //    printf("error[%d]=%.6f  residual[%d]=%.6f   pull[%d]=%.6f\n", i, error[i], i, residual[i], i, pull[i]);
 319   }
 320
 321   Int_t fPdgParticle[10] = {
 322     kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
 323     -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
 324
 325   for(Int_t iPart=0; iPart<kNParticles-2; iPart++){
 326
 327     // identified ones
 328     for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 329       if(pdg==fPdgParticle[iPart] && (esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1])) {
 330         fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);  // in microns
 331         fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);   // in microns
 332         fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
 333         fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
 334       }
 335       else
 336         continue;
 337     }
 338   }
 339
 340   // for charged particles
 341   for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 342     if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
 343       Int_t iPart = 10;
 344       if(charge>0) iPart = 11;
 345       fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
 346       fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
 347       fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
 348       fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
 349     }
 350     else
 351       continue;
 352   }
 353
 354 }
 355