PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.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 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------
  19 // Class for HF corrections as a function of many variables and steps
  20 // For Lc->V0+bachelor
  21 //
  22 // - Based on AliCFVertexingHFCascade -
  23 //
  24 // Contact : A.De Caro - decaro@sa.infn.it
  25 //           Centro 'E.Fermi' - Rome (Italy)
  26 //           INFN and University of Salerno (Italy)
  27 //
  28 //-----------------------------------------------------------------------
  29
  30 #include "TDatabasePDG.h"
  31 #include "TClonesArray.h"
  32 #include "AliAODv0.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliAODRecoDecayHF.h"
  35 #include "AliAODRecoCascadeHF.h"
  36 #include "AliCFTaskVertexingHF.h"
  37 #include "AliCFContainer.h"
  38 #include "AliCFVertexingHF.h"
  39 #include "AliCFVertexingHFLctoV0bachelor.h"
  40
  41 #include <Riostream.h>
  42
  43 using std::cout;
  44 using std::endl;
  45
  46 ClassImp(AliCFVertexingHFLctoV0bachelor)
  47
  48 //_________________________________________
  49 AliCFVertexingHFLctoV0bachelor::AliCFVertexingHFLctoV0bachelor():
  50 fGenLcOption(0)
  51 {
  52   // standard constructor
  53 }
  54
  55 //_____________________________________
  56 AliCFVertexingHFLctoV0bachelor::AliCFVertexingHFLctoV0bachelor(TClonesArray *mcArray, UShort_t originDselection, Int_t lcDecay):
  57 AliCFVertexingHF(mcArray, originDselection),
  58   fGenLcOption(lcDecay)
  59 {
  60   // standard constructor
  61
  62   SetNProngs(3);
  63   fPtAccCut=new Float_t[fProngs];
  64   fEtaAccCut=new Float_t[fProngs];
  65   for(Int_t iP=0; iP<fProngs; iP++){
  66     fPtAccCut[iP]=0.1;
  67     fEtaAccCut[iP]=0.9;
  68   }
  69
  70 }
  71
  72
  73 //_____________________________________
  74 AliCFVertexingHFLctoV0bachelor& AliCFVertexingHFLctoV0bachelor::operator=(const AliCFVertexingHFLctoV0bachelor& c)
  75 {
  76   // operator =
  77
  78   if  (this != &c) {
  79     AliCFVertexingHF::operator=(c);
  80   }
  81
  82   return *this;
  83
  84 }
  85
  86 //__________________________________________
  87 Bool_t AliCFVertexingHFLctoV0bachelor::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
  88 {
  89   // set the AliAODRecoDecay candidate
  90
  91   Bool_t bSignAssoc = kFALSE;
  92
  93   fRecoCandidate = recoCand;
  94   if (!fRecoCandidate) {
  95     AliError("fRecoCandidate not found, problem in assignement\n");
  96     return bSignAssoc;
  97   }
  98
  99   if (fRecoCandidate->GetPrimaryVtx()) AliDebug(4,"fReco Candidate has a pointer to PrimVtx\n");
 100
 101   AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate;
 102   if ( !(lcV0bachelor->Getv0()) ) {
 103     AliDebug(1,"It is not a Lc->V0+bachelor candidate");
 104     return bSignAssoc;
 105   }
 106
 107   Int_t pdgCand = 4122;
 108   Int_t mcLabel = -1;
 109   Int_t mcLabelK0S = -1;
 110   Int_t mcLabelLambda = -1;
 111
 112   // Lc->K0S+p and cc
 113   Int_t pdgDgLctoV0bachelor[2]={2212,310}; // first bachelor, second V0
 114   Int_t pdgDgV0toDaughters[2]={211,211};
 115   mcLabelK0S = lcV0bachelor->MatchToMC(pdgCand,pdgDgLctoV0bachelor[1],pdgDgLctoV0bachelor,pdgDgV0toDaughters,fmcArray,kTRUE);
 116
 117   // Lc->Lambda+pi and cc
 118   pdgDgLctoV0bachelor[0]=211, pdgDgLctoV0bachelor[1]=3122; // first bachelor, second V0
 119   pdgDgV0toDaughters[0]=2212,  pdgDgV0toDaughters[1]=211;
 120   mcLabelLambda = lcV0bachelor->MatchToMC(pdgCand,pdgDgLctoV0bachelor[1],pdgDgLctoV0bachelor,pdgDgV0toDaughters,fmcArray,kTRUE);
 121
 122   if (mcLabelK0S!=-1 && mcLabelLambda!=-1)
 123     AliDebug(2,"Strange: current Lc->V0+bachelor candidate has two MC different labels!");
 124
 125   if (fGenLcOption==kCountAllLctoV0) {
 126     if (mcLabelK0S!=-1) mcLabel=mcLabelK0S;
 127     if (mcLabelLambda!=-1) mcLabel=mcLabelLambda;
 128   }
 129   else if (fGenLcOption==kCountK0Sp) {
 130     if (mcLabelK0S!=-1) mcLabel=mcLabelK0S;
 131     if (mcLabelLambda!=-1) {
 132       mcLabel=-1;
 133       fFake = 0;    // fake candidate
 134       if (fFakeSelection==1) return bSignAssoc;
 135     }
 136   }
 137   else if (fGenLcOption==kCountLambdapi) {
 138     if (mcLabelLambda!=-1) mcLabel=mcLabelLambda;
 139     if (mcLabelK0S!=-1) {
 140       mcLabel=-1;
 141       fFake = 0;    // fake candidate
 142       if (fFakeSelection==1) return bSignAssoc;
 143     }
 144   }
 145
 146   if (mcLabel==-1) {
 147     AliDebug(4,"No mcLabel found for current candidate");
 148     return bSignAssoc;
 149   }
 150   AliDebug(1,Form("Found mcLabel (%d) for current candidate",mcLabel));
 151
 152   if (fRecoCandidate->NumberOfFakeDaughters()>0){
 153     fFake = 0;    // fake candidate
 154     if (fFakeSelection==1) return bSignAssoc;
 155   }
 156   if (fRecoCandidate->NumberOfFakeDaughters()==0){
 157     fFake = 2;    // non-fake candidate
 158     if (fFakeSelection==2) return bSignAssoc;
 159   }
 160
 161   SetMCLabel(mcLabel); // fmcLabel=mcLabel
 162   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 163   if (!fmcPartCandidate){
 164     AliDebug(3,"No MC object for current candidate");
 165     return bSignAssoc;
 166   }
 167
 168   bSignAssoc = kTRUE;
 169   return bSignAssoc;
 170
 171 }
 172
 173 //______________________________________________
 174 Bool_t AliCFVertexingHFLctoV0bachelor::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
 175 {
 176   //
 177   // collecting all the necessary info (pt, y, invMassV0, cosPAwrtPVxV0, onTheFlyStatusV0) from MC particle
 178   // (additional infos: pTbachelor, pTV0pos, pTV0neg, phi, dcaV0, cTV0, cT, cosPA)
 179   //
 180
 181   Bool_t bGenValues = kFALSE;
 182
 183   if (fmcPartCandidate->GetNDaughters()!=2) {
 184     AliDebug(2,"Lc MC particle doesn't decay in 2 daughters");
 185     return bGenValues;
 186   }
 187
 188   Int_t daughter0lc = fmcPartCandidate->GetDaughter(0);
 189   Int_t daughter1lc = fmcPartCandidate->GetDaughter(1);
 190   if (daughter0lc<0 || daughter1lc<0) {
 191     AliDebug(2,"Lc daughters are not in MC array");
 192     return bGenValues;
 193   }
 194
 195   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0lc));
 196   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1lc));
 197   if (!mcPartDaughter0 || !mcPartDaughter1) {
 198     AliDebug(2,"Problems in the MC Daughters\n");
 199     return bGenValues;
 200   }
 201
 202   if ( fGenLcOption==kCountLambdapi &&
 203        !(TMath::Abs(mcPartDaughter0->GetPdgCode())==3122 &&
 204          TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 205        !(TMath::Abs(mcPartDaughter1->GetPdgCode())==3122 &&
 206          TMath::Abs(mcPartDaughter0->GetPdgCode())==211) ) return bGenValues;
 207   if ( fGenLcOption==kCountK0Sp &&
 208        !(TMath::Abs(mcPartDaughter0->GetPdgCode())==2212 &&
 209          TMath::Abs(mcPartDaughter1->GetPdgCode())==311) &&
 210        !(TMath::Abs(mcPartDaughter1->GetPdgCode())==2212 &&
 211          TMath::Abs(mcPartDaughter0->GetPdgCode())==311) ) return bGenValues;
 212
 213   if ( (TMath::Abs(mcPartDaughter0->GetPdgCode())==311   &&
 214         TMath::Abs(mcPartDaughter1->GetPdgCode())==2212) ||
 215        (TMath::Abs(mcPartDaughter0->GetPdgCode())==3122  &&
 216         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) )
 217     bGenValues = FillVectorFromMCarray(mcPartDaughter1,mcPartDaughter0,vectorMC);
 218   else if ( (TMath::Abs(mcPartDaughter1->GetPdgCode())==311   &&
 219              TMath::Abs(mcPartDaughter0->GetPdgCode())==2212) ||
 220             (TMath::Abs(mcPartDaughter1->GetPdgCode())==3122  &&
 221              TMath::Abs(mcPartDaughter0->GetPdgCode())==211) )
 222     bGenValues = FillVectorFromMCarray(mcPartDaughter0,mcPartDaughter1,vectorMC);
 223
 224   if (!bGenValues)
 225     AliDebug(2,"There is something wrong in filling MC vector");
 226
 227   return bGenValues;
 228
 229 }
 230
 231 //____________________________________________
 232 Bool_t AliCFVertexingHFLctoV0bachelor::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 233 {
 234   // read the variables for the container
 235
 236   Bool_t bFillRecoValues = kFALSE;
 237
 238   //Get the Lc and the V0 from Lc
 239   AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate;
 240
 241   AliAODTrack* bachelor = (AliAODTrack*)lcV0bachelor->GetBachelor();
 242   AliAODv0* v0toDaughters = (AliAODv0*)lcV0bachelor->Getv0();
 243   if (!lcV0bachelor || !bachelor || !v0toDaughters) {
 244     AliDebug(2,"No V0 or bachelor in this reco candidate, skipping!");
 245     return bFillRecoValues;
 246   }
 247
 248   Bool_t onTheFlyStatus = v0toDaughters->GetOnFlyStatus();
 249   AliAODTrack* v0positiveTrack = (AliAODTrack*)lcV0bachelor->Getv0PositiveTrack();
 250   AliAODTrack* v0negativeTrack = (AliAODTrack*)lcV0bachelor->Getv0NegativeTrack();
 251   if (!v0positiveTrack || !v0negativeTrack) {
 252     AliDebug(2,"No V0daughters in this reco candidate, skipping!");
 253     return bFillRecoValues;
 254   }
 255
 256   Double_t pt = lcV0bachelor->Pt();
 257   Double_t rapidity = lcV0bachelor->Y(4122);
 258
 259   Double_t cosPAwrtPrimVtxV0 = lcV0bachelor->CosV0PointingAngle();
 260
 261   Double_t pTbachelor = bachelor->Pt();
 262   Double_t pTV0pos = v0positiveTrack->Pt();
 263   Double_t pTV0neg = v0negativeTrack->Pt();
 264   Double_t phi = lcV0bachelor->Phi();
 265   Double_t dcaV0 = v0toDaughters->GetDCA();
 266   Double_t cTLc = lcV0bachelor->Ct(4122); // wrt PrimVtx
 267   //Double_t dcaLc = lcV0bachelor->GetDCA();
 268   Double_t cosPointingAngleLc = lcV0bachelor->CosPointingAngle();
 269
 270   Double_t cTV0 = 0.;
 271   AliAODVertex *vtx0 = (AliAODVertex*)lcV0bachelor->GetPrimaryVtx();
 272   if (!vtx0) {
 273     AliDebug(2,"Candidate has not primary vtx");
 274   } else {
 275     Double_t primVtxPos[3] = {0.,0.,0.}; vtx0->GetXYZ(primVtxPos);
 276     if (fGenLcOption==kCountK0Sp) {
 277       cTV0 = v0toDaughters->Ct(310,primVtxPos);
 278     } else if (fGenLcOption==kCountLambdapi) {
 279       cTV0 = v0toDaughters->Ct(3122,primVtxPos);
 280     }
 281   }
 282
 283   Double_t invMassV0 = 0.;
 284   if (fGenLcOption==kCountLambdapi) {
 285
 286     Short_t bachelorCharge = bachelor->Charge();
 287     if (bachelorCharge==1) {
 288       invMassV0 = v0toDaughters->MassLambda();
 289     } else if (bachelorCharge==-1) {
 290       invMassV0 = v0toDaughters->MassAntiLambda();
 291     }
 292
 293   } else if (fGenLcOption==kCountK0Sp) {
 294
 295     invMassV0 = v0toDaughters->MassK0Short();
 296
 297   }
 298
 299   vectorReco[0]  = pt;
 300   vectorReco[1]  = rapidity;
 301   vectorReco[2]  = phi;
 302   vectorReco[3]  = cosPAwrtPrimVtxV0;
 303   vectorReco[4]  = onTheFlyStatus;
 304   vectorReco[5]  = fCentValue;
 305   vectorReco[6]  = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
 306   vectorReco[7]  = fMultiplicity;
 307
 308   if (fConfiguration==AliCFTaskVertexingHF::kSnail) {
 309     vectorReco[8]  = pTbachelor;
 310     vectorReco[9]  = pTV0pos;
 311     vectorReco[10] = pTV0neg;
 312     vectorReco[11] = invMassV0;
 313     vectorReco[12] = dcaV0;
 314     vectorReco[13] = cTV0*1.E4; // in micron
 315     vectorReco[14] = cTLc*1.E4; // in micron
 316     vectorReco[15] = cosPointingAngleLc;
 317   }
 318
 319   bFillRecoValues = kTRUE;
 320
 321   return bFillRecoValues;
 322 }
 323
 324 //_____________________________________________________________
 325 Bool_t AliCFVertexingHFLctoV0bachelor::CheckMCChannelDecay() const
 326 {
 327   // check the required decay channel
 328
 329   Bool_t checkCD = kFALSE;
 330
 331   if (fmcPartCandidate->GetNDaughters()!=2) {
 332     AliDebug(2, Form("The MC particle doesn't decay in 2 particles, skipping!!"));
 333     return checkCD;
 334   }
 335
 336   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 337   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 338   if (daughter0<0 || daughter1<0){
 339     AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 340     return checkCD;
 341   }
 342   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 343   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 344   if (!mcPartDaughter0 || !mcPartDaughter1) {
 345     AliDebug(2,"Problems in the MC Daughters\n");
 346     return checkCD;
 347   }
 348
 349   // Lc -> Lambda + pion AND cc
 350   if (fGenLcOption==kCountLambdapi) {
 351
 352     if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==3122 &&
 353           TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 354         !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211  &&
 355           TMath::Abs(mcPartDaughter1->GetPdgCode())==3122)) {
 356       AliDebug(2, "The Lc MC doesn't decay in Lambda+pion (or cc), skipping!!");
 357       return checkCD;
 358     }
 359
 360     if (TMath::Abs(mcPartDaughter0->GetPdgCode())==3122) {
 361       mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); // the bachelor
 362       mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); // the V0
 363     }
 364     if (!mcPartDaughter0 || !mcPartDaughter1) {
 365       AliDebug(2,"Problems in the MC Daughters\n");
 366       return checkCD;
 367     }
 368
 369     if (mcPartDaughter1->GetNDaughters()!=2) {
 370       AliDebug(2, "The Lambda MC particle doesn't decay in 2 particles, skipping!!");
 371       return checkCD;
 372     }
 373
 374     Int_t daughter1D0 = mcPartDaughter1->GetDaughter(0);
 375     Int_t daughter1D1 = mcPartDaughter1->GetDaughter(1);
 376     if (daughter1D0<0 || daughter1D1<0) {
 377       AliDebug(2, Form("The Lambda MC particle doesn't have correct daughters, skipping!!"));
 378       return checkCD;
 379     }
 380
 381     AliAODMCParticle* mcPartDaughter1D0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1D0));
 382     AliAODMCParticle* mcPartDaughter1D1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1D1));
 383     if(!mcPartDaughter1D0 || !mcPartDaughter1D1) {
 384       AliError("The Lambda daughter particle not found in MC array");
 385       return checkCD;
 386     }
 387
 388     if (!(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==211   &&
 389           TMath::Abs(mcPartDaughter1D1->GetPdgCode())==2212) &&
 390         !(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==2212  &&
 391           TMath::Abs(mcPartDaughter1D1->GetPdgCode())==211)) {
 392       AliDebug(2, "The Lambda MC doesn't decay in pi+proton (or cc), skipping!!");
 393       return checkCD;
 394     }
 395
 396   } else if (fGenLcOption==kCountK0Sp) { // Lc -> K0bar + proton AND cc
 397
 398     if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==311   &&
 399           TMath::Abs(mcPartDaughter1->GetPdgCode())==2212) &&
 400         !(TMath::Abs(mcPartDaughter0->GetPdgCode())==2212  &&
 401           TMath::Abs(mcPartDaughter1->GetPdgCode())==311)) {
 402       AliDebug(2, "The Lc MC doesn't decay in K0+proton (or cc), skipping!!");
 403       return checkCD;
 404     }
 405
 406     if (TMath::Abs(mcPartDaughter0->GetPdgCode())==311) {
 407       mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); // the bachelor
 408       mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); // the V0
 409     }
 410     if (!mcPartDaughter0 || !mcPartDaughter1) {
 411       AliDebug(2,"Problems in the MC Daughters after swapping V0 and bachelor\n");
 412       return checkCD;
 413     }
 414
 415     Int_t daughter = mcPartDaughter1->GetDaughter(0);
 416     if (daughter<0) {
 417       AliDebug(2, Form("The K0/K0bar MC particle doesn't have correct daughter, skipping!!"));
 418       return checkCD;
 419     }
 420
 421     AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter));
 422     if(!mcPartDaughter){
 423       AliError("The K0/K0bar daughter particle not found in MC array");
 424       return checkCD;
 425     }
 426
 427     if (!(TMath::Abs(mcPartDaughter->GetPdgCode())==310)) {
 428       AliDebug(2, "The K0/K0bar MC doesn't go in K0S, skipping!!");
 429       return checkCD;
 430     }
 431
 432     if (mcPartDaughter->GetNDaughters()!=2) {
 433       AliDebug(2, "The K0S MC doesn't decay in 2 particles, skipping!!");
 434       return checkCD;
 435     }
 436
 437     Int_t daughterD0 = mcPartDaughter->GetDaughter(0);
 438     Int_t daughterD1 = mcPartDaughter->GetDaughter(1);
 439     if (daughterD0<0 || daughterD1<0) {
 440       AliDebug(2, Form("The K0S MC particle doesn't have correct daughters, skipping!!"));
 441       return checkCD;
 442     }
 443
 444     AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD0));
 445     AliAODMCParticle* mcPartDaughterD1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD1));
 446     if (!mcPartDaughterD0 || !mcPartDaughterD1) {
 447       AliError("Daughter particle not found in MC array");
 448       return checkCD;
 449     }
 450
 451     if (! ( TMath::Abs(mcPartDaughterD0->GetPdgCode())==211 &&
 452             TMath::Abs(mcPartDaughterD1->GetPdgCode())==211 ) ) {
 453       AliDebug(2, "The K0S MC doesn't decay in pi+ pi-, skipping!!");
 454       return checkCD;
 455     }
 456
 457   }
 458
 459   checkCD = kTRUE;
 460   return checkCD;
 461
 462 }
 463
 464 //_____________________________________________________________
 465 Double_t AliCFVertexingHFLctoV0bachelor::GetEtaProng(Int_t iProng) const
 466 {
 467   //
 468   // getting eta of the prong - overload the mother class method
 469   //
 470
 471   Double_t etaProng =-9999;
 472
 473   if (!fRecoCandidate) {
 474     AliDebug(2,"No reco candidate selected");
 475     return etaProng;
 476   }
 477
 478   AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate;
 479   AliAODTrack* bachelor = (AliAODTrack*)lcV0bachelor->GetBachelor();
 480   AliAODTrack* v0Pos = (AliAODTrack*)lcV0bachelor->Getv0PositiveTrack();
 481   AliAODTrack* v0Neg = (AliAODTrack*)lcV0bachelor->Getv0NegativeTrack();
 482   if (!(lcV0bachelor->Getv0()) || !bachelor || !v0Pos || !v0Neg) {
 483     AliDebug(2,"No V0 for this reco candidate selected");
 484     return etaProng;
 485   }
 486
 487   if (iProng==0) etaProng = bachelor->Eta();
 488   else if (iProng==1) etaProng = v0Pos->Eta();
 489   else if (iProng==2) etaProng = v0Neg->Eta();
 490
 491   AliDebug(4,Form("Eta value for prong number %1d = %f",iProng,etaProng));
 492
 493   return etaProng;
 494
 495 }
 496
 497 //_____________________________________________________________
 498
 499 Double_t AliCFVertexingHFLctoV0bachelor::GetPtProng(Int_t iProng) const
 500 {
 501   //
 502   // getting pt of the prong
 503   //
 504
 505   Double_t ptProng=-9999.;
 506
 507   if (!fRecoCandidate) {
 508     AliDebug(2,"No reco candidate selected");
 509     return ptProng;
 510   }
 511
 512   AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate;
 513   AliAODTrack* bachelor = (AliAODTrack*)lcV0bachelor->GetBachelor();
 514   AliAODTrack* v0Pos = (AliAODTrack*)lcV0bachelor->Getv0PositiveTrack();
 515   AliAODTrack* v0Neg = (AliAODTrack*)lcV0bachelor->Getv0NegativeTrack();
 516   if (!(lcV0bachelor->Getv0()) || !bachelor || !v0Pos || !v0Neg) {
 517     AliDebug(2,"No V0 for this reco candidate selected");
 518     return ptProng;
 519   }
 520
 521   if (iProng==0) ptProng = bachelor->Pt();
 522   else if (iProng==1) ptProng = v0Pos->Pt();
 523   else if (iProng==2) ptProng = v0Neg->Pt();
 524
 525   AliDebug(4,Form("Pt value for prong number %1d = %f",iProng,ptProng));
 526
 527   return ptProng;
 528
 529 }
 530
 531 //_____________________________________________________________
 532
 533 Double_t AliCFVertexingHFLctoV0bachelor::Ctau(AliAODMCParticle *mcPartCandidate)
 534 {
 535
 536   Double_t cTau = 999999.;
 537
 538   Int_t daughterD0 = mcPartCandidate->GetDaughter(0);
 539   Int_t daughterD1 = mcPartCandidate->GetDaughter(1);
 540   if (daughterD0<0 || daughterD1<0) {
 541     AliDebug(2, Form("The Lc MC particle doesn't have correct daughters, skipping!!"));
 542     return cTau;
 543   }
 544
 545   AliAODMCParticle *mcPartDaughter0 = (AliAODMCParticle*)fmcArray->At(mcPartCandidate->GetDaughter(0));
 546   AliAODMCParticle *mcPartDaughter1 = (AliAODMCParticle*)fmcArray->At(mcPartCandidate->GetDaughter(1));
 547   if (!mcPartDaughter0 || !mcPartDaughter1) {
 548     AliDebug(2,"The candidate daughter particles not found in MC array");
 549     return cTau;
 550   }
 551
 552   Double_t vtx1[3] = {0,0,0};   // primary vertex
 553   Bool_t hasProdVertex = mcPartCandidate->XvYvZv(vtx1);  // cm
 554
 555   Double_t vtx1daughter[3] = {0,0,0};   // secondary vertex
 556   Bool_t v0Vertex = mcPartDaughter0->XvYvZv(vtx1daughter);  //cm
 557   Double_t vtx2daughter[3] = {0,0,0};   // secondary vertex
 558   Bool_t bachVertex = hasProdVertex && mcPartDaughter1->XvYvZv(vtx2daughter);  //cm
 559
 560   if (!hasProdVertex || !v0Vertex || !bachVertex) {
 561     AliDebug(2,"At least one of Prim.vtx, V0vtx, BachelorVtx doesn't exist!");
 562     return cTau;
 563   }
 564
 565   if (TMath::Abs(vtx1daughter[0]-vtx2daughter[0])>1E-5 ||
 566       TMath::Abs(vtx1daughter[1]-vtx2daughter[1])>1E-5 ||
 567       TMath::Abs(vtx1daughter[2]-vtx2daughter[2])>1E-5) {
 568     AliDebug(2,"Bachelor and V0 haven't common vtx!");
 569     return cTau;
 570   }
 571
 572   Double_t decayLength = 0.;
 573   for (Int_t ii=0; ii<3; ii++) decayLength += (vtx1daughter[ii]-vtx1[ii])*(vtx1daughter[ii]-vtx1[ii]);
 574   decayLength = TMath::Sqrt(decayLength);
 575
 576   cTau = decayLength * mcPartCandidate->M()/mcPartCandidate->P();
 577
 578   AliDebug(2,Form(" cTau(4122)=%f",cTau));
 579
 580   return cTau;
 581
 582 }
 583
 584 //------------
 585 Bool_t AliCFVertexingHFLctoV0bachelor::SetLabelArray()
 586 {
 587   //
 588   // setting the label arrays
 589   //
 590
 591   Bool_t checkCD = kFALSE;
 592
 593   if (fmcPartCandidate->GetNDaughters()!=2) {
 594     AliDebug(2, Form("The MC particle doesn't have 2 daughters, skipping!!"));
 595     return checkCD;
 596   }
 597
 598   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
 599   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
 600   if (daughter0<0 || daughter1<0){
 601     AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 602     return checkCD;
 603   }
 604
 605   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
 606   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
 607   if (!mcPartDaughter0 || !mcPartDaughter1) {
 608     AliDebug(2,"Problems in the MC Daughters\n");
 609     return checkCD;
 610   }
 611
 612
 613   fLabelArray = new Int_t[fProngs];
 614
 615   if (fGenLcOption==kCountLambdapi) { // Lc -> Lambda + pion OR cc
 616
 617     if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==3122 &&
 618           TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 619         !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211  &&
 620           TMath::Abs(mcPartDaughter1->GetPdgCode())==3122)) {
 621       AliDebug(2, "The Lc MC doesn't decay in Lambda+pion (or cc), skipping!!");
 622       delete [] fLabelArray;
 623       fLabelArray = 0x0;
 624       return checkCD;
 625     }
 626
 627     // it is Lc -> Lambda + pion OR cc
 628     if (TMath::Abs(mcPartDaughter0->GetPdgCode())==3122) {
 629       mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); // the bachelor
 630       mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); // the V0
 631       Int_t daughterTemp = daughter0;
 632       daughter0 = daughter1; // the bachelor label
 633       daughter1 = daughterTemp; // the V0 label
 634     }
 635
 636     if (mcPartDaughter1->GetNDaughters()!=2) {
 637       AliDebug(2, "The Lambda MC particle doesn't decay in 2 particles, skipping!!");
 638       delete [] fLabelArray;
 639       fLabelArray = 0x0;
 640       return checkCD;
 641     }
 642
 643     Int_t daughter1D0 = mcPartDaughter1->GetDaughter(0);
 644     Int_t daughter1D1 = mcPartDaughter1->GetDaughter(1);
 645     if (daughter1D0<0 || daughter1D1<0) {
 646       AliDebug(2, Form("The Lambda MC particle doesn't have correct daughters, skipping!!"));
 647       delete [] fLabelArray;
 648       fLabelArray = 0x0;
 649       return checkCD;
 650     }
 651
 652     AliAODMCParticle* mcPartDaughter1D0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1D0));
 653     AliAODMCParticle* mcPartDaughter1D1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1D1));
 654     if (!mcPartDaughter1D0 || !mcPartDaughter1D1) {
 655       AliError("The Lambda daughter particles not found in MC array");
 656       delete [] fLabelArray;
 657       fLabelArray = 0x0;
 658       return checkCD;
 659     }
 660
 661     if (!(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==211   &&
 662           TMath::Abs(mcPartDaughter1D1->GetPdgCode())==2212) &&
 663         !(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==2212  &&
 664           TMath::Abs(mcPartDaughter1D1->GetPdgCode())==211)) {
 665       AliDebug(2, "The Lambda MC doesn't decay in pi+proton (or cc), skipping!!");
 666       delete [] fLabelArray;
 667       fLabelArray = 0x0;
 668       return checkCD;
 669     }
 670
 671     // Lambda -> p+pi OR cc
 672
 673     fLabelArray[0] = daughter0;//mcPartDaughter0->GetLabel(); // bachelor
 674
 675     if (fmcPartCandidate->Charge()>0) {
 676
 677       if (mcPartDaughter1D0->GetPdgCode()==2212) {
 678         fLabelArray[1] = daughter1D0;//mcPartDaughter1D0->GetLabel(); // proton
 679         fLabelArray[2] = daughter1D1;//mcPartDaughter1D1->GetLabel(); // pion
 680       } else if (mcPartDaughter1D1->GetPdgCode()==2212) {
 681         fLabelArray[1] = daughter1D1;//mcPartDaughter1D1->GetLabel(); // proton
 682         fLabelArray[2] = daughter1D0;//mcPartDaughter1D0->GetLabel(); // pion
 683       }
 684
 685     } else if (fmcPartCandidate->Charge()<0) {
 686
 687       if (mcPartDaughter1D0->GetPdgCode()==211) {
 688         fLabelArray[1] = daughter1D0;//mcPartDaughter1D0->GetLabel(); // pion
 689         fLabelArray[2] = daughter1D1;//mcPartDaughter1D1->GetLabel(); // proton
 690       } else if (mcPartDaughter1D1->GetPdgCode()==211) {
 691         fLabelArray[1] = daughter1D1;//mcPartDaughter1D1->GetLabel(); // pion
 692         fLabelArray[2] = daughter1D0;//mcPartDaughter1D0->GetLabel(); // proton
 693       }
 694
 695     }
 696
 697   } else if (fGenLcOption==kCountK0Sp) { // Lc -> K0bar + proton OR cc
 698
 699     if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==311   &&
 700           TMath::Abs(mcPartDaughter1->GetPdgCode())==2212) &&
 701         !(TMath::Abs(mcPartDaughter0->GetPdgCode())==2212  &&
 702           TMath::Abs(mcPartDaughter1->GetPdgCode())==311)) {
 703       AliDebug(2, "The Lc MC doesn't decay in K0bar+proton (or cc), skipping!!");
 704       delete [] fLabelArray;
 705       fLabelArray = 0x0;
 706       return checkCD;
 707     }
 708
 709     if (TMath::Abs(mcPartDaughter0->GetPdgCode())==311) {
 710       mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); // the bachelor
 711       mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); // the V0
 712       Int_t daughterTemp = daughter0;
 713       daughter0 = daughter1; // the bachelor label
 714       daughter1 = daughterTemp; // the V0 label
 715     }
 716
 717     Int_t daughter = mcPartDaughter1->GetDaughter(0);
 718     if (daughter<0) {
 719       AliDebug(2, Form("The K0/K0bar MC particle doesn't have correct daughter, skipping!!"));
 720       delete [] fLabelArray;
 721       fLabelArray = 0x0;
 722       return checkCD;
 723     }
 724
 725     AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter));
 726     if (!mcPartDaughter) {
 727       AliError("The K0/K0bar daughter particle not found in MC array");
 728       delete [] fLabelArray;
 729       fLabelArray = 0x0;
 730       return checkCD;
 731     }
 732
 733     if (!(TMath::Abs(mcPartDaughter->GetPdgCode())==310)) {
 734       AliDebug(2, "The K0/K0bar MC doesn't go in K0S, skipping!!");
 735       delete [] fLabelArray;
 736       fLabelArray = 0x0;
 737       return checkCD;
 738     }
 739
 740     if (mcPartDaughter->GetNDaughters()!=2) {
 741       AliDebug(2, "The K0S MC doesn't decay in 2 particles, skipping!!");
 742       delete [] fLabelArray;
 743       fLabelArray = 0x0;
 744       return checkCD;
 745     }
 746
 747     Int_t daughterD0 = mcPartDaughter->GetDaughter(0);
 748     Int_t daughterD1 = mcPartDaughter->GetDaughter(1);
 749     if (daughterD0<0 || daughterD1<0) {
 750       AliDebug(2, Form("The K0S MC particle doesn't have correct daughters, skipping!!"));
 751       delete [] fLabelArray;
 752       fLabelArray = 0x0;
 753       return checkCD;
 754     }
 755
 756     AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD0));
 757     AliAODMCParticle* mcPartDaughterD1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD1));
 758     if (!mcPartDaughterD0 || !mcPartDaughterD1) {
 759       AliError("The K0S daughter particles not found in MC array");
 760       delete [] fLabelArray;
 761       fLabelArray = 0x0;
 762       return checkCD;
 763     }
 764
 765     if (! ( TMath::Abs(mcPartDaughterD0->GetPdgCode())==211 &&
 766             TMath::Abs(mcPartDaughterD1->GetPdgCode())==211 ) ) {
 767       AliDebug(2, "The K0S MC doesn't decay in pi+ pi-, skipping!!");
 768       delete [] fLabelArray;
 769       fLabelArray = 0x0;
 770       return checkCD;
 771     }
 772
 773     // K0S -> pi+ pi-
 774
 775     fLabelArray[0] = daughter0;//mcPartDaughter0->GetLabel(); // bachelor
 776
 777     if (mcPartDaughterD0->GetPdgCode()==211) {
 778       fLabelArray[1] = daughterD0;//mcPartDaughterD0->GetLabel(); // pi+
 779       fLabelArray[2] = daughterD1;//mcPartDaughterD1->GetLabel(); // pi-
 780       AliDebug(2,Form(" daughter0=%d ------ daughter1=%d ------ dg0->GetLabel()=%d ------ dg1->GetLabel()=%d ",daughterD0,daughterD1,mcPartDaughterD0->GetLabel(),mcPartDaughterD1->GetLabel()));
 781     } else if (mcPartDaughterD1->GetPdgCode()==211) {
 782       fLabelArray[1] = daughterD1;//mcPartDaughterD1->GetLabel(); // pi+
 783       fLabelArray[2] = daughterD0;//mcPartDaughterD0->GetLabel(); // pi-
 784       AliDebug(2,Form(" daughter0=%d ------ daughter1=%d ------ dg0->GetLabel()=%d ------ dg1->GetLabel()=%d ",daughterD1,daughterD0,mcPartDaughterD1->GetLabel(),mcPartDaughterD0->GetLabel()));
 785     }
 786   }
 787
 788   AliDebug(2,Form(" label0=%d, label1=%d, label2=%d",fLabelArray[0],fLabelArray[1],fLabelArray[2]));
 789
 790   SetAccCut(); // setting the pt and eta acceptance cuts
 791
 792   checkCD = kTRUE;
 793   return checkCD;
 794
 795 }
 796 //____________________________________________
 797 Bool_t AliCFVertexingHFLctoV0bachelor::FillVectorFromMCarray(AliAODMCParticle *mcPartDaughterBachelor,
 798                                                              AliAODMCParticle *mcPartDaughterK0,
 799                                                              Double_t *vectorMC)
 800 {
 801   // fill the vector
 802
 803   Bool_t bGenValues = kFALSE;
 804
 805   AliAODMCParticle *mcPartV0DaughterPos = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[1]));
 806   AliAODMCParticle *mcPartV0DaughterNeg = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[2]));
 807   AliAODMCParticle *mcPartDaughterV0 = 0x0;
 808
 809   if(!mcPartV0DaughterPos && !mcPartV0DaughterNeg) return bGenValues;
 810
 811   if (TMath::Abs(mcPartDaughterK0->GetPdgCode())==311) {
 812     Int_t daughterK0 = mcPartDaughterK0->GetDaughter(0);
 813     if (daughterK0<0) {
 814       AliDebug(2, Form("The K0/K0bar particle doesn't have correct daughter, skipping!!"));
 815       return bGenValues;
 816     }
 817     mcPartDaughterV0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterK0));
 818     if (!mcPartDaughterV0) {
 819       AliDebug(2,"The K0/K0bar daughter particle not found in MC array");
 820       return bGenValues;
 821     }
 822     if (TMath::Abs(mcPartDaughterV0->GetPdgCode())!=310) {
 823       AliDebug(2,"The K0/K0bar daughter particle is not a K0S");
 824       return bGenValues;
 825     }
 826   } else if (TMath::Abs(mcPartDaughterK0->GetPdgCode())==3122) {
 827     mcPartDaughterV0 = dynamic_cast<AliAODMCParticle*>(mcPartDaughterK0);
 828     if (!mcPartDaughterV0) {
 829       AliDebug(2,"The Lambda particle not found in MC array");
 830       return bGenValues;
 831     }
 832   }
 833
 834   if (!mcPartDaughterV0) {
 835     AliDebug(2,"V0 particle not found in MC array");
 836     return bGenValues;
 837   }
 838
 839   Double_t cTLc = Ctau(fmcPartCandidate); // by default wrt Primary Vtx
 840   Double_t pTbach = mcPartDaughterBachelor->Pt(); // get the bachelor pT
 841
 842   Double_t vtx1[3] = {0,0,0};   // primary vertex
 843   Bool_t hasPrimVtx = fmcPartCandidate->XvYvZv(vtx1);  // cm
 844
 845   // getting vertex from daughters
 846   Double_t vtx1daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 847   Bool_t hasSecVtx1 = mcPartDaughterBachelor->XvYvZv(vtx1daughter0);  //cm
 848   Double_t vtx1daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 849   Bool_t hasSecVtx2 = mcPartDaughterV0->XvYvZv(vtx1daughter1);  //cm
 850   if (!hasPrimVtx || !hasSecVtx1 || !hasSecVtx2) {
 851     AliDebug(2,"At least one of Prim.vtx, V0vtx, BachelorVtx doesn't exist!");
 852     //return bGenValues;
 853   }
 854
 855   if (TMath::Abs(vtx1daughter0[0]-vtx1daughter1[0])>1E-5 ||
 856       TMath::Abs(vtx1daughter0[1]-vtx1daughter1[1])>1E-5 ||
 857       TMath::Abs(vtx1daughter0[2]-vtx1daughter1[2])>1E-5) {
 858     AliError("Daughters have different secondary vertex, skipping the track");
 859     //return bGenValues;
 860   }
 861
 862   // getting the momentum from the daughters
 863   Double_t px1[2] = {mcPartDaughterBachelor->Px(), mcPartDaughterV0->Px()};
 864   Double_t py1[2] = {mcPartDaughterBachelor->Py(), mcPartDaughterV0->Py()};
 865   Double_t pz1[2] = {mcPartDaughterBachelor->Pz(), mcPartDaughterV0->Pz()};
 866
 867   Int_t nprongs = 2;
 868   Short_t charge = mcPartDaughterBachelor->Charge();
 869   Double_t d0[2] = {0.,0.};
 870   AliAODRecoDecayHF* decayLc = new AliAODRecoDecayHF(vtx1,vtx1daughter0,nprongs,charge,px1,py1,pz1,d0);
 871   Double_t cosPAwrtPrimVtxLc = decayLc->CosPointingAngle();
 872   delete decayLc;
 873
 874   // getting vertex from daughters
 875   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 876   Bool_t hasSecVtx3 = mcPartV0DaughterPos->XvYvZv(vtx2daughter0);  //cm
 877   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 878   Bool_t hasSecVtx4 = mcPartV0DaughterNeg->XvYvZv(vtx2daughter1);  //cm
 879   if (!hasSecVtx3 || !hasSecVtx4) {
 880     AliDebug(2,"At least one of V0Posvtx, V0Negtx doesn't exist!");
 881     //return bGenValues;
 882   }
 883
 884   if (TMath::Abs(vtx2daughter0[0]-vtx2daughter1[0])>1E-5 ||
 885       TMath::Abs(vtx2daughter0[1]-vtx2daughter1[1])>1E-5 ||
 886       TMath::Abs(vtx2daughter0[2]-vtx2daughter1[2])>1E-5) {
 887     AliError("Daughters have different secondary vertex, skipping the track");
 888     //return bGenValues;
 889   }
 890
 891   // getting the momentum from the daughters
 892   Double_t px[2] = {mcPartV0DaughterPos->Px(), mcPartV0DaughterNeg->Px()};
 893   Double_t py[2] = {mcPartV0DaughterPos->Py(), mcPartV0DaughterNeg->Py()};
 894   Double_t pz[2] = {mcPartV0DaughterPos->Pz(), mcPartV0DaughterNeg->Pz()};
 895
 896   nprongs = 2;
 897   charge = 0;
 898   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
 899   Double_t cosPAwrtPrimVtxV0 = decay->CosPointingAngle();
 900   Double_t cTV0 = 0.; //ct
 901   if (fGenLcOption==kCountK0Sp) {
 902     cTV0 = decay->Ct(310); // by default wrt Primary Vtx
 903   } else if (fGenLcOption==kCountLambdapi) {
 904     cTV0 = decay->Ct(3122); // by default wrt Primary Vtx
 905   }
 906
 907   Double_t invMass = 0.; //invMass
 908   if (fGenLcOption==kCountK0Sp) {
 909     invMass = decay->InvMass2Prongs(0,1,211,211);
 910   } else if (fGenLcOption==kCountLambdapi) {
 911     if (fmcPartCandidate->GetPdgCode() == 4122)
 912       invMass = decay->InvMass2Prongs(0,1,2212,211);
 913     else if (fmcPartCandidate->GetPdgCode() ==-4122)
 914       invMass = decay->InvMass2Prongs(0,1,211,2212);
 915   }
 916   delete decay;
 917
 918   vectorMC[0]  = fmcPartCandidate->Pt();
 919   vectorMC[1]  = fmcPartCandidate->Y() ;
 920   vectorMC[2]  = fmcPartCandidate->Phi();
 921   vectorMC[3]  = cosPAwrtPrimVtxV0;
 922   vectorMC[4]  = 0; // dummy value x MC, onTheFlyStatus
 923   vectorMC[5]  = fCentValue; // reconstructed centrality
 924   vectorMC[6]  = 1; // dummy value x MC, fFake
 925   vectorMC[7]  = fMultiplicity; // reconstructed multiplicity
 926
 927   if (fConfiguration==AliCFTaskVertexingHF::kSnail) {
 928     vectorMC[8]  = pTbach;
 929     vectorMC[9]  = mcPartV0DaughterPos->Pt();
 930     vectorMC[10] = mcPartV0DaughterNeg->Pt();
 931     vectorMC[11] = invMass;
 932     vectorMC[12] = 0; // dummy value x MC, V0 DCA
 933     vectorMC[13] = cTV0*1.E4; // in micron
 934     vectorMC[14] = cTLc*1.E4; // in micron
 935     vectorMC[15] = cosPAwrtPrimVtxLc;
 936   }
 937
 938   bGenValues = kTRUE;
 939   return bGenValues;
 940
 941 }