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