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