PWGHF/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 #include "AliCFTaskVertexingHF.h"
  34 #include "AliPIDResponse.h"
  35 #include "AliPID.h"
  36
  37 ClassImp(AliCFVertexingHFCascade)
  38
  39
  40 //_________________________________________
  41 AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
  42 AliCFVertexingHF(mcArray, originDselection),
  43   fPDGcascade(0),
  44   fPDGbachelor(0),
  45   fPDGneutrDaugh(0),
  46   fPDGneutrDaughForMC(0),
  47   fPDGneutrDaughPositive(0),
  48   fPDGneutrDaughNegative(0),
  49   fPrimVtx(0x0),
  50   fUseCutsForTMVA(kFALSE)
  51 {
  52   // standard constructor
  53
  54   SetNProngs(3);
  55   fPtAccCut = new Float_t[fProngs];
  56   fEtaAccCut = new Float_t[fProngs];
  57   // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
  58   fPtAccCut[0] = 0.;
  59   fEtaAccCut[0] = 0.;
  60   for(Int_t iP=1; iP<fProngs; iP++){
  61     fPtAccCut[iP] = 0.1;
  62     fEtaAccCut[iP] = 0.9;
  63   }
  64
  65 }
  66
  67
  68 //_____________________________________
  69 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
  70 {
  71   // operator =
  72
  73   if  (this != &c) {
  74
  75     AliCFVertexingHF::operator=(c);
  76
  77   }
  78   return *this;
  79 }
  80
  81 //__________________________________________
  82 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
  83 {
  84   // set the AliAODRecoDecay candidate
  85
  86   Bool_t bSignAssoc = kFALSE;
  87
  88   fRecoCandidate = recoCand;
  89   AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
  90
  91   if (!fRecoCandidate) {
  92     AliError("fRecoCandidate not found, problem in assignement\n");
  93     return bSignAssoc;
  94   }
  95
  96   if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
  97
  98   //Int_t pdgCand = 413;
  99
 100   Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
 101   Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
 102
 103   Int_t nentries = fmcArray->GetEntriesFast();
 104
 105   AliDebug(3,Form("nentries = %d\n", nentries));
 106
 107   Bool_t isV0 = kFALSE;
 108   if (fPDGcascade == 4122) {
 109     isV0 = kTRUE;
 110     pdgDgCascade[0] = fPDGbachelor;
 111     pdgDgCascade[1] = fPDGneutrDaugh;
 112   }
 113   AliDebug(3, Form("calling MatchToMC with: fPDGcascade = %d, fPDGneutrDaugh = %d, pdgDgCascade[0] = %d, pdgDgCascade[1] = %d, pdgDgNeutrDaugh[0] = %d, pdgDgNeutrDaugh[1] = %d, fmcArray = %p", fPDGcascade, fPDGneutrDaugh, pdgDgCascade[0], pdgDgCascade[1], pdgDgNeutrDaugh[0], pdgDgNeutrDaugh[1], fmcArray));
 114  Int_t mcLabel = cascade->MatchToMC(fPDGcascade, fPDGneutrDaugh, pdgDgCascade, pdgDgNeutrDaugh, fmcArray, isV0);
 115
 116   if (mcLabel == -1) return bSignAssoc;
 117
 118   if (fRecoCandidate->NumberOfFakeDaughters()>0){
 119     fFake = 0;    // fake candidate
 120     if (fFakeSelection == 1) return bSignAssoc;
 121   }
 122   if (fRecoCandidate->NumberOfFakeDaughters()==0){
 123     fFake = 2;    // non-fake candidate
 124     if (fFakeSelection == 2) return bSignAssoc;
 125   }
 126
 127   SetMCLabel(mcLabel);
 128   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 129
 130   if (!fmcPartCandidate){
 131     AliDebug(3,"No part candidate");
 132     return bSignAssoc;
 133   }
 134
 135   bSignAssoc = kTRUE;
 136   return bSignAssoc;
 137 }
 138
 139 //______________________________________________
 140 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
 141 {
 142   //
 143   // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 144   //
 145
 146   Bool_t bGenValues = kFALSE;
 147
 148   Int_t daughter0cascade = fmcPartCandidate->GetDaughter(0);
 149   Int_t daughter1cascade = fmcPartCandidate->GetDaughter(1);
 150
 151   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0cascade));
 152   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1cascade));
 153   AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
 154   AliAODMCParticle* mcPartDaughterBachelor = NULL;
 155
 156   // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
 157   // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
 158   // charge of the daughters to decide which is which
 159   if (mcPartDaughter0->Charge()/3 == 0){
 160     mcPartDaughterNeutrDaugh = mcPartDaughter0;
 161     mcPartDaughterBachelor =  mcPartDaughter1;
 162   }
 163   else {
 164     mcPartDaughterNeutrDaugh = mcPartDaughter1;
 165     mcPartDaughterBachelor =  mcPartDaughter0;
 166   }
 167
 168   if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
 169
 170   Double_t vtx1[3] = {0,0,0};   // primary vertex
 171   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 172   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 173   fmcPartCandidate->XvYvZv(vtx1);  // cm
 174
 175   //Daughters of the neutral particle of the cascade
 176   Int_t daughter0 = mcPartDaughterNeutrDaugh->GetDaughter(0); // this is the positive
 177   Int_t daughter1 = mcPartDaughterNeutrDaugh->GetDaughter(1); // this is the negative
 178
 179   AliAODMCParticle* mcPartNeutrDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 180   AliAODMCParticle* mcPartNeutrDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 181
 182   if (!mcPartNeutrDaughter0 || !mcPartNeutrDaughter1) return kFALSE;
 183
 184   // getting vertex from daughters
 185   mcPartNeutrDaughter0->XvYvZv(vtx2daughter0);  // cm
 186   mcPartNeutrDaughter1->XvYvZv(vtx2daughter1);  //cm
 187   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) {
 188     AliError("Daughters have different secondary vertex, skipping the track");
 189     return bGenValues;
 190   }
 191
 192   Int_t nprongs = 2;
 193   Short_t charge = 0;
 194   // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
 195   AliAODMCParticle* positiveDaugh = mcPartNeutrDaughter0;
 196   AliAODMCParticle* negativeDaugh = mcPartNeutrDaughter1;
 197   if (mcPartNeutrDaughter0->GetPdgCode() < 0 && mcPartNeutrDaughter1->GetPdgCode() > 0){
 198     // inverting in case the positive daughter is the second one
 199     positiveDaugh = mcPartNeutrDaughter1;
 200     negativeDaugh = mcPartNeutrDaughter0;
 201   }
 202
 203   // getting the momentum from the daughters
 204   Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
 205   Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
 206   Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
 207
 208   Double_t d0[2] = {0.,0.};
 209
 210   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1, vtx2daughter0, nprongs, charge, px, py, pz, d0);
 211
 212   Double_t cosThetaStar = 0.;
 213   Double_t cosThetaStarNeutrDaugh = 0.;
 214   Double_t cosThetaStarNeutrDaughBar = 0.;
 215   cosThetaStarNeutrDaugh = decay->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
 216   cosThetaStarNeutrDaughBar = decay->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughNegative, fPDGneutrDaughPositive);
 217   if (mcPartDaughterNeutrDaugh->GetPdgCode() == fPDGneutrDaughForMC){  // neutral particle
 218     AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
 219     cosThetaStar = cosThetaStarNeutrDaugh;
 220   }
 221   else if (mcPartDaughterNeutrDaugh->GetPdgCode() == -fPDGneutrDaughForMC){  // neutral particle bar
 222     AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 223     cosThetaStar = cosThetaStarNeutrDaughBar;
 224   }
 225   else{
 226     AliWarning(Form("There are problems!! particle was expected to be either with pdg = %d or its antiparticle with pdg = %d, while we have a %d, check...", fPDGneutrDaughForMC, -fPDGneutrDaughForMC, mcPartDaughterNeutrDaugh->GetPdgCode()));
 227     delete decay;
 228     return vectorMC;
 229   }
 230   if (cosThetaStar < -1 || cosThetaStar > 1) {
 231     AliWarning(Form("Invalid value for cosine Theta star %f, returning", cosThetaStar));
 232     delete decay;
 233     return bGenValues;
 234   }
 235
 236   Double_t vectorNeutrDaugh[2] = {0.,0.};
 237
 238   // evaluate the correct cascade
 239   if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
 240     AliDebug(2, "Error! the Neutral Daughter MC doesn't have correct daughters!!");
 241     delete decay;
 242     return bGenValues;
 243   }
 244
 245   //ct
 246   Double_t cT = decay->Ct(fPDGneutrDaugh);
 247   // // get the pT of the daughters
 248   // Double_t pTNeutrDaugh= 0.;
 249   // Double_t pTBachelor = 0.;
 250
 251   // if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
 252   //   pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
 253   //   pTBachelor = mcPartDaughterBachelor->Pt();
 254   // }
 255
 256   AliDebug(3, Form("The candidate has pt = %f, y = %f", fmcPartCandidate->Pt(), fmcPartCandidate->Y()));
 257
 258   switch (fConfiguration){
 259   case AliCFTaskVertexingHF::kSnail:
 260     vectorMC[0] = fmcPartCandidate->Pt();
 261     vectorMC[1] = fmcPartCandidate->Y() ;
 262     vectorMC[2] = cosThetaStar ;
 263     vectorMC[3] = vectorNeutrDaugh[0];
 264     vectorMC[4] = vectorNeutrDaugh[1];
 265     vectorMC[5] = cT*1.E4 ;  // in micron
 266     vectorMC[6] = 0.;   // dummy value, meaningless in MC
 267     vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
 268     vectorMC[8] = 1.01;    // dummy value, meaningless in MC
 269     vectorMC[9] = fmcPartCandidate->Phi();
 270     vectorMC[10] = fzMCVertex;    // z of reconstructed of primary vertex
 271     vectorMC[11] = fCentValue; // reconstructed centrality
 272     vectorMC[12] = 1.;           // always filling with 1 at MC level
 273     vectorMC[13] = 1.01; // dummy value for cosPointingXY  multiplicity
 274     vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
 275     vectorMC[15] = fMultiplicity; // reconstructed multiplicity
 276     break;
 277   case AliCFTaskVertexingHF::kCheetah:
 278     vectorMC[0] = fmcPartCandidate->Pt();
 279     vectorMC[1] = fmcPartCandidate->Y() ;
 280     vectorMC[2] = cT*1.E4; // in micron
 281     vectorMC[3] = fmcPartCandidate->Phi();
 282     vectorMC[4] = fzMCVertex;
 283     vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
 284     vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level
 285     vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
 286     break;
 287   }
 288
 289   delete decay;
 290   bGenValues = kTRUE;
 291   return bGenValues;
 292 }
 293 //____________________________________________
 294 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 295 {
 296   // read the variables for the container
 297
 298   Bool_t bFillRecoValues = kFALSE;
 299
 300   //Get the cascade and the neutral particle from it
 301   AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
 302   AliAODRecoDecay* neutrDaugh = NULL;
 303   if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
 304   else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
 305   else {
 306     return kFALSE;
 307   }
 308
 309   //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
 310   //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
 311
 312   Double_t pt =  cascade->Pt();
 313   Double_t rapidity =  cascade->Y(fPDGcascade);
 314   // Double_t invMass = 0.;
 315   Double_t cosThetaStar = 9999.;
 316   Double_t pTneutrDaughPos = 0.;
 317   Double_t pTneutrDaughNeg = 0.;
 318   Double_t dca = neutrDaugh->GetDCA();
 319   // Double_t d0neutrDaughPos = 0.;
 320   // Double_t d0neutrDaughNeg = 0.;
 321   Double_t d0xd0 = neutrDaugh->Prodd0d0();
 322   Double_t cosPointingAngle = neutrDaugh->CosPointingAngle(fPrimVtx);
 323   Double_t phi = cascade->Phi();
 324   Double_t cosPointingAngleXY = neutrDaugh->CosPointingAngleXY(fPrimVtx);
 325   Double_t normDecayLengthXY = neutrDaugh->NormalizedDecayLengthXY(fPrimVtx);
 326
 327   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 328
 329   // UInt_t pdgDaughCascade[2] = { static_cast<UInt_t>(fPDGbachelor),  static_cast<UInt_t>(fPDGneutrDaugh) };    // bachelor is first daughter of cascade
 330   // UInt_t pdgDaughBarCascade[2] = { static_cast<UInt_t>(fPDGneutrDaugh),  static_cast<UInt_t>(fPDGbachelor) }; // bachelor is second daughter in case of a cascade-bar
 331
 332   if (pdgCode > 0){
 333     cosThetaStar = neutrDaugh->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
 334     pTneutrDaughPos = neutrDaugh->PtProng(0);
 335     pTneutrDaughNeg = neutrDaugh->PtProng(1);
 336     // d0neutrDaughPos = neutrDaugh->Getd0Prong(0);
 337     // d0neutrDaughNeg = neutrDaugh->Getd0Prong(1);
 338     // invMass = neutrDaugh->InvMass(2, pdgDaughCascade);
 339   }
 340   else {
 341     cosThetaStar = neutrDaugh->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
 342     pTneutrDaughPos = neutrDaugh->PtProng(1);
 343     pTneutrDaughNeg = neutrDaugh->PtProng(0);
 344     // d0neutrDaughPos = neutrDaugh->Getd0Prong(1);
 345     // d0neutrDaughNeg = neutrDaugh->Getd0Prong(0);
 346     // invMass = neutrDaugh->InvMass(2, pdgDaughBarCascade);
 347   }
 348
 349   Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
 350
 351   switch (fConfiguration){
 352   case AliCFTaskVertexingHF::kSnail:
 353     vectorReco[0] = pt;
 354     vectorReco[1] = rapidity;
 355     vectorReco[2] = cosThetaStar;
 356     vectorReco[3] = pTneutrDaughPos;
 357     vectorReco[4] = pTneutrDaughNeg;
 358     vectorReco[5] = cT*1.E4;  // in micron
 359     vectorReco[6] = dca*1.E4;  // in micron
 360     vectorReco[7] = d0xd0*1.E8;  // in micron^2
 361     vectorReco[8] = cosPointingAngle;  // in micron
 362     vectorReco[9] = phi;
 363     vectorReco[10] = fzPrimVertex;    // z of reconstructed of primary vertex
 364     vectorReco[11] = fCentValue;
 365     vectorReco[12] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
 366     vectorReco[13] = cosPointingAngleXY;
 367     vectorReco[14] = normDecayLengthXY; // in cm
 368     vectorReco[15] = fMultiplicity; // reconstructed multiplicity
 369     break;
 370   case AliCFTaskVertexingHF::kCheetah:
 371     vectorReco[0] = pt;
 372     vectorReco[1] = rapidity;
 373     vectorReco[2] = cT*1.E4; // in micron
 374     vectorReco[3] = phi;
 375     vectorReco[4] = fzPrimVertex;
 376     vectorReco[5] = fCentValue;
 377     vectorReco[6] = fFake;
 378     vectorReco[7] = fMultiplicity;
 379     break;
 380   }
 381
 382   bFillRecoValues = kTRUE;
 383
 384   return bFillRecoValues;
 385 }
 386
 387
 388 //_____________________________________________________________
 389 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
 390 {
 391   // check the required decay channel
 392
 393   Bool_t checkCD = kFALSE;
 394
 395   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 396   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 397   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 398   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 399
 400   if (!mcPartDaughter0 || !mcPartDaughter1) {
 401     AliDebug (2,"Problems in the MC Daughters\n");
 402     return checkCD;
 403   }
 404
 405   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGneutrDaughForMC &&
 406         TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGbachelor) &&
 407       !(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGbachelor &&
 408         TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGneutrDaughForMC)) {
 409     AliDebug(2, Form("The cascade MC doesn't come from a the decay under study, skipping!! (Pdg codes of daughters = %d, %d)", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
 410     return checkCD;
 411   }
 412
 413   // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
 414   AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
 415
 416   // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
 417   // charge of teh daughters to decide which is which
 418   AliDebug(3, Form("Charge0 = %d, Charge1 = %d", mcPartDaughter0->Charge()/3, mcPartDaughter1->Charge()/3));
 419   if (mcPartDaughter0->Charge()/3 != 0){
 420     mcPartDaughterNeutrDaugh = mcPartDaughter1;
 421   }
 422   else {
 423     mcPartDaughterNeutrDaugh = mcPartDaughter0;
 424   }
 425
 426   Double_t vectorNeutrDaugh[2] ={0., 0.};
 427
 428   // We are looking at a cascade ...evaluate the correct cascade
 429   if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
 430     AliDebug(2, "Error! the NeutrDaugh MC doesn't have correct daughters!!");
 431     return checkCD;
 432   }
 433
 434   checkCD = kTRUE;
 435   return checkCD;
 436
 437 }
 438
 439 //__________________________________________
 440 Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
 441 {
 442   //
 443   // check wether D0 is decaing into kpi
 444   //
 445
 446   Bool_t isHadronic = kFALSE;
 447   AliDebug(2, Form("neutralDaugh = %p, pdg = %d", neutralDaugh, neutralDaugh->GetPdgCode()));
 448
 449   if (fPDGcascade == 4122) {
 450     Int_t labelresonanceDaugh = neutralDaugh->GetDaughter(0);
 451     AliAODMCParticle* resonanceDaugh = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelresonanceDaugh));
 452     if (!resonanceDaugh){
 453       return kFALSE;
 454     }
 455     else {
 456       AliDebug(3, Form("The daughter of the resonant particle is a %d (we are looking for a %d)", resonanceDaugh->GetPdgCode(), fPDGneutrDaugh));
 457       if (TMath::Abs(resonanceDaugh->GetPdgCode()) != fPDGneutrDaugh){
 458         return kFALSE;
 459       }
 460       else {
 461         neutralDaugh = resonanceDaugh;
 462       }
 463     }
 464   }
 465
 466   Int_t daughterNeutrDaugh0 = neutralDaugh->GetDaughter(0);
 467   Int_t daughterNeutrDaugh1 = neutralDaugh->GetDaughter(1);
 468
 469   AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
 470   if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
 471     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 472     return isHadronic;
 473   }
 474
 475   Int_t numberOfExpectedDaughters = 2;
 476   if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
 477     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 478     return isHadronic;
 479   }
 480
 481   AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
 482   AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
 483   if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
 484     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
 485     return isHadronic;
 486   }
 487
 488   AliDebug(3, Form("Daughter 0 has pdg = %d, daughter 1 has pdg = %d", mcPartDaughterNeutrDaugh0->GetPdgCode(), mcPartDaughterNeutrDaugh1->GetPdgCode()));
 489   if (!(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughPositive &&
 490         TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughNegative) &&
 491       !(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughNegative &&
 492         TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughPositive)) {
 493     AliDebug(2, "The neutral particle (MC) doesn't come from the required decay, skipping!!");
 494     return isHadronic;
 495   }
 496
 497   Double_t sumPxDau = mcPartDaughterNeutrDaugh0->Px()+mcPartDaughterNeutrDaugh1->Px();
 498   Double_t sumPyDau = mcPartDaughterNeutrDaugh0->Py()+mcPartDaughterNeutrDaugh1->Py();
 499   Double_t sumPzDau = mcPartDaughterNeutrDaugh0->Pz()+mcPartDaughterNeutrDaugh1->Pz();
 500   Double_t pxMother = neutralDaugh->Px();
 501   Double_t pyMother = neutralDaugh->Py();
 502   Double_t pzMother = neutralDaugh->Pz();
 503   AliDebug(3, Form("pxMother = %f, pyMother = %f, pzMother = %f", pxMother, pyMother, pzMother));
 504   AliDebug(3, Form("sumPxDau = %f, sumPyDau = %f, sumPzDau = %f", sumPxDau, sumPyDau, sumPzDau));
 505   if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>0.00001 ||
 506      TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>0.00001 ||
 507      TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>0.00001){
 508     AliDebug(2, "Momentum conservation violated, skipping!!");
 509     return isHadronic;
 510   }
 511
 512   Double_t pTNeutrDaughPositive = 0;
 513   Double_t pTNeutrDaughNegative = 0;
 514
 515   if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
 516     pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
 517     pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
 518   }
 519   else {
 520     pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
 521     pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
 522   }
 523
 524   isHadronic = kTRUE;
 525
 526   vectorNeutrDaugh[0] = pTNeutrDaughPositive;
 527   vectorNeutrDaugh[1] = pTNeutrDaughNegative;
 528
 529   return isHadronic;
 530
 531 }
 532
 533 //___________________________________________________________
 534
 535 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
 536 {
 537   //
 538   // setting the pt cut to be used in the Acceptance steps (MC+Reco)
 539   //
 540
 541   AliDebug(3, "The 3rd element of the pt cut array will correspond to the cut applied to the soft pion - please check that it is correct");
 542   if (fProngs>0){
 543     for (Int_t iP=0; iP<fProngs; iP++){
 544       fPtAccCut[iP]=ptAccCut[iP];
 545     }
 546   }
 547   return;
 548 }
 549
 550
 551
 552 //___________________________________________________________
 553
 554 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
 555 {
 556   //
 557   // setting the eta cut to be used in the Acceptance steps (MC+Reco)
 558   //
 559
 560   AliDebug(3, "The 3rd element of the eta cut array will correspond to the cut applied to the soft pion - please check that it is correct");
 561   if (fProngs>0){
 562     for (Int_t iP=0; iP<fProngs; iP++){
 563       fEtaAccCut[iP] = etaAccCut[iP];
 564     }
 565   }
 566   return;
 567 }
 568 //___________________________________________________________
 569
 570 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 571 {
 572   //
 573   // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 574   //
 575
 576   AliDebug(3, "The 3rd element of the pt and cut array will correspond to the cut applied to the soft pion - please check that they are correct");
 577   if (fProngs>0){
 578     for (Int_t iP=0; iP<fProngs; iP++){
 579       fPtAccCut[iP]=ptAccCut[iP];
 580       fEtaAccCut[iP]=etaAccCut[iP];
 581     }
 582   }
 583   return;
 584 }
 585
 586 //___________________________________________________________
 587
 588 void AliCFVertexingHFCascade::SetAccCut()
 589 {
 590   //
 591   // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 592   //
 593
 594   Int_t bachelorPosition = 2;
 595   if (fPDGcascade == 4122) bachelorPosition = 0;
 596   AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[bachelorPosition]));  // should be the soft pion...
 597   if(!mcPartDaughter) return;
 598   Int_t mother =  mcPartDaughter->GetMother();
 599   AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
 600   if(!mcMother) return;
 601
 602   if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
 603     AliFatal("Apparently the soft pion is not in the third position, causing a crash!!");
 604   }
 605   if (fProngs>0){
 606     for (Int_t iP=0; iP<fProngs-1; iP++){
 607       fPtAccCut[iP]=0.1;
 608       fEtaAccCut[iP]=0.9;
 609     }
 610     fPtAccCut[2]=0.06;  // soft pion
 611     fEtaAccCut[2]=0.9;  // soft pion
 612   }
 613   return;
 614 }
 615
 616 //_____________________________________________________________
 617 Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const
 618 {
 619   //
 620   // getting eta of the prong - overload the mother class method
 621   //
 622
 623   if (fRecoCandidate){
 624
 625     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
 626
 627     Double_t etaProng =-9999;
 628     AliAODRecoDecay* neutrDaugh=0;
 629     if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
 630     else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
 631     if (iProng==0) etaProng = neutrDaugh->EtaProng(0);
 632     if (iProng==1) etaProng = neutrDaugh->EtaProng(1);
 633     if (iProng==2) etaProng = cascade->EtaProng(1);
 634
 635     return etaProng;
 636
 637   }
 638   return 999999;
 639 }
 640 //_____________________________________________________________
 641 Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const
 642 {
 643   //
 644   // getting pt of the prong
 645   //
 646
 647   if (fRecoCandidate){
 648
 649     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
 650     Double_t ptProng= -9999;
 651     AliAODRecoDecay* neutrDaugh=0;
 652     if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
 653     else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
 654     if (iProng == 0) ptProng = neutrDaugh->PtProng(0);
 655     if (iProng == 1) ptProng = neutrDaugh->PtProng(1);
 656     if (iProng == 2) ptProng = cascade->PtProng(1);
 657
 658     //  Double_t ptProng = fRecoCandidate->PtProng(iProng);
 659     return ptProng;
 660
 661   }
 662   return 999999;
 663
 664 }
 665 //_____________________________________________________________
 666 Bool_t AliCFVertexingHFCascade::CheckAdditionalCuts(AliPIDResponse* pidResponse) const {
 667
 668   // function to check whether the candidate passes the additional cuts defined in the task to get the
 669   // invariant mass spectra; these cuts are NOT pt-dependent
 670
 671   if (fPDGcascade == 4122){
 672     // the method is implemented only in this case so far
 673     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
 674     AliAODv0 * v0part = cascade->Getv0();
 675     AliAODTrack *bachelor = (AliAODTrack*)cascade->GetBachelor();
 676
 677     Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
 678     Double_t nSigmaTPCpr=-999.;
 679     Double_t nSigmaTOFpr=-999.;
 680     nSigmaTPCpr = pidResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
 681     nSigmaTOFpr = pidResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
 682     Double_t ptArm = v0part->PtArmV0();
 683     Double_t invmassK0s = v0part->MassK0Short();
 684     Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
 685
 686     Bool_t cutsForTMVA = ((nSigmaTOFpr > -800) && (nSigmaTOFpr < 3)) &&
 687       ((ptArm > 0.07) && (ptArm < 0.105)) &&
 688       ((TMath::Abs(invmassK0s - mK0SPDG)) < 0.01);
 689
 690     if (!fUseCutsForTMVA) cutsForTMVA = kTRUE;
 691
 692     Bool_t cutsForInvMassTask = !(onFlyV0) &&
 693       (cascade->CosV0PointingAngle()>0.99) &&
 694       (TMath::Abs(nSigmaTPCpr) <= 3) &&
 695       (v0part->Getd0Prong(0) < 20) &&
 696       (v0part->Getd0Prong(1) < 20);
 697
 698     if (cutsForTMVA && cutsForInvMassTask) {
 699       // K0Smass cut
 700       // eta cut
 701       // TOF PID cut
 702       // Arm cut
 703       return kTRUE;
 704     }
 705   }
 706
 707   return kFALSE;
 708
 709 }