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         fPtAccCut=new Float_t[fProngs];
  47         fEtaAccCut=new Float_t[fProngs];
  48         for(Int_t iP=0; iP<fProngs; iP++){
  49                 fPtAccCut[iP]=0.1;
  50                 fEtaAccCut[iP]=0.9;
  51         }
  52
  53 }
  54
  55
  56 //_____________________________________
  57 AliCFVertexingHF2Prong& AliCFVertexingHF2Prong::operator=(const AliCFVertexingHF2Prong& c)
  58 {
  59         //
  60         // copy constructor
  61         //
  62
  63         if  (this != &c) {
  64                 AliCFVertexingHF::operator=(c);
  65         }
  66         return *this;
  67 }
  68
  69 //__________________________________________
  70 Bool_t AliCFVertexingHF2Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
  71 {
  72         //
  73         // setting the recontructed candidate
  74         //
  75
  76         Bool_t bSignAssoc = kFALSE;
  77         fRecoCandidate = recoCand;
  78         if (!fRecoCandidate) {
  79                 AliError("fRecoCandidate not found, problem in assignement\n");
  80                 return bSignAssoc;
  81         }
  82
  83         if (fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
  84         if (recoCand->GetPrimaryVtx()) AliDebug(3,"Reco Cand has a pointer to PrimVtx\n");
  85
  86         Int_t pdgCand = 421;
  87         Int_t pdgDgD0toKpi[2]={321,211};
  88         Int_t nentries = fmcArray->GetEntriesFast();
  89
  90         AliDebug(3,Form("nentries = %d\n", nentries));
  91
  92         Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,2,pdgDgD0toKpi);
  93         if (mcLabel == -1) return bSignAssoc;
  94
  95         if (fRecoCandidate->NumberOfFakeDaughters()>0){
  96                 fFake = 0;    // fake candidate
  97                 if (fFakeSelection==1) return bSignAssoc;
  98         }
  99         if (fRecoCandidate->NumberOfFakeDaughters()==0){
 100                 fFake = 2;    // non-fake candidate
 101                 if (fFakeSelection==2) return bSignAssoc;
 102         }
 103
 104         SetMCLabel(mcLabel);
 105         fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 106         if (!fmcPartCandidate){
 107                 AliDebug(3,"No part candidate");
 108                 return bSignAssoc;
 109         }
 110         bSignAssoc = kTRUE;
 111         return bSignAssoc;
 112 }
 113
 114 //______________________________________________
 115 Bool_t AliCFVertexingHF2Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
 116 {
 117         //
 118         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 119         //
 120
 121         Bool_t bGenValues = kFALSE;
 122         Double_t vtx1[3] = {0,0,0};   // primary vertex
 123         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 124         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 125         fmcPartCandidate->XvYvZv(vtx1);  // cm
 126
 127         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 128         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 129         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 130         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 131         if(!mcPartDaughter0 || !mcPartDaughter1) return bGenValues;
 132
 133         // getting vertex from daughters
 134         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 135         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 136         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) {
 137                 AliError("Daughters have different secondary vertex, skipping the track");
 138                 return bGenValues;
 139         }
 140
 141         Int_t nprongs = 2;
 142         Short_t charge = 0;
 143         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
 144         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
 145         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
 146         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
 147                 // inverting in case the positive daughter is the second one
 148                 positiveDaugh = mcPartDaughter1;
 149                 negativeDaugh = mcPartDaughter0;
 150         }
 151         // getting the momentum from the daughters
 152         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
 153         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
 154         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
 155
 156         Double_t d0[2] = {0.,0.};
 157
 158         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
 159
 160         Double_t cosThetaStar = 0.;
 161         Double_t cosThetaStarD0 = 0.;
 162         Double_t cosThetaStarD0bar = 0.;
 163         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
 164         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
 165         if (fmcPartCandidate->GetPdgCode() == 421){  // D0
 166                 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 167                 cosThetaStar = cosThetaStarD0;
 168         }
 169         else if (fmcPartCandidate->GetPdgCode() == -421){  // D0bar{
 170                 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 171                 cosThetaStar = cosThetaStarD0bar;
 172         }
 173         else{
 174                 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
 175                 return vectorMC;
 176         }
 177         if (cosThetaStar < -1 || cosThetaStar > 1) {
 178                 AliWarning("Invalid value for cosine Theta star, returning");
 179                 return bGenValues;
 180         }
 181
 182         //ct
 183         Double_t cT = decay->Ct(421);
 184         // get the pT of the daughters
 185         Double_t pTpi = 0.;
 186         Double_t pTK = 0.;
 187
 188         if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
 189                 pTpi = mcPartDaughter0->Pt();
 190                 pTK = mcPartDaughter1->Pt();
 191         }
 192         else {
 193                 pTpi = mcPartDaughter1->Pt();
 194                 pTK = mcPartDaughter0->Pt();
 195         }
 196
 197         vectorMC[0] = fmcPartCandidate->Pt();
 198         vectorMC[1] = fmcPartCandidate->Y() ;
 199         vectorMC[2] = cosThetaStar ;
 200         vectorMC[3] = pTpi ;
 201         vectorMC[4] = pTK ;
 202         vectorMC[5] = cT*1.E4 ;  // in micron
 203         vectorMC[6] = 0.;   // dummy value, meaningless in MC
 204         vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
 205         vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
 206         vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
 207         vectorMC[10] = 1.01;    // dummy value, meaningless in MC
 208         vectorMC[11] = fmcPartCandidate->Phi();
 209         vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
 210         vectorMC[13] = fCentValue;   //reconstructed centrality
 211         vectorMC[14] = 1.;           // always filling with 1 at MC level
 212
 213         delete decay;
 214         bGenValues = kTRUE;
 215         return bGenValues;
 216 }
 217 //____________________________________________
 218 Bool_t AliCFVertexingHF2Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 219 {
 220         //
 221         // Getting the reconstructed values from the candidate
 222         //
 223
 224         Bool_t bFillRecoValues=kFALSE;
 225
 226         AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)fRecoCandidate;
 227
 228         if (d0toKpi->GetPrimaryVtx())AliDebug(3,"d0toKpi has primary vtx\n");
 229         if (fRecoCandidate->GetPrimaryVtx())AliDebug(3,"fRecoCandidate has primary vtx\n");
 230
 231         Double_t pt = d0toKpi->Pt();
 232         Double_t rapidity = d0toKpi->YD0();
 233         Double_t invMass=0.;
 234         Double_t cosThetaStar = 9999.;
 235         Double_t pTpi = 0.;
 236         Double_t pTK = 0.;
 237         Double_t dca = d0toKpi->GetDCA();
 238         Double_t d0pi = 0.;
 239         Double_t d0K = 0.;
 240         Double_t d0xd0 = d0toKpi->Prodd0d0();
 241         Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
 242         Double_t phi = d0toKpi->Phi();
 243         Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 244
 245         if (pdgCode > 0){
 246                 cosThetaStar = d0toKpi->CosThetaStarD0();
 247                 pTpi = d0toKpi->PtProng(0);
 248                 pTK = d0toKpi->PtProng(1);
 249                 d0pi = d0toKpi->Getd0Prong(0);
 250                 d0K = d0toKpi->Getd0Prong(1);
 251                 invMass=d0toKpi->InvMassD0();
 252         }
 253         else {
 254                 cosThetaStar = d0toKpi->CosThetaStarD0bar();
 255                 pTpi = d0toKpi->PtProng(1);
 256                 pTK = d0toKpi->PtProng(0);
 257                 d0pi = d0toKpi->Getd0Prong(1);
 258                 d0K = d0toKpi->Getd0Prong(0);
 259                 invMass= d0toKpi->InvMassD0bar();
 260         }
 261
 262         Double_t cT = d0toKpi->CtD0();
 263
 264         vectorReco[0] = pt;
 265         vectorReco[1] = rapidity;
 266         vectorReco[2] = cosThetaStar;
 267         vectorReco[3] = pTpi;
 268         vectorReco[4] = pTK;
 269         vectorReco[5] = cT*1.E4;  // in micron
 270         vectorReco[6] = dca*1.E4;  // in micron
 271         vectorReco[7] = d0pi*1.E4;  // in micron
 272         vectorReco[8] = d0K*1.E4;  // in micron
 273         vectorReco[9] = d0xd0*1.E8;  // in micron^2
 274         vectorReco[10] = cosPointingAngle;  // in micron
 275         vectorReco[11] = phi;
 276         vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
 277         vectorReco[13] = fCentValue; //reconstructed centrality
 278         vectorReco[14] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
 279
 280         bFillRecoValues = kTRUE;
 281
 282         return bFillRecoValues;
 283 }
 284
 285 //_____________________________________________________________
 286 Bool_t AliCFVertexingHF2Prong::CheckMCChannelDecay() const
 287 {
 288         //
 289         // checking the MC decay channel
 290         //
 291         Bool_t checkCD = kFALSE;
 292         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 293         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 294         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 295         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 296
 297         if (!mcPartDaughter0 || !mcPartDaughter1) {
 298                 AliDebug (2,"Problems in the MC Daughters\n");
 299                 return checkCD;
 300         }
 301
 302         if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
 303               TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 304             !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 305               TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
 306                 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 307                 return checkCD;
 308         }
 309
 310         checkCD = kTRUE;
 311         return checkCD;
 312
 313 }
 314