PWG3/vertexingHF/AliCFVertexingHFCascade.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2011, 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 steps
  18 // For D* and other cascades
  19 //
  20 // Author : A.Grelli a.grelli@uu.nl  UTECHT
  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 #include "AliAODRecoCascadeHF.h"
  31 #include "AliCFVertexingHFCascade.h"
  32 #include "AliCFContainer.h"
  33
  34 ClassImp(AliCFVertexingHFCascade)
  35
  36
  37 //_________________________________________
  38   AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
  39   AliCFVertexingHF(mcArray, originDselection)
  40 {
  41   // standard constructor
  42
  43   SetNProngs(3);
  44 }
  45
  46
  47 //_____________________________________
  48 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
  49 {
  50   // operator =
  51
  52   if  (this != &c) {
  53
  54     AliCFVertexingHF::operator=(c);
  55
  56   }
  57   return *this;
  58 }
  59
  60 //__________________________________________
  61 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
  62 {
  63   // set the AliAODRecoDecay candidate
  64
  65   Bool_t bSignAssoc = kFALSE;
  66
  67   fRecoCandidate = recoCand;
  68   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)recoCand;
  69
  70   if (!fRecoCandidate) {
  71     AliError("fRecoCandidate not found, problem in assignement\n");
  72     return bSignAssoc;
  73   }
  74
  75   if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
  76   //if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
  77
  78   //Int_t pdgCand = 413;
  79
  80   Int_t pdgDgDStartoD0pi[2]={421,211};
  81   Int_t pdgDgD0toKpi[2]={321,211};
  82
  83   Int_t nentries = fmcArray->GetEntriesFast();
  84
  85   AliDebug(3,Form("nentries = %d\n", nentries));
  86
  87   Int_t mcLabel =  dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
  88
  89   if (mcLabel < 0) return bSignAssoc;
  90   SetMCLabel(mcLabel);
  91   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
  92
  93   if (!fmcPartCandidate){
  94     AliDebug(3,"No part candidate");
  95     return bSignAssoc;
  96   }
  97
  98   bSignAssoc = kTRUE;
  99   return bSignAssoc;
 100 }
 101
 102 //______________________________________________
 103 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
 104 {
 105   //
 106   // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 107   //
 108
 109         Bool_t bGenValues = kFALSE;
 110
 111
 112         Int_t daughter0ds = fmcPartCandidate->GetDaughter(0);
 113         Int_t daughter1ds = fmcPartCandidate->GetDaughter(1);
 114
 115         //the D0
 116         AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0ds));
 117         AliAODMCParticle* mcPartDaughterPis = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1ds));
 118
 119         Double_t vtx1[3] = {0,0,0};   // primary vertex
 120         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 121         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 122         fmcPartCandidate->XvYvZv(vtx1);  // cm
 123
 124         //D0 daughters
 125         Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
 126         Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);
 127
 128         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
 129         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis
 130
 131         // getting vertex from daughters
 132         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 133         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 134         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) {
 135           AliError("Daughters have different secondary vertex, skipping the track");
 136           return bGenValues;
 137         }
 138
 139         Int_t nprongs = 2;
 140         Short_t charge = 0;
 141         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
 142         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
 143         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
 144         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
 145                 // inverting in case the positive daughter is the second one
 146                 positiveDaugh = mcPartDaughter1;
 147                 negativeDaugh = mcPartDaughter0;
 148         }
 149         // getting the momentum from the daughters
 150         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
 151         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
 152         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
 153
 154         Double_t d0[2] = {0.,0.};
 155
 156         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
 157
 158         Double_t cosThetaStar = 0.;
 159         Double_t cosThetaStarD0 = 0.;
 160         Double_t cosThetaStarD0bar = 0.;
 161         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
 162         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
 163         if (mcPartDaughterD0->GetPdgCode() == 421){  // D0
 164           AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 165           cosThetaStar = cosThetaStarD0;
 166         }
 167         else if (mcPartDaughterD0->GetPdgCode() == -421){  // D0bar{
 168           AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 169           cosThetaStar = cosThetaStarD0bar;
 170         }
 171         else{
 172           AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
 173           return vectorMC;
 174         }
 175         if (cosThetaStar < -1 || cosThetaStar > 1) {
 176           AliWarning("Invalid value for cosine Theta star, returning");
 177           return bGenValues;
 178         }
 179
 180         Double_t vectorD0[2] ={0.,0.};
 181
 182         // evaluate the correct cascade
 183         if (!EvaluateIfD0toKpi(mcPartDaughterD0,vectorD0)) {
 184           AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 185           return bGenValues;
 186         }
 187
 188         //ct
 189         Double_t cT = decay->Ct(421);
 190         // get the pT of the daughters
 191         Double_t pTD0 = 0.;
 192         Double_t pTpi = 0.;
 193
 194         if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
 195                 pTD0 = mcPartDaughterD0->Pt();
 196                 pTpi = mcPartDaughterPis->Pt();
 197         }
 198
 199
 200         vectorMC[0] = fmcPartCandidate->Pt();
 201         vectorMC[1] = fmcPartCandidate->Y() ;
 202         vectorMC[2] = cosThetaStar ;
 203         vectorMC[3] = vectorD0[0];
 204         vectorMC[4] = vectorD0[1];
 205         vectorMC[5] = cT*1.E4 ;  // in micron
 206         vectorMC[6] = 0.;   // dummy value, meaningless in MC
 207         vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
 208         vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
 209         vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
 210         vectorMC[10] = 1.01;    // dummy value, meaningless in MC
 211         vectorMC[11] = fmcPartCandidate->Phi();
 212         vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
 213
 214
 215         delete decay;
 216         bGenValues = kTRUE;
 217         return bGenValues;
 218 }
 219 //____________________________________________
 220 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 221 {
 222   // read the variables for the container
 223
 224   Bool_t bFillRecoValues=kFALSE;
 225
 226   //Get the D* and the D0 from D*
 227   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
 228   AliAODRecoDecayHF2Prong* d0toKpi = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 229
 230
 231   if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi has primary vtx\n");
 232   if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
 233
 234   Double_t pt =  dstarD0pi->Pt();
 235   Double_t rapidity =  dstarD0pi->YDstar();
 236   Double_t invMass=0.;
 237   Double_t cosThetaStar = 9999.;
 238   Double_t pTpi = 0.;
 239   Double_t pTK = 0.;
 240   Double_t dca = d0toKpi->GetDCA();
 241   Double_t d0pi = 0.;
 242   Double_t d0K = 0.;
 243   Double_t d0xd0 = d0toKpi->Prodd0d0();
 244   Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
 245   Double_t phi = dstarD0pi->Phi();
 246
 247   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 248
 249   if (pdgCode > 0){
 250     cosThetaStar = d0toKpi->CosThetaStarD0();
 251     pTpi = d0toKpi->PtProng(0);
 252     pTK = d0toKpi->PtProng(1);
 253     d0pi = d0toKpi->Getd0Prong(0);
 254     d0K = d0toKpi->Getd0Prong(1);
 255     invMass=d0toKpi->InvMassD0();
 256   }
 257   else {
 258     cosThetaStar = d0toKpi->CosThetaStarD0bar();
 259     pTpi = d0toKpi->PtProng(1);
 260     pTK = d0toKpi->PtProng(0);
 261     d0pi = d0toKpi->Getd0Prong(1);
 262     d0K = d0toKpi->Getd0Prong(0);
 263     invMass= d0toKpi->InvMassD0bar();
 264   }
 265
 266   Double_t cT = d0toKpi->CtD0();
 267
 268   vectorReco[0] = pt;
 269   vectorReco[1] = rapidity;
 270   vectorReco[2] = cosThetaStar;
 271   vectorReco[3] = pTpi;
 272   vectorReco[4] = pTK;
 273   vectorReco[5] = cT*1.E4;  // in micron
 274   vectorReco[6] = dca*1.E4;  // in micron
 275   vectorReco[7] = d0pi*1.E4;  // in micron
 276   vectorReco[8] = d0K*1.E4;  // in micron
 277   vectorReco[9] = d0xd0*1.E8;  // in micron^2
 278   vectorReco[10] = cosPointingAngle;  // in micron
 279   vectorReco[11] = phi;
 280   vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
 281
 282   bFillRecoValues = kTRUE;
 283
 284   return bFillRecoValues;
 285 }
 286
 287
 288 //_____________________________________________________________
 289 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
 290 {
 291   // check the required decay channel
 292
 293   Bool_t checkCD = kFALSE;
 294
 295
 296   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 297   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 298   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 299   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 300
 301   if (!mcPartDaughter0 || !mcPartDaughter1) {
 302     AliDebug (2,"Problems in the MC Daughters\n");
 303     return checkCD;
 304   }
 305
 306   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
 307         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 308       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 309         TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
 310     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 311     return checkCD;
 312   }
 313
 314   Double_t vectorD0[2] ={0.,0.};
 315
 316   // D* is a cascade ...evaluate the correct cascade
 317   if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
 318     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 319     return checkCD;
 320   }
 321
 322   checkCD = kTRUE;
 323   return checkCD;
 324
 325 }
 326
 327 //__________________________________________
 328 Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
 329 {
 330   //
 331   // chack wether D0 is decaing into kpi
 332   //
 333
 334   Bool_t isHadronic = kFALSE;
 335
 336   Int_t daughterD00 = neutralDaugh->GetDaughter(0);
 337   Int_t daughterD01 = neutralDaugh->GetDaughter(1);
 338
 339   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
 340   if (daughterD00 == 0 || daughterD01 == 0) {
 341     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 342     return isHadronic;
 343   }
 344
 345   if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
 346     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 347     return isHadronic;
 348   }
 349
 350   AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD00));
 351   AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD01));
 352   if (!mcPartDaughterD00 || !mcPartDaughterD01) {
 353     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
 354     return isHadronic;
 355   }
 356
 357   if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
 358         TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
 359       !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
 360         TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
 361     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 362     return isHadronic;
 363   }
 364
 365   Double_t pTD0pi = 0;
 366   Double_t pTD0K = 0;
 367
 368
 369   if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
 370     pTD0pi = mcPartDaughterD00->Pt();
 371     pTD0K = mcPartDaughterD01->Pt();
 372   }
 373   else {
 374     pTD0pi = mcPartDaughterD01->Pt();
 375     pTD0K  = mcPartDaughterD00->Pt();
 376   }
 377
 378   isHadronic = kTRUE;
 379
 380   vectorD0[0] = pTD0pi;
 381   vectorD0[1] = pTD0K;
 382
 383   return isHadronic;
 384
 385 }
 386
 387
 388
 389