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         // constructor
  43         //
  44
  45         SetNProngs(2);
  46 }
  47
  48
  49 //_____________________________________
  50 AliCFVertexingHF2Prong& AliCFVertexingHF2Prong::operator=(const AliCFVertexingHF2Prong& c)
  51 {
  52         //
  53         // copy constructor
  54         //
  55
  56         if  (this != &c) {
  57                 AliCFVertexingHF::operator=(c);
  58         }
  59         return *this;
  60 }
  61
  62 //__________________________________________
  63 Bool_t AliCFVertexingHF2Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
  64 {
  65         //
  66         // setting the recontructed candidate
  67         //
  68
  69         Bool_t bSignAssoc = kFALSE;
  70         fRecoCandidate = recoCand;
  71         if (!fRecoCandidate) {
  72                 AliError("fRecoCandidate not found, problem in assignement\n");
  73                 return bSignAssoc;
  74         }
  75
  76         if (fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
  77         if (recoCand->GetPrimaryVtx()) AliDebug(3,"Reco Cand has a pointer to PrimVtx\n");
  78
  79         Int_t pdgCand = 421;
  80         Int_t pdgDgD0toKpi[2]={321,211};
  81         Int_t nentries = fmcArray->GetEntriesFast();
  82
  83         AliDebug(3,Form("nentries = %d\n", nentries));
  84
  85         Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,2,pdgDgD0toKpi);
  86         if (mcLabel == -1) return bSignAssoc;
  87         SetMCLabel(mcLabel);
  88         fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
  89         if (!fmcPartCandidate){
  90                 AliDebug(3,"No part candidate");
  91                 return bSignAssoc;
  92         }
  93         bSignAssoc = kTRUE;
  94         return bSignAssoc;
  95 }
  96
  97 //______________________________________________
  98 Bool_t AliCFVertexingHF2Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
  99 {
 100         //
 101         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 102         //
 103
 104         Bool_t bGenValues = kFALSE;
 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         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 111         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 112         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 113         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 114         if(!mcPartDaughter0 || !mcPartDaughter1) return bGenValues;
 115
 116         // getting vertex from daughters
 117         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 118         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 119         if (TMath::Abs(vtx2daughter0[0] - vtx2daughter1[0]) > 1E-5 || TMath::Abs(vtx2daughter0[1] - vtx2daughter1[1]) > 1E-5 || TMath::Abs(vtx2daughter0[2] - vtx2daughter1[2])>1E-5) {
 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         vectorMC[13] = fCentValue;   //reconstructed centrality
 194
 195         delete decay;
 196         bGenValues = kTRUE;
 197         return bGenValues;
 198 }
 199 //____________________________________________
 200 Bool_t AliCFVertexingHF2Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 201 {
 202         //
 203         // Getting the reconstructed values from the candidate
 204         //
 205
 206         Bool_t bFillRecoValues=kFALSE;
 207
 208         AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)fRecoCandidate;
 209
 210         if (d0toKpi->GetPrimaryVtx())AliDebug(3,"d0toKpi has primary vtx\n");
 211         if (fRecoCandidate->GetPrimaryVtx())AliDebug(3,"fRecoCandidate has primary vtx\n");
 212
 213         Double_t pt = d0toKpi->Pt();
 214         Double_t rapidity = d0toKpi->YD0();
 215         Double_t invMass=0.;
 216         Double_t cosThetaStar = 9999.;
 217         Double_t pTpi = 0.;
 218         Double_t pTK = 0.;
 219         Double_t dca = d0toKpi->GetDCA();
 220         Double_t d0pi = 0.;
 221         Double_t d0K = 0.;
 222         Double_t d0xd0 = d0toKpi->Prodd0d0();
 223         Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
 224         Double_t phi = d0toKpi->Phi();
 225         Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 226
 227         if (pdgCode > 0){
 228                 cosThetaStar = d0toKpi->CosThetaStarD0();
 229                 pTpi = d0toKpi->PtProng(0);
 230                 pTK = d0toKpi->PtProng(1);
 231                 d0pi = d0toKpi->Getd0Prong(0);
 232                 d0K = d0toKpi->Getd0Prong(1);
 233                 invMass=d0toKpi->InvMassD0();
 234         }
 235         else {
 236                 cosThetaStar = d0toKpi->CosThetaStarD0bar();
 237                 pTpi = d0toKpi->PtProng(1);
 238                 pTK = d0toKpi->PtProng(0);
 239                 d0pi = d0toKpi->Getd0Prong(1);
 240                 d0K = d0toKpi->Getd0Prong(0);
 241                 invMass= d0toKpi->InvMassD0bar();
 242         }
 243
 244         Double_t cT = d0toKpi->CtD0();
 245
 246         vectorReco[0] = pt;
 247         vectorReco[1] = rapidity;
 248         vectorReco[2] = cosThetaStar;
 249         vectorReco[3] = pTpi;
 250         vectorReco[4] = pTK;
 251         vectorReco[5] = cT*1.E4;  // in micron
 252         vectorReco[6] = dca*1.E4;  // in micron
 253         vectorReco[7] = d0pi*1.E4;  // in micron
 254         vectorReco[8] = d0K*1.E4;  // in micron
 255         vectorReco[9] = d0xd0*1.E8;  // in micron^2
 256         vectorReco[10] = cosPointingAngle;  // in micron
 257         vectorReco[11] = phi;
 258         vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
 259         vectorReco[13] = fCentValue; //reconstructed centrality
 260
 261         bFillRecoValues = kTRUE;
 262
 263         return bFillRecoValues;
 264 }
 265
 266 //_____________________________________________________________
 267 Bool_t AliCFVertexingHF2Prong::CheckMCChannelDecay() const
 268 {
 269         //
 270         // checking the MC decay channel
 271         //
 272         Bool_t checkCD = kFALSE;
 273         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 274         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 275         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 276         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 277
 278         if (!mcPartDaughter0 || !mcPartDaughter1) {
 279                 AliDebug (2,"Problems in the MC Daughters\n");
 280                 return checkCD;
 281         }
 282
 283         if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
 284               TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 285             !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 286               TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
 287                 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 288                 return checkCD;
 289         }
 290
 291         checkCD = kTRUE;
 292         return checkCD;
 293
 294 }
 295