PWG3/hfe/AliHFEpriVtx.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 /* $Id$ */
  17
  18 //
  19 // QA class of primary vertex study for Heavy Flavor electrons
  20 // this has functionality to reject electrons from primary vertex
  21 // and check primary vertex characteristics
  22 //
  23 // Authors:
  24 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
  25 //
  26
  27
  28 #include <TH2F.h>
  29 #include <TParticle.h>
  30
  31 #include "AliESDEvent.h"
  32 #include <AliESDtrack.h>
  33 #include <AliMCEvent.h>
  34
  35 #include <AliLog.h>
  36 #include <AliKFParticle.h>
  37 #include <AliKFVertex.h>
  38
  39 #include "AliHFEpriVtx.h"
  40
  41
  42 ClassImp(AliHFEpriVtx)
  43
  44 //_______________________________________________________________________________________________
  45 AliHFEpriVtx::AliHFEpriVtx():
  46         fESD1(0x0)
  47         ,fMCEvent(0x0)
  48         ,fNtrackswoPid(0)
  49         ,fHNtrackswoPid(0x0)
  50         ,fNESDprimVtxContributor(0x0)
  51         ,fNESDprimVtxIndices(0x0)
  52         ,fDiffDCAvsPt(0x0)
  53         ,fDiffDCAvsNt(0x0)
  54         ,fNsectrk2prim(0)
  55         ,fPVxRe(-999.)
  56         ,fPVyRe(-999.)
  57         ,fPVzRe(-999.)
  58 {
  59         //
  60         // Default constructor
  61         //
  62
  63         Init();
  64
  65 }
  66
  67 //_______________________________________________________________________________________________
  68 AliHFEpriVtx::AliHFEpriVtx(const AliHFEpriVtx &p):
  69          TObject(p)
  70         ,fESD1(0x0)
  71         ,fMCEvent(0x0)
  72         ,fNtrackswoPid(p.fNtrackswoPid)
  73         ,fHNtrackswoPid(0x0)
  74         ,fNESDprimVtxContributor(0x0)
  75         ,fNESDprimVtxIndices(0x0)
  76         ,fDiffDCAvsPt(0x0)
  77         ,fDiffDCAvsNt(0x0)
  78         ,fNsectrk2prim(p.fNsectrk2prim)
  79         ,fPVxRe(p.fPVxRe)
  80         ,fPVyRe(p.fPVyRe)
  81         ,fPVzRe(p.fPVzRe)
  82 {
  83         //
  84         // Copy constructor
  85         //
  86 }
  87
  88 //_______________________________________________________________________________________________
  89 AliHFEpriVtx&
  90 AliHFEpriVtx::operator=(const AliHFEpriVtx &)
  91 {
  92         //
  93         // Assignment operator
  94         //
  95
  96         AliInfo("Not yet implemented.");
  97         return *this;
  98 }
  99
 100 //_______________________________________________________________________________________________
 101 AliHFEpriVtx::~AliHFEpriVtx()
 102 {
 103         //
 104         // Destructor
 105         //
 106
 107         AliInfo("Analysis Done.");
 108 }
 109
 110 //__________________________________________
 111 void AliHFEpriVtx::Init()
 112 {
 113         //
 114         // initialize counters
 115         //
 116
 117         fNtrackswoPid = 0;
 118         for (int i=0; i<10; i++){
 119                 fPrimVtx[i].fNtrackCount = 0 ;
 120                 fPrimVtx[i].fNprimVtxContributorCount = 0 ;
 121         }
 122 }
 123
 124 //_______________________________________________________________________________________________
 125 void AliHFEpriVtx::CreateHistograms(TString hnopt)
 126 {
 127         //
 128         // create histograms
 129         //
 130
 131         fkSourceLabel[kAll]="all";
 132         fkSourceLabel[kDirectCharm]="directCharm";
 133         fkSourceLabel[kDirectBeauty]="directBeauty";
 134         fkSourceLabel[kBeautyCharm]="beauty2charm";
 135         fkSourceLabel[kGamma]="gamma";
 136         fkSourceLabel[kPi0]="pi0";
 137         fkSourceLabel[kElse]="others";
 138         fkSourceLabel[kBeautyGamma]="beauty22gamma";
 139         fkSourceLabel[kBeautyPi0]="beauty22pi0";
 140         fkSourceLabel[kBeautyElse]="beauty22others";
 141
 142
 143         TString hname;
 144         for (Int_t isource = 0; isource < 10; isource++ ){
 145
 146            hname=hnopt+"ntracks_"+fkSourceLabel[isource];
 147            fPrimVtx[isource].fNtracks = new TH1F(hname,hname,50,0,50);
 148            hname=hnopt+"nPrimVtxContributor_"+fkSourceLabel[isource];
 149            fPrimVtx[isource].fNprimVtxContributor = new TH1F(hname,hname,100,0,100);
 150            hname=hnopt+"PtElec_"+fkSourceLabel[isource];
 151            fPrimVtx[isource].fPtElec = new TH1F(hname,hname,250,0,50);
 152            hname=hnopt+"PtElecContributor_"+fkSourceLabel[isource];
 153            fPrimVtx[isource].fPtElecContributor = new TH1F(hname,hname,250,0,50);
 154
 155         }
 156
 157         hname=hnopt+"ntrackswopid";
 158         fHNtrackswoPid = new TH1F(hname,hname,50,0,50);
 159         hname=hnopt+"nESDprimVtxContributor";
 160         fNESDprimVtxContributor = new TH1I(hname,hname,100,0,100);
 161         hname=hnopt+"nESDprimVtxIndices";
 162         fNESDprimVtxIndices= new TH1I(hname,hname,100,0,100);
 163         hname=hnopt+"diffDCAvsPt";
 164         fDiffDCAvsPt = new TH2F(hname,hname,250,0,50,500,0,1);
 165         hname=hnopt+"diffDCAvsNt";
 166         fDiffDCAvsNt = new TH2F(hname,hname,100,0,100,500,0,1);
 167
 168 }
 169
 170 //_______________________________________________________________________________________________
 171 void AliHFEpriVtx::CountNtracks(Int_t sourcePart, Int_t recpid, Double_t recprob)
 172 {
 173         //
 174         // count number of tracks passed certain cuts
 175         //
 176
 177         fNtrackswoPid++;
 178
 179         if (!recpid && recprob>0.5)
 180         fPrimVtx[kAll].fNtrackCount++;
 181         if(sourcePart<0) return;
 182         fPrimVtx[sourcePart].fNtrackCount++;
 183
 184 }
 185
 186 //_______________________________________________________________________________________________
 187 void AliHFEpriVtx::FillNtracks()
 188 {
 189         //
 190         // count number of tracks passed certain cuts
 191         //
 192
 193         fHNtrackswoPid->Fill(fNtrackswoPid);
 194         for (int i=0; i<10; i++){
 195           fPrimVtx[i].fNtracks->Fill(fPrimVtx[i].fNtrackCount);
 196         }
 197
 198 }
 199
 200 //_______________________________________________________________________________________________
 201 Int_t AliHFEpriVtx::GetMCPID(AliESDtrack const *track)
 202 {
 203         //
 204         // get MC pid
 205         //
 206
 207         AliMCParticle *mctrack = NULL;
 208         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return 0;
 209         TParticle *mcpart = mctrack->Particle();
 210
 211         if ( !mcpart ) return 0;
 212         Int_t pdgCode = mcpart->GetPdgCode();
 213
 214         return pdgCode;
 215 }
 216
 217 //_______________________________________________________________________________________________
 218 void AliHFEpriVtx::GetNPriVxtContributor()
 219 {
 220         //
 221         // count number of primary vertex contributor
 222         //
 223
 224         const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 225         fNESDprimVtxContributor->Fill(primvtx->GetNContributors());
 226         fNESDprimVtxIndices->Fill(primvtx->GetNIndices());
 227 }
 228
 229 //_______________________________________________________________________________________________
 230 void AliHFEpriVtx::CountPriVxtElecContributor(AliESDtrack *ESDelectron, Int_t sourcePart, Int_t recpid, Double_t recprob)
 231 {
 232         //
 233         // count number of electrons contributing to the primary vertex
 234         //
 235
 236
 237         if (recpid || recprob<0.5) return;
 238
 239         // get track id of our selected electron
 240         Int_t elecTrkID = ESDelectron->GetID();
 241
 242         AliMCParticle *mctrack = NULL;
 243         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(ESDelectron->GetLabel()))))) return;
 244         TParticle *mcpart = mctrack->Particle();
 245
 246
 247         if(!mcpart){
 248           AliDebug(1, "no mc particle, return\n");
 249           return;
 250         }
 251
 252         AliKFParticle::SetField(fESD1->GetMagneticField());
 253         AliKFParticle kfElectron(*ESDelectron,11);
 254
 255         // prepare kfprimary vertex
 256         AliKFVertex kfESDprimary;
 257
 258         // Reconstructed Primary Vertex (with ESD tracks)
 259         const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 260 //         Int_t nt = primvtx->GetNContributors();
 261         Int_t n=primvtx->GetNIndices();
 262
 263         if (n>0 && primvtx->GetStatus()){
 264
 265                 kfESDprimary = AliKFVertex(*primvtx);
 266                 Double_t dcaBFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
 267
 268                 UShort_t *priIndex = primvtx->GetIndices();
 269
 270                 fPrimVtx[kAll].fPtElec->Fill(mcpart->Pt());
 271                 if(sourcePart>=0) fPrimVtx[sourcePart].fPtElec->Fill(mcpart->Pt());
 272
 273                 for (Int_t i=0;i<n;i++){
 274
 275                         Int_t idx = Int_t(priIndex[i]);
 276                         if (idx == elecTrkID){
 277                                 fPrimVtx[kAll].fNprimVtxContributorCount++;
 278                                 fPrimVtx[kAll].fPtElecContributor->Fill(mcpart->Pt());
 279                                 if(sourcePart<0) continue;
 280                                 fPrimVtx[sourcePart].fNprimVtxContributorCount++;
 281                                 fPrimVtx[sourcePart].fPtElecContributor->Fill(mcpart->Pt());
 282
 283                                 kfESDprimary -= kfElectron;
 284                                 Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
 285                                 fDiffDCAvsPt->Fill(mcpart->Pt(),dcaBFrmElec-dcaAFrmElec);
 286                                 fDiffDCAvsNt->Fill(n,dcaBFrmElec-dcaAFrmElec);
 287                         }
 288                 }
 289         }
 290
 291 }
 292
 293 //_______________________________________________________________________________________________
 294 void AliHFEpriVtx::FillNprimVtxContributor() const
 295 {
 296         //
 297         // Fill histogram with number of electrons contributing to the primary vertex
 298         //
 299
 300         for (int i=0; i<10; i++){
 301           fPrimVtx[i].fNprimVtxContributor->Fill(fPrimVtx[i].fNprimVtxContributorCount);
 302         }
 303
 304 }
 305
 306 //_______________________________________________________________________________________________
 307 Double_t AliHFEpriVtx::GetDistanceFromRecalVertexXY(AliESDtrack * const ESDelectron)
 308 {
 309         //
 310         // return recalculated DCA after removing input track from the primary vertex
 311         //
 312
 313         // get track id of our selected electron
 314         Int_t elecTrkID = ESDelectron->GetID();
 315
 316         AliKFParticle::SetField(fESD1->GetMagneticField());
 317         AliKFParticle kfElectron(*ESDelectron,11);
 318
 319         // prepare kfprimary vertex
 320         AliKFVertex kfESDprimary;
 321
 322         // Reconstructed Primary Vertex (with ESD tracks)
 323         const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 324         Int_t n=primvtx->GetNIndices();
 325
 326         if (n>0 && primvtx->GetStatus()){
 327
 328                 kfESDprimary = AliKFVertex(*primvtx);
 329                 UShort_t *priIndex = primvtx->GetIndices();
 330
 331                 for (Int_t i=0;i<n;i++){
 332
 333                         Int_t idx = Int_t(priIndex[i]);
 334                         if (idx == elecTrkID){
 335                                 kfESDprimary -= kfElectron;
 336                                 Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
 337
 338                                 return dcaAFrmElec;
 339                         }
 340                 }
 341         }
 342               return -1;
 343
 344 }
 345
 346 void AliHFEpriVtx::RecalcPrimvtx(Int_t nkftrk, const Int_t * const trkid, const AliKFParticle * const kftrk)
 347 {
 348         //
 349         // recalculate primary vertex after removing the input track
 350         //
 351
 352         const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 353
 354         AliKFVertex kfESDprimary;
 355         Int_t n = primvtx->GetNIndices();
 356         fNsectrk2prim = 0;
 357         fPVxRe = -999.;
 358         fPVyRe = -999.;
 359         fPVyRe = -999.;
 360
 361         if (n>0 && primvtx->GetStatus()){
 362           kfESDprimary = AliKFVertex(*primvtx);
 363           UShort_t *priIndex = primvtx->GetIndices();
 364           for (Int_t j=0; j<nkftrk; j++){
 365             for (Int_t i=0;i<n;i++){
 366               Int_t idx = Int_t(priIndex[i]);
 367               if (idx == trkid[j]){
 368                 kfESDprimary -= kftrk[j];
 369                 fNsectrk2prim++;
 370               }
 371             }
 372           }
 373         }
 374
 375         fPVxRe = kfESDprimary.GetX();
 376         fPVyRe = kfESDprimary.GetY();
 377         fPVzRe = kfESDprimary.GetZ();
 378
 379 }
 380
 381
 382 //_______________________________________________________________________________________________
 383 void AliHFEpriVtx::RecalcPrimvtx(AliESDtrack * const ESDelectron)
 384 {
 385         //
 386         // recalculate primary vertex after removing the input track
 387         //
 388
 389         // get track id of our selected electron
 390         Int_t elecTrkID = ESDelectron->GetID();
 391
 392         AliKFParticle::SetField(fESD1->GetMagneticField());
 393         AliKFParticle kfElectron(*ESDelectron,11);
 394
 395         const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 396
 397         AliKFVertex kfESDprimary;
 398         Int_t n = primvtx->GetNIndices();
 399         fPVxRe = -999.;
 400         fPVyRe = -999.;
 401         fPVyRe = -999.;
 402
 403         if (n>0 && primvtx->GetStatus()){
 404           kfESDprimary = AliKFVertex(*primvtx);
 405           UShort_t *priIndex = primvtx->GetIndices();
 406           for (Int_t i=0;i<n;i++){
 407             Int_t idx = Int_t(priIndex[i]);
 408             if (idx == elecTrkID){
 409               kfESDprimary -= kfElectron;
 410             }
 411           }
 412         }
 413
 414         fPVxRe = kfESDprimary.GetX();
 415         fPVyRe = kfESDprimary.GetY();
 416         fPVzRe = kfESDprimary.GetZ();
 417
 418 }