PWG3/vertexingHF/AliCFVertexingHF2Prong.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 //-----------------------------------------------------------------------
  17 // Class for HF corrections as a function of many variables and step
  18 // Author : C. Zampolli, CERN
  19 // D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
  20 // Base class for HF Unfolding - agrelli@uu.nl
  21 //-----------------------------------------------------------------------
  22
  23 #include "AliAODRecoDecayHF2Prong.h"
  24 #include "AliAODMCParticle.h"
  25 #include "AliAODEvent.h"
  26 #include "TClonesArray.h"
  27 #include "AliCFVertexingHF.h"
  28 #include "AliESDtrack.h"
  29 #include "TDatabasePDG.h"
  30
  31 #include "AliCFVertexingHF2Prong.h"
  32 #include "AliCFContainer.h"
  33
  34 ClassImp(AliCFVertexingHF2Prong)
  35
  36
  37 //_________________________________________
  38   AliCFVertexingHF2Prong::AliCFVertexingHF2Prong(TClonesArray *mcArray, UShort_t originDselection):
  39   AliCFVertexingHF(mcArray, originDselection)
  40
  41 {
  42
  43   SetNProngs(2);
  44 }
  45
  46
  47 //_____________________________________
  48 AliCFVertexingHF2Prong& AliCFVertexingHF2Prong::operator=(const AliCFVertexingHF2Prong& c){
  49
  50   if  (this != &c) {
  51
  52           AliCFVertexingHF::operator=(c);
  53
  54   }
  55     return *this;
  56
  57
  58 }
  59
  60 //__________________________________________
  61 Bool_t AliCFVertexingHF2Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
  62
  63   Bool_t bSignAssoc = kFALSE;
  64
  65   fRecoCandidate = recoCand;
  66
  67   if (!fRecoCandidate) {
  68     printf("fRecoCandidate not found, problem in assignement\n");
  69     return bSignAssoc;
  70   }
  71
  72   if (fRecoCandidate->GetPrimaryVtx()) printf("fReco Candidate has a pointer to PrimVtx\n");
  73   if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
  74
  75   Int_t pdgCand = 421;
  76   Int_t pdgDgD0toKpi[2]={321,211};
  77   Int_t nentries = fmcArray->GetEntriesFast();
  78
  79   printf ("nentries = %d\n", nentries);
  80
  81   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,2,pdgDgD0toKpi);
  82
  83   if (mcLabel == -1) return bSignAssoc;
  84   SetMCLabel(mcLabel);
  85   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
  86
  87   if (!fmcPartCandidate){
  88     AliDebug(3,"No part candidate");
  89     return bSignAssoc;
  90   }
  91
  92
  93   bSignAssoc = kTRUE;
  94   return bSignAssoc;
  95 }
  96
  97 //______________________________________________
  98 Bool_t AliCFVertexingHF2Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
  99         //
 100         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 101         //
 102
 103         Bool_t bGenValues = kFALSE;
 104
 105         Double_t vtx1[3] = {0,0,0};   // primary vertex
 106         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 107         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 108         fmcPartCandidate->XvYvZv(vtx1);  // cm
 109
 110
 111         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 112         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 113         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 114         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 115
 116         // getting vertex from daughters
 117         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 118         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 119         if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
 120           AliError("Daughters have different secondary vertex, skipping the track");
 121           return bGenValues;
 122         }
 123
 124         Int_t nprongs = 2;
 125         Short_t charge = 0;
 126         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
 127         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
 128         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
 129         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
 130                 // inverting in case the positive daughter is the second one
 131                 positiveDaugh = mcPartDaughter1;
 132                 negativeDaugh = mcPartDaughter0;
 133         }
 134         // getting the momentum from the daughters
 135         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
 136         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
 137         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
 138
 139         Double_t d0[2] = {0.,0.};
 140
 141         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
 142
 143         Double_t cosThetaStar = 0.;
 144         Double_t cosThetaStarD0 = 0.;
 145         Double_t cosThetaStarD0bar = 0.;
 146         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
 147         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
 148         if (fmcPartCandidate->GetPdgCode() == 421){  // D0
 149           AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 150           cosThetaStar = cosThetaStarD0;
 151         }
 152         else if (fmcPartCandidate->GetPdgCode() == -421){  // D0bar{
 153           AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 154           cosThetaStar = cosThetaStarD0bar;
 155         }
 156         else{
 157           AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
 158           return vectorMC;
 159         }
 160         if (cosThetaStar < -1 || cosThetaStar > 1) {
 161           AliWarning("Invalid value for cosine Theta star, returning");
 162           return bGenValues;
 163         }
 164
 165         //ct
 166         Double_t cT = decay->Ct(421);
 167         // get the pT of the daughters
 168         Double_t pTpi = 0.;
 169         Double_t pTK = 0.;
 170
 171         if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
 172                 pTpi = mcPartDaughter0->Pt();
 173                 pTK = mcPartDaughter1->Pt();
 174         }
 175         else {
 176                 pTpi = mcPartDaughter1->Pt();
 177                 pTK = mcPartDaughter0->Pt();
 178         }
 179
 180         vectorMC[0] = fmcPartCandidate->Pt();
 181         vectorMC[1] = fmcPartCandidate->Y() ;
 182         vectorMC[2] = cosThetaStar ;
 183         vectorMC[3] = pTpi ;
 184         vectorMC[4] = pTK ;
 185         vectorMC[5] = cT*1.E4 ;  // in micron
 186         vectorMC[6] = 0.;   // dummy value, meaningless in MC
 187         vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
 188         vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
 189         vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
 190         vectorMC[10] = 1.01;    // dummy value, meaningless in MC
 191         vectorMC[11] = fmcPartCandidate->Phi();
 192         vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
 193         bGenValues = kTRUE;
 194         return bGenValues;
 195 }
 196
 197
 198
 199 //____________________________________________
 200 Bool_t AliCFVertexingHF2Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 201 {
 202   Bool_t bFillRecoValues=kFALSE;
 203
 204   AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)fRecoCandidate;
 205
 206   if (d0toKpi->GetPrimaryVtx())printf("d0toKpi has primary vtx\n");
 207   if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidate has primary vtx\n");
 208
 209   Double_t pt = d0toKpi->Pt();
 210   Double_t rapidity = d0toKpi->YD0();
 211   Double_t invMass=0.;
 212   Double_t cosThetaStar = 9999.;
 213   Double_t pTpi = 0.;
 214   Double_t pTK = 0.;
 215   Double_t dca = d0toKpi->GetDCA();
 216   Double_t d0pi = 0.;
 217   Double_t d0K = 0.;
 218   Double_t d0xd0 = d0toKpi->Prodd0d0();
 219   Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
 220   Double_t phi = d0toKpi->Phi();
 221   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 222
 223
 224   if (pdgCode > 0){
 225     cosThetaStar = d0toKpi->CosThetaStarD0();
 226     pTpi = d0toKpi->PtProng(0);
 227     pTK = d0toKpi->PtProng(1);
 228     d0pi = d0toKpi->Getd0Prong(0);
 229     d0K = d0toKpi->Getd0Prong(1);
 230     invMass=d0toKpi->InvMassD0();
 231   }
 232   else {
 233     cosThetaStar = d0toKpi->CosThetaStarD0bar();
 234     pTpi = d0toKpi->PtProng(1);
 235     pTK = d0toKpi->PtProng(0);
 236     d0pi = d0toKpi->Getd0Prong(1);
 237     d0K = d0toKpi->Getd0Prong(0);
 238     invMass= d0toKpi->InvMassD0bar();
 239   }
 240
 241   Double_t cT = d0toKpi->CtD0();
 242
 243
 244   vectorReco[0] = pt;
 245   vectorReco[1] = rapidity;
 246   vectorReco[2] = cosThetaStar;
 247   vectorReco[3] = pTpi;
 248   vectorReco[4] = pTK;
 249   vectorReco[5] = cT*1.E4;  // in micron
 250   vectorReco[6] = dca*1.E4;  // in micron
 251   vectorReco[7] = d0pi*1.E4;  // in micron
 252   vectorReco[8] = d0K*1.E4;  // in micron
 253   vectorReco[9] = d0xd0*1.E8;  // in micron^2
 254   vectorReco[10] = cosPointingAngle;  // in micron
 255   vectorReco[11] = phi;
 256   vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
 257
 258   bFillRecoValues = kTRUE;
 259
 260
 261   return bFillRecoValues;
 262 }
 263
 264
 265 //_____________________________________________________________
 266 Bool_t AliCFVertexingHF2Prong::CheckMCChannelDecay() const
 267 {
 268   Bool_t checkCD = kFALSE;
 269   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 270   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 271   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 272   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 273
 274   if (!mcPartDaughter0 || !mcPartDaughter1) {
 275     AliDebug (2,"Problems in the MC Daughters\n");
 276     return checkCD;
 277   }
 278
 279   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
 280         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 281       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 282         TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
 283     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 284     return checkCD;
 285   }
 286
 287   checkCD = kTRUE;
 288   return checkCD;
 289
 290 }
 291
 292 //_______________________________________________