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