Updating classes to be able to use the common AliCFVertexingHFCascade.h/cxx class...
authorzconesa <zaida.conesa.del.valle@cern.ch>
Mon, 6 Oct 2014 13:33:23 +0000 (15:33 +0200)
committerzconesa <zaida.conesa.del.valle@cern.ch>
Mon, 6 Oct 2014 13:33:59 +0000 (15:33 +0200)
PWGHF/vertexingHF/AliAODRecoCascadeHF.cxx
PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx [changed mode: 0755->0644]
PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.h [changed mode: 0755->0644]
PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx
PWGHF/vertexingHF/AliCFTaskVertexingHF.h [changed mode: 0755->0644]
PWGHF/vertexingHF/AliCFVertexingHF.cxx
PWGHF/vertexingHF/AliCFVertexingHFCascade.cxx
PWGHF/vertexingHF/AliCFVertexingHFCascade.h

index ddfe39c..c89209c 100644 (file)
@@ -130,7 +130,7 @@ Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
   if ( isV0 &&
        ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
         (pdgDg[1]==211 && pdgDg[0]==3122) ) ) {
-    AliWarning("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping!");
+    AliWarning(Form("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping! (pdgDg[0] = %d, (pdgDg[1] = %d)", pdgDg[0], pdgDg[1]));
     return -1;
   }
 
old mode 100755 (executable)
new mode 100644 (file)
index 0321495..fbd1506
@@ -187,7 +187,8 @@ AliAnalysisTaskSE(),
   fCutKFDeviationFromVtxV0(0.),
   fCurrentEvent(-1),
   fBField(0),
-  fKeepingOnlyPYTHIABkg(kFALSE)
+  fKeepingOnlyPYTHIABkg(kFALSE),
+  fHistoMCLcK0Sp(0x0)
 {
   //
   // Default ctor
@@ -307,7 +308,8 @@ AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Cha
   fCutKFDeviationFromVtxV0(0.),
   fCurrentEvent(-1),
   fBField(0),
-  fKeepingOnlyPYTHIABkg(kFALSE)
+  fKeepingOnlyPYTHIABkg(kFALSE),
+  fHistoMCLcK0Sp(0x0)
 
 {
   //
@@ -415,6 +417,10 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
     AliError("fOutput not available");
     return;
   }
+
+  
+  AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0Sp->GetEntries()));
+
   fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
   if (!fOutputKF) {     
     AliError("fOutputKF not available");
@@ -605,6 +611,10 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
   //fOutput->Add(fVariablesTreeSgn);
   //fOutput->Add(fVariablesTreeBkg);
 
+  const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
+
+  fHistoMCLcK0Sp = new TH1F("fHistoMCLcK0Sp", "fHistoMCLcK0Sp", 14, ptbins);
+
   fOutput->Add(fHistoEvents);
   fOutput->Add(fHistoLc);
   fOutput->Add(fHistoLcOnTheFly);
@@ -614,6 +624,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
   fOutput->Add(fHistoCodesBkg);
   fOutput->Add(fHistoLcpKpiBeforeCuts);
   fOutput->Add(fHistoBackground);
+  fOutput->Add(fHistoMCLcK0Sp);
 
   PostData(1, fOutput);
   PostData(4, fVariablesTreeSgn);
@@ -824,7 +835,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
   }
 
   fCurrentEvent++;
-  Printf("Processing event = %d", fCurrentEvent);
+  AliDebug(2, Form("Processing event = %d", fCurrentEvent));
   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
   TClonesArray *arrayLctopKos=0;
 
@@ -859,23 +870,6 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
 
   fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
 
-  // AOD primary vertex
-  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
-  if (!fVtx1) return;
-  if (fVtx1->GetNContributors()<1) return;
-
-  fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
-
-  if ( !fIsEventSelected ) {
-    fHistoEvents->Fill(0);
-    return; // don't take into account not selected events 
-  }
-  fHistoEvents->Fill(1);
-
-  // Setting magnetic field for KF vertexing
-  fBField = aodEvent->GetMagneticField();
-  AliKFParticle::SetField(fBField);
-
   // mc analysis 
   TClonesArray *mcArray = 0;
   AliAODMCHeader *mcHeader=0;
@@ -893,8 +887,27 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
       AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
       return;
     }
+    //Printf("Filling MC histo");
+    FillMCHisto(mcArray);
   }
 
+  // AOD primary vertex
+  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+  if (!fVtx1) return;
+  if (fVtx1->GetNContributors()<1) return;
+
+  fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
+
+  if ( !fIsEventSelected ) {
+    fHistoEvents->Fill(0);
+    return; // don't take into account not selected events 
+  }
+  fHistoEvents->Fill(1);
+
+  // Setting magnetic field for KF vertexing
+  fBField = aodEvent->GetMagneticField();
+  AliKFParticle::SetField(fBField);
+
   Int_t nSelectedAnal = 0;
   if (fIsK0sAnalysis) {
     MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
@@ -910,6 +923,103 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
   PostData(6, fOutputKF);
 
 }
+//-------------------------------------------------------------------------------
+void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){
+
+  // method to fill MC histo: how many Lc --> K0S + p are there at MC level
+  for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
+    AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+    if (!mcPart){
+      AliError("Failed casting particle from MC array!, Skipping particle");
+      AliInfo("Failed casting particle from MC array!, Skipping particle");
+      continue;
+    }
+    Int_t pdg = mcPart->GetPdgCode();
+    if (TMath::Abs(pdg) != 4122){
+      AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
+      continue;
+    }
+    AliInfo(Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
+    Int_t labeldaugh0 = mcPart->GetDaughter(0);
+    Int_t labeldaugh1 = mcPart->GetDaughter(1);
+    if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
+      AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
+      continue;
+    }
+    else if (labeldaugh1 - labeldaugh0 == 1){
+      AliInfo(Form("Step 1 ok: The MC particle has correct daughters!!"));
+      AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
+      AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
+      Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
+      Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
+      AliAODMCParticle* bachelorMC = daugh0;
+      AliAODMCParticle* v0MC = daugh1;
+      AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
+      if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){ 
+       // we are in the case of Lc --> K0 + p; now we have to check if the K0 decays in K0S, and if this goes in pi+pi-
+       /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it:
+       if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
+         bachelorMC = daugh1;
+         v0MC = daugh0;
+       }
+       AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
+       if (v0MC->GetNDaughters() != 1) { 
+         AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
+         continue;
+       }
+       else { // So far: Lc --> K0 + p, K0 with 1 daughter 
+         AliInfo("Step 2 ok: The K0 does decay in 1 body only! ");
+         Int_t labelK0daugh = v0MC->GetDaughter(0);
+         AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
+         if(!partK0S){
+           AliError("Error while casting particle! returning a NULL array");
+           AliInfo("Error while casting particle! returning a NULL array");
+           continue;
+         }
+         else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
+           if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
+             AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
+             continue;
+           }
+           else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
+             AliInfo("Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
+             Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0);
+             Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1);
+             AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
+             AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
+             if (!daughK0S0 || ! daughK0S1){
+               AliDebug(2, "Could not access K0S daughters, continuing...");
+               continue;
+             }
+             else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
+               AliInfo("Step 4 ok: Could access K0S daughters, continuing...");
+               Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
+               Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
+               if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
+                 AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
+                 AliInfo("The K0S does not decay in pi+pi-, continuing");
+               }
+               else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
+                 if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
+                   AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
+                   fHistoMCLcK0Sp->Fill(mcPart->Pt());
+                 }
+                 else {
+                   AliDebug(2, "not in fiducial acceptance! Skipping");
+                   continue;
+                 }
+               }
+             }
+           }
+         }
+       }
+      }
+    }
+  } // closing loop over mcArray
+
+  return; 
+
+}
 
 //-------------------------------------------------------------------------------
 void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
@@ -976,7 +1086,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *a
       // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
       fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
       if (fmcLabelLc>=0) {
-       AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
+       AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
 
        AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
        if(partLc){
old mode 100755 (executable)
new mode 100644 (file)
index d66ddad..060730d
@@ -133,6 +133,8 @@ class AliAnalysisTaskSELc2V0bachelorTMVA : public AliAnalysisTaskSE
                        Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
                        Double_t* distances, Double_t* armPolKF);
 
+  void FillMCHisto(TClonesArray *mcArray);
+
   AliAnalysisTaskSELc2V0bachelorTMVA(const AliAnalysisTaskSELc2V0bachelorTMVA &source);
   AliAnalysisTaskSELc2V0bachelorTMVA& operator=(const AliAnalysisTaskSELc2V0bachelorTMVA& source); 
   
@@ -253,8 +255,9 @@ class AliAnalysisTaskSELc2V0bachelorTMVA : public AliAnalysisTaskSE
   Int_t fCurrentEvent;              // current event number - for debug purposes
   Double_t fBField;                   // magnetic field of current event
   Bool_t fKeepingOnlyPYTHIABkg;       // flag to allow to use only PYTHIA tracks for background
+  TH1F* fHistoMCLcK0Sp;               // histo with MC Lc --> K0S + p 
 
-  ClassDef(AliAnalysisTaskSELc2V0bachelorTMVA, 4); // class for Lc->p K0
+  ClassDef(AliAnalysisTaskSELc2V0bachelorTMVA, 5); // class for Lc->p K0
 };
 
 #endif
index 6eeffa4..7944c8b 100644 (file)
@@ -131,7 +131,9 @@ AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
   fRefMult(9.26),
   fZvtxCorrectedNtrkEstimator(kFALSE),
   fIsPPData(kFALSE),
-  fIsPPbData(kFALSE)
+  fIsPPbData(kFALSE),
+  fUseAdditionalCuts(kFALSE),
+  fUseCutsForTMVA(kFALSE)
 {
   //
   //Default ctor
@@ -190,7 +192,9 @@ AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts
   fRefMult(9.26),
   fZvtxCorrectedNtrkEstimator(kFALSE),
   fIsPPData(kFALSE),
-  fIsPPbData(kFALSE)
+  fIsPPbData(kFALSE),
+  fUseAdditionalCuts(kFALSE),
+  fUseCutsForTMVA(kFALSE)
 {
   //
   // Constructor. Initialization of Inputs and Outputs
@@ -565,8 +569,11 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
   }
 
   AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
-  if (!aodVtx) return;
-       
+  if (!aodVtx) {
+    AliDebug(3, "The event was skipped due to missing vertex");
+    return;
+  }
+
   if (!arrayBranch) {
     AliError("Could not find array of HF vertices");
     return;
@@ -617,10 +624,27 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
   }
   case 21:{ 
     cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(413);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(211);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(421);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(421);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(321);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
     break;
   }
   case 22:{
-    cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption); // Lc -> K0S+proton
+    // Lc ->  K0S+proton
+    //    cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption); 
+    cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(4122);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(2212);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(310);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(311);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(211);
+    ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
+    if (fUseAdditionalCuts) ((AliCFVertexingHFCascade*)cfVtxHF)->SetUseCutsForTMVA(fUseCutsForTMVA);
     break;
   }
   case 31:
@@ -735,6 +759,7 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
   if (fUseMCVertex) fCuts->SetUseMCVertex(); 
 
   if (fCentralitySelection){ // keep only the requested centrality
+
     if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
       delete[] containerInput;
       delete[] containerInputMC;
@@ -759,7 +784,6 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
   Double_t multiplicity = nTracklets; // set to the Ntracklet estimator
   if(fMultiplicityEstimator==kVZERO) { multiplicity = vzeroMult; }
 
-
   cfVtxHF->SetMultiplicity(multiplicity);
 
   //  printf("Multiplicity estimator %d, value %2.2f\n",fMultiplicityEstimator,multiplicity);
@@ -779,26 +803,29 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
       AliDebug(2,"Check the MC-level cuts - not desidered particle");
       continue;  // 0 stands for MC level
     }
+    else {
+      AliDebug(3, Form("\n\n---> COOL! we found a particle (particle %d)!!! with PDG code = %d \n\n", iPart, mcPart->GetPdgCode()));
+    }
     cfVtxHF->SetMCCandidateParam(iPart);
         
          
     if (!(cfVtxHF->SetLabelArray())){
-      AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
+      AliDebug(2,Form("Impossible to set the label array for particle %d (decaychannel = %d)", iPart, fDecayChannel));
       continue;
     }             
 
     //check the candiate family at MC level
     if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
-      AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
+      AliDebug(2,Form("Check on the family wrong for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
       continue;
     }
     else{
-      AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
+      AliDebug(2,Form("Check on the family OK for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
     }
                
     //Fill the MC container
     Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
-    AliDebug(2,Form("mcContainerFilled = %d)",mcContainerFilled));
+    AliDebug(2, Form("particle = %d mcContainerFilled = %d",iPart, mcContainerFilled));
     if (mcContainerFilled) {
       if (fUseWeight){
        if (fFuncWeight){ // using user-defined function
@@ -811,12 +838,19 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
        }
        AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
       }
-      if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
+      if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) {
+       AliDebug(3, Form("Not in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0], containerInputMC[1]));
+       continue;
+      }
+      else{
+       AliDebug(3, Form("YES!! in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0],containerInputMC[1]));
+      }
+
       //MC Limited Acceptance
       if (TMath::Abs(containerInputMC[1]) < 0.5) {
        fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
        AliDebug(3,"MC Lim Acc container filled\n");
-      }            
+      }
                        
       //MC 
       fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
@@ -928,6 +962,7 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
     AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
 
     Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
+    AliDebug(3, Form("CF task: RecoContFilled for candidate %d is %d", iCandid, (Int_t)recoContFilled));
     if (recoContFilled){
 
       // weight according to pt
@@ -949,6 +984,7 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
       }                
       //Reco Step
       Bool_t recoStep = cfVtxHF->RecoStep();
+      if (recoStep) AliDebug(2, Form("CF task: Reco step for candidate %d is %d", iCandid, (Int_t)recoStep));
       Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
                        
 
@@ -1004,9 +1040,13 @@ void AliCFTaskVertexingHF::UserExec(Option_t *)
            else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
              Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
              if (keepLctoV0bachelor) recoAnalysisCuts=3;
-           }
-
-                                                   
+             AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+             AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+             AliPIDResponse* pidResponse = inputHandler->GetPIDResponse();
+             if (fUseAdditionalCuts){
+               if (!((AliCFVertexingHFCascade*)cfVtxHF)->CheckAdditionalCuts(pidResponse)) recoAnalysisCuts = -1;
+             }
+           }                                       
 
            fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object 
            Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
old mode 100755 (executable)
new mode 100644 (file)
index 369298e..4e4d971
@@ -234,6 +234,11 @@ public:
 
        Bool_t ProcessLctoV0Bachelor(Int_t returnCodeDs) const;
 
+       void SetUseAdditionalCuts(Bool_t flag) { fDecayChannel == 22 ? fUseAdditionalCuts = flag : fUseAdditionalCuts = kFALSE;}
+       Bool_t GetUseAdditionalCuts() const { return fUseAdditionalCuts; }
+       void SetUseCutsForTMVA(Bool_t useCutsForTMVA) { fDecayChannel == 22 ? fUseCutsForTMVA = useCutsForTMVA : fUseAdditionalCuts = kFALSE;}
+       Bool_t GetUseCutsForTMVA() const {return fUseCutsForTMVA;}
+
 protected:
        AliCFManager   *fCFManager;   //  pointer to the CF manager
        TH1I *fHistEventsProcessed;   //! simple histo for monitoring the number of events processed
@@ -287,8 +292,11 @@ protected:
        Bool_t fZvtxCorrectedNtrkEstimator; // flag to use the z-vtx corrected (if not use uncorrected) multiplicity estimator
        Bool_t fIsPPData; // flag for pp data (not checking centrality)
        Bool_t fIsPPbData; // flag for pPb data (used for multiplicity corrections)
+       Bool_t fUseAdditionalCuts;  // flag to use additional cuts needed for Lc --> K0S + p, TMVA
+       Bool_t fUseCutsForTMVA;     // flag to use additional cuts needed for Lc --> K0S + p, TMVA
+                                   // these are the pre-selection cuts for the TMVA
    
-       ClassDef(AliCFTaskVertexingHF,20); // class for HF corrections as a function of many variables
+       ClassDef(AliCFTaskVertexingHF,21); // class for HF corrections as a function of many variables
 };
 
 #endif
index 249c0bb..cde20ed 100644 (file)
@@ -302,18 +302,20 @@ Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClones
 
        Int_t pdgGranma = CheckOrigin();
 
+       AliDebug(3, Form("pdgGranma = %d", pdgGranma));
+
        if (pdgGranma == -99999){
-               AliDebug(2,"This particle does not have a quark in his genealogy\n");
-               return kFALSE;
+         AliDebug(2, "This particle does not have a quark in his genealogy\n");
+         return kFALSE;
        }
        if (pdgGranma == -9999){
-               AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");   
-               return kFALSE;
+         AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n"); 
+         return kFALSE;
        }       
        
        if (pdgGranma == -999){
-               AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");  
-               return kFALSE;
+         AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");  
+         return kFALSE;
        }       
        
        if (!CheckMCDaughters()) {
@@ -321,8 +323,8 @@ Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClones
          return kFALSE;
        }
        if (!CheckMCChannelDecay()) {
-               AliDebug(3,"CheckMCChannelDecay false");
-               return kFALSE;
+         AliDebug(3, "CheckMCChannelDecay false");
+         return kFALSE;
        }
        return kTRUE;
 }
@@ -354,18 +356,25 @@ Int_t AliCFVertexingHF::CheckOrigin() const
                        }
                        if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
                        mother = mcGranma->GetMother();
+                       AliDebug(3, Form("mother = %d", mother));
                }else{
                        AliError("Failed casting the mother particle!");
                        break;
                }
        }
 
-       if(fRejectIfNoQuark && !isQuarkFound) return -99999;
+       if(fRejectIfNoQuark && !isQuarkFound) {
+         return -99999;
+       }
        if(isFromB){
-         if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
+         if (!fKeepDfromB) {
+           return -9999; //skip particle if come from a B meson.
+         }
        }
        else{
-         if (fKeepDfromBOnly) return -999;
+         if (fKeepDfromBOnly) {
+           return -999;
+         }
        }
        return pdgGranma;
 }
@@ -382,22 +391,27 @@ Bool_t AliCFVertexingHF::CheckMCDaughters()const
        
        Int_t label0 = fmcPartCandidate->GetDaughter(0);
        Int_t label1 = fmcPartCandidate->GetDaughter(1);
-       AliDebug(3,Form("label0 = %d, label1 = %d",label0,label1));
+       AliDebug(3,Form("label0 = %d, label1 = %d", label0, label1));
        if (label1<=0 || label0 <= 0){
-               AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
-               return checkDaughters;  
+         AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
+         return checkDaughters;  
        }
        
        if (fLabelArray == 0x0) {
-               return checkDaughters;
+         return checkDaughters;
        }  
-
+       
        for (Int_t iProng = 0; iProng<fProngs; iProng++){
-               mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));    
-               if (!mcPartDaughter) {
-                       AliWarning("At least one Daughter Particle not found in tree, skipping"); 
-                       return checkDaughters;  
-               }
+         AliDebug(3, Form("fLabelArray[%d] = %d", iProng, fLabelArray[iProng]));
+       }
+       AliDebug(3, Form("Entries in MC array = %d (fast  = %d)", fmcArray->GetEntries(), fmcArray->GetEntriesFast()));
+       for (Int_t iProng = 0; iProng<fProngs; iProng++){
+         AliDebug(3, Form("fLabelArray[%d] = %d", iProng, fLabelArray[iProng]));
+         mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));    
+         if (!mcPartDaughter) {
+           AliWarning("At least one Daughter Particle not found in tree, skipping"); 
+           return checkDaughters;  
+         }
        }
        
        checkDaughters = kTRUE;
@@ -422,6 +436,9 @@ Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
                }               
                mcContainerFilled = kTRUE;              
        }
+       else {
+         AliDebug(3, "We could not fill the array for the container");
+       }
        delete [] vectorMC;
        vectorMC = 0x0;
        return mcContainerFilled;       
@@ -621,11 +638,11 @@ Bool_t AliCFVertexingHF::RecoStep()
        Int_t pdgGranma = CheckOrigin();
        
        if (pdgGranma == -99999){
-               AliDebug(2,"This particle does not have a quark in his genealogy\n");
+                       AliDebug(2,"This particle does not have a quark in his genealogy\n");
                return bRecoStep;
        }
        if (pdgGranma == -9999){
-               AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only prompt charm particles\n");
+               AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only prompt charm particles\n");               
                return bRecoStep;
        }
 
@@ -771,9 +788,13 @@ Bool_t AliCFVertexingHF::SetLabelArray()
        AliAODMCParticle *mcPartDaughter;
        Int_t label0 = fmcPartCandidate->GetDaughter(0);
        Int_t label1 = fmcPartCandidate->GetDaughter(1);
-       AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
+       //AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
+       AliAODMCParticle* tmp0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0));
+       AliAODMCParticle* tmp1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label1));
+
+       AliDebug(2, Form("label0 = %d (pdg = %d), label1 = %d (pdg = %d)", label0, tmp0->GetPdgCode(), label1, tmp1->GetPdgCode()));
        if (label1<=0 || label0 <= 0){
-               AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
+               AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
                delete [] fLabelArray; 
                fLabelArray = 0x0;  
                return bLabelArray;
@@ -796,6 +817,7 @@ Bool_t AliCFVertexingHF::SetLabelArray()
        }
        // resonant decay channel
        else if (label1 - label0 == fProngs-2 && fProngs > 2){
+         AliDebug(3, "In the resonance decay channel");
                Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
                Int_t foundDaughters = 0;
                for(Int_t iDau=0; iDau<fProngs-1; iDau++){
@@ -808,9 +830,11 @@ Bool_t AliCFVertexingHF::SetLabelArray()
                           return bLabelArray;
                         }  
                        Int_t pdgCode=TMath::Abs(part->GetPdgCode());
+                       AliDebug(3, Form("Prong %d had pdg = %d", iDau, pdgCode));
                        if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
                                if (part) {
                                        fLabelArray[foundDaughters] = part->GetLabel();
+                                       AliDebug(3, Form("part found at %d has label = %d", iLabelDau, part->GetLabel()));
                                        foundDaughters++;
                                }
                                else{
@@ -822,6 +846,7 @@ Bool_t AliCFVertexingHF::SetLabelArray()
                        }
                        // added K0S case - Start
                        else if (pdgCode==311) {
+                         AliDebug(3, "K0S case");
                          if (part->GetNDaughters()!=1) {
                            delete [] fLabelArray; 
                            fLabelArray = 0x0;  
@@ -835,19 +860,22 @@ Bool_t AliCFVertexingHF::SetLabelArray()
                            fLabelArray = 0x0;
                            return bLabelArray;
                          }                 
-                         Int_t nDauRes=partK0S->GetNDaughters();
-                         if(nDauRes!=2 || partK0S->GetPdgCode()!=310) {
-                           AliDebug(2,"No K0S on no 2-body decay");
+                         Int_t nDauRes = partK0S->GetNDaughters();
+                         AliDebug(3, Form("nDauRes = %d", nDauRes));
+                         if (nDauRes!=2 || partK0S->GetPdgCode() != 310) {
+                           AliDebug(2, "No K0S on no 2-body decay");
                            delete [] fLabelArray;
                            fLabelArray = 0x0;
                            return bLabelArray;
                          }
                          Int_t labelFirstDauRes = partK0S->GetDaughter(0);
-                         AliDebug(2,Form(" Found K0S (%d)",labelK0Dau));
+                         AliDebug(2, Form(" Found K0S (%d)", labelK0Dau));
                          for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
                            Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
                            AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
+                           AliDebug(3, Form("daughter = %d, pointer = %p", iLabelDauRes, dauRes));
                            if (dauRes){
+                             AliDebug(3, Form("PDG code = %d", dauRes->GetPdgCode()));
                              if (TMath::Abs(dauRes->GetPdgCode())!=211) {
                                AliDebug(2,"K0S doesn't decay in 2 charged pions!");
                                delete [] fLabelArray; 
@@ -856,6 +884,7 @@ Bool_t AliCFVertexingHF::SetLabelArray()
                              }
                              else {
                                fLabelArray[foundDaughters] = dauRes->GetLabel();
+
                                foundDaughters++;
                              }
                            }
@@ -869,44 +898,56 @@ Bool_t AliCFVertexingHF::SetLabelArray()
                        }
                        // added K0S case - End
                        else{
-                               Int_t nDauRes=part->GetNDaughters();
-                               if(nDauRes!=2) {
-                                       delete [] fLabelArray; 
-                                       fLabelArray = 0x0;  
-                                       return bLabelArray;
-                               }
-                               Int_t labelFirstDauRes = part->GetDaughter(0); 
-                               for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
-                                       Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
-                                       AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
-                                       if (dauRes){
-                                               fLabelArray[foundDaughters] = dauRes->GetLabel();
-                                               foundDaughters++;
-                                       }
-                                       else{
-                                               AliError("Error while casting resonant daughter! returning a NULL array");
-                                               delete [] fLabelArray; 
-                                               fLabelArray = 0x0;  
-                                               return bLabelArray;
-                                       }
-                               }
+                         Int_t nDauRes=part->GetNDaughters();
+                         AliDebug(3, Form("nDauRes = %d", nDauRes));
+                         if(nDauRes!=2) {
+                           AliDebug(3, Form("nDauRes = %d, different from 2", nDauRes));
+                           Int_t labelFirstDauResTest = part->GetDaughter(0);
+                           for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
+                             Int_t iLabelDauResTest = labelFirstDauResTest+iDauRes;
+                             AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauResTest));
+                             if (dauRes){
+                               AliDebug(3, Form("pdg of daugh %d = %d", iDauRes, dauRes->GetPdgCode()));
+                             }
+                           }
+                           delete [] fLabelArray; 
+                           fLabelArray = 0x0;  
+                           return bLabelArray;                     
+                         }
+                         Int_t labelFirstDauRes = part->GetDaughter(0); 
+                         for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
+                           Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
+                           AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
+                           if (dauRes){
+                             fLabelArray[foundDaughters] = dauRes->GetLabel();
+                             foundDaughters++;
+                           }
+                           else{
+                             AliError("Error while casting resonant daughter! returning a NULL array");
+                             delete [] fLabelArray; 
+                             fLabelArray = 0x0;  
+                             return bLabelArray;
+                           }
+                         }
                        }
                }
                if (foundDaughters != fProngs){
-                       delete [] fLabelArray; 
-                       fLabelArray = 0x0;  
-                       return bLabelArray;
+                 AliDebug(3, Form("foundDaughters (%d) != fProngs (%d)", foundDaughters, fProngs));
+                 delete [] fLabelArray; 
+                 fLabelArray = 0x0;  
+                 return bLabelArray;
                }
        }
        // wrong correspondance label <--> prongs
        else{
-               delete [] fLabelArray; 
-               fLabelArray = 0x0;  
-               return bLabelArray;
-       }
-       SetAccCut();   // setting the pt and eta acceptance cuts
-       bLabelArray = kTRUE;
-       return bLabelArray;
+         delete [] fLabelArray; 
+         fLabelArray = 0x0;  
+         return bLabelArray;
+       }
+        AliDebug(3, "Setting AccCuts");
+        SetAccCut();   // setting the pt and eta acceptance cuts
+        bLabelArray = kTRUE;
+        return bLabelArray;
 }
 
 //___________________________________________________________
index ff1aa13..a216bc2 100644 (file)
@@ -16,7 +16,7 @@
 //-----------------------------------------------------------------------
 // Class for HF corrections as a function of many variables and steps
 // For D* and other cascades
-// 
+//
 // Author : A.Grelli a.grelli@uu.nl  UTECHT
 //-----------------------------------------------------------------------
 
 #include "AliCFVertexingHFCascade.h"
 #include "AliCFContainer.h"
 #include "AliCFTaskVertexingHF.h"
+#include "AliPIDResponse.h"
+#include "AliPID.h"
 
 ClassImp(AliCFVertexingHFCascade)
 
 
 //_________________________________________
-  AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
-  AliCFVertexingHF(mcArray, originDselection)
-{      
+AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
+AliCFVertexingHF(mcArray, originDselection),
+  fPDGcascade(0),
+  fPDGbachelor(0),
+  fPDGneutrDaugh(0),
+  fPDGneutrDaughPositive(0),
+  fPDGneutrDaughNegative(0),
+  fPrimVtx(0x0)
+{
   // standard constructor
 
   SetNProngs(3);
-  fPtAccCut=new Float_t[fProngs];
-  fEtaAccCut=new Float_t[fProngs];
+  fPtAccCut = new Float_t[fProngs];
+  fEtaAccCut = new Float_t[fProngs];
   // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
-  fPtAccCut[0]=0.;
-  fEtaAccCut[0]=0.;
+  fPtAccCut[0] = 0.;
+  fEtaAccCut[0] = 0.;
   for(Int_t iP=1; iP<fProngs; iP++){
-         fPtAccCut[iP]=0.1;
-         fEtaAccCut[iP]=0.9;
+    fPtAccCut[iP] = 0.1;
+    fEtaAccCut[iP] = 0.9;
   }
 
 }
@@ -63,7 +71,7 @@ AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexing
   if  (this != &c) {
 
     AliCFVertexingHF::operator=(c);
-   
+
   }
   return *this;
 }
@@ -72,40 +80,46 @@ AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexing
 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
 {
   // set the AliAODRecoDecay candidate
-  
+
   Bool_t bSignAssoc = kFALSE;
+
   fRecoCandidate = recoCand;
-  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)recoCand;
+  AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
 
   if (!fRecoCandidate) {
     AliError("fRecoCandidate not found, problem in assignement\n");
     return bSignAssoc;
   }
-  
+
   if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
-  //if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
-  
+
   //Int_t pdgCand = 413;
 
-  Int_t pdgDgDStartoD0pi[2]={421,211};
-  Int_t pdgDgD0toKpi[2]={321,211};
+  Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
+  Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
 
   Int_t nentries = fmcArray->GetEntriesFast();
 
   AliDebug(3,Form("nentries = %d\n", nentries));
  
-  Int_t mcLabel =  dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray); 
+  Bool_t isV0 = kFALSE;
+  if (fPDGcascade == 4122) {
+    isV0 = kTRUE;
+    pdgDgCascade[0] = fPDGbachelor;
+    pdgDgCascade[1] = fPDGneutrDaugh;
+  }
+  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));
+ Int_t mcLabel = cascade->MatchToMC(fPDGcascade, fPDGneutrDaugh, pdgDgCascade, pdgDgNeutrDaugh, fmcArray, isV0); 
   
   if (mcLabel == -1) return bSignAssoc;
 
   if (fRecoCandidate->NumberOfFakeDaughters()>0){
-         fFake = 0;    // fake candidate
-         if (fFakeSelection==1) return bSignAssoc;
+    fFake = 0;    // fake candidate
+    if (fFakeSelection == 1) return bSignAssoc;
   }
   if (fRecoCandidate->NumberOfFakeDaughters()==0){
-         fFake = 2;    // non-fake candidate
-         if (fFakeSelection==2) return bSignAssoc;
+    fFake = 2;    // non-fake candidate
+    if (fFakeSelection == 2) return bSignAssoc;
   }
   
   SetMCLabel(mcLabel);
@@ -126,224 +140,245 @@ Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vecto
   // 
   // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
   //
-       
-       Bool_t bGenValues = kFALSE;
-       
-       
-       Int_t daughter0ds = fmcPartCandidate->GetDaughter(0);
-       Int_t daughter1ds = fmcPartCandidate->GetDaughter(1);
 
-        //the D0
-       AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0ds));
-       AliAODMCParticle* mcPartDaughterPis = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1ds));
+  Bool_t bGenValues = kFALSE;
 
-       if(!mcPartDaughterD0 || !mcPartDaughterPis) return kFALSE;
+  Int_t daughter0cascade = fmcPartCandidate->GetDaughter(0);
+  Int_t daughter1cascade = fmcPartCandidate->GetDaughter(1);
 
-       Double_t vtx1[3] = {0,0,0};   // primary vertex         
-       Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
-       Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
-       fmcPartCandidate->XvYvZv(vtx1);  // cm
+  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0cascade));
+  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1cascade));
+  AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
+  AliAODMCParticle* mcPartDaughterBachelor = NULL;
 
-       //D0 daughters
-       Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
-       Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);
+  // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
+  // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the 
+  // charge of the daughters to decide which is which
+  if (mcPartDaughter0->Charge()/3 == 0){
+    mcPartDaughterNeutrDaugh = mcPartDaughter0;
+    mcPartDaughterBachelor =  mcPartDaughter1;
+  }
+  else {
+    mcPartDaughterNeutrDaugh = mcPartDaughter1;
+    mcPartDaughterBachelor =  mcPartDaughter0;
+  }
 
-       AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
-       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis
+  if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
 
-       if(!mcPartDaughter0 || !mcPartDaughter1) return kFALSE;
+  Double_t vtx1[3] = {0,0,0};   // primary vertex
+  Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
+  Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
+  fmcPartCandidate->XvYvZv(vtx1);  // cm
 
-       // getting vertex from daughters
-       mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
-       mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
-       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) {
-         AliError("Daughters have different secondary vertex, skipping the track");
-         return bGenValues;
-       }
-       
-       Int_t nprongs = 2;
-       Short_t charge = 0;
-       // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
-       AliAODMCParticle* positiveDaugh = mcPartDaughter0;
-       AliAODMCParticle* negativeDaugh = mcPartDaughter1;
-       if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
-               // inverting in case the positive daughter is the second one
-               positiveDaugh = mcPartDaughter1;
-               negativeDaugh = mcPartDaughter0;
-       }
-       // getting the momentum from the daughters
-       Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};            
-       Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};            
-       Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
-       
-       Double_t d0[2] = {0.,0.};               
-       
-       AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
+  //Daughters of the neutral particle of the cascade
+  Int_t daughter0 = mcPartDaughterNeutrDaugh->GetDaughter(0); // this is the positive
+  Int_t daughter1 = mcPartDaughterNeutrDaugh->GetDaughter(1); // this is the negative
+
+  AliAODMCParticle* mcPartNeutrDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
+  AliAODMCParticle* mcPartNeutrDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
+
+  if (!mcPartNeutrDaughter0 || !mcPartNeutrDaughter1) return kFALSE;
+
+  // getting vertex from daughters
+  mcPartNeutrDaughter0->XvYvZv(vtx2daughter0);  // cm
+  mcPartNeutrDaughter1->XvYvZv(vtx2daughter1);  //cm
+  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) {
+    AliError("Daughters have different secondary vertex, skipping the track");
+    return bGenValues;
+  }
+
+  Int_t nprongs = 2;
+  Short_t charge = 0;
+  // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
+  AliAODMCParticle* positiveDaugh = mcPartNeutrDaughter0;
+  AliAODMCParticle* negativeDaugh = mcPartNeutrDaughter1;
+  if (mcPartNeutrDaughter0->GetPdgCode() < 0 && mcPartNeutrDaughter1->GetPdgCode() > 0){
+    // inverting in case the positive daughter is the second one
+    positiveDaugh = mcPartNeutrDaughter1;
+    negativeDaugh = mcPartNeutrDaughter0;
+  }
+
+  // getting the momentum from the daughters
+  Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};         
+  Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};         
+  Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
+
+  Double_t d0[2] = {0.,0.};
+
+  AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1, vtx2daughter0, nprongs, charge, px, py, pz, d0);
        
-       Double_t cosThetaStar = 0.;
-       Double_t cosThetaStarD0 = 0.;
-       Double_t cosThetaStarD0bar = 0.;
-       cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
-       cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
-       if (mcPartDaughterD0->GetPdgCode() == 421){  // D0
-         AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
-         cosThetaStar = cosThetaStarD0;
-       }
-       else if (mcPartDaughterD0->GetPdgCode() == -421){  // D0bar{
-         AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
-         cosThetaStar = cosThetaStarD0bar;
-       }
-       else{
-         AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
-         delete decay;
-         return vectorMC;
-       }
-       if (cosThetaStar < -1 || cosThetaStar > 1) {
-         AliWarning("Invalid value for cosine Theta star, returning");
-         delete decay;
-         return bGenValues;
-       }
+  Double_t cosThetaStar = 0.;
+  Double_t cosThetaStarNeutrDaugh = 0.;
+  Double_t cosThetaStarNeutrDaughBar = 0.;
+  cosThetaStarNeutrDaugh = decay->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
+  cosThetaStarNeutrDaughBar = decay->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughNegative, fPDGneutrDaughPositive);
+  if (mcPartDaughterNeutrDaugh->GetPdgCode() == fPDGneutrDaughForMC){  // neutral particle
+    AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
+    cosThetaStar = cosThetaStarNeutrDaugh;
+  }
+  else if (mcPartDaughterNeutrDaugh->GetPdgCode() == -fPDGneutrDaughForMC){  // neutral particle bar
+    AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
+    cosThetaStar = cosThetaStarNeutrDaughBar;
+  }
+  else{
+    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()));
+    delete decay;
+    return vectorMC;
+  }
+  if (cosThetaStar < -1 || cosThetaStar > 1) {
+    AliWarning(Form("Invalid value for cosine Theta star %f, returning", cosThetaStar));
+    delete decay;
+    return bGenValues;
+  }
        
-       Double_t vectorD0[2] ={0.,0.};
-
-       // evaluate the correct cascade
-       if (!EvaluateIfD0toKpi(mcPartDaughterD0,vectorD0)) {
-         AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
-         delete decay;
-         return bGenValues;  
-       }       
-
-       //ct
-       Double_t cT = decay->Ct(421);
-       // get the pT of the daughters
-       Double_t pTD0 = 0.;
-       Double_t pTpi = 0.;
+  Double_t vectorNeutrDaugh[2] = {0.,0.};
+
+  // evaluate the correct cascade
+  if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
+    AliDebug(2, "Error! the Neutral Daughter MC doesn't have correct daughters!!");
+    delete decay;
+    return bGenValues;  
+  }    
+
+  //ct
+  Double_t cT = decay->Ct(fPDGneutrDaugh);
+  // get the pT of the daughters
+  Double_t pTNeutrDaugh= 0.;
+  Double_t pTBachelor = 0.;
        
-       if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
-               pTD0 = mcPartDaughterD0->Pt();
-               pTpi = mcPartDaughterPis->Pt();
-       }
+  if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
+    pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
+    pTBachelor = mcPartDaughterBachelor->Pt();
+  }
 
-       
-       switch (fConfiguration){
-       case AliCFTaskVertexingHF::kSnail:
-               vectorMC[0] = fmcPartCandidate->Pt();
-               vectorMC[1] = fmcPartCandidate->Y() ;
-               vectorMC[2] = cosThetaStar ;
-               vectorMC[3] = vectorD0[0]; 
-               vectorMC[4] = vectorD0[1];
-               vectorMC[5] = cT*1.E4 ;  // in micron
-               vectorMC[6] = 0.;   // dummy value, meaningless in MC
-               vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
-               vectorMC[8] = 1.01;    // dummy value, meaningless in MC
-               vectorMC[9] = fmcPartCandidate->Phi(); 
-               vectorMC[10] = fzMCVertex;    // z of reconstructed of primary vertex
-               vectorMC[11] = fCentValue; // reconstructed centrality
-               vectorMC[12] = 1.;           // always filling with 1 at MC level 
-               vectorMC[13] = 1.01; // dummy value for cosPointingXY  multiplicity
-               vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
-               vectorMC[15] = fMultiplicity; // reconstructed multiplicity
-               break;
-       case AliCFTaskVertexingHF::kCheetah:
-               vectorMC[0] = fmcPartCandidate->Pt();
-               vectorMC[1] = fmcPartCandidate->Y() ;
-               vectorMC[2] = cT*1.E4; // in micron
-               vectorMC[3] = fmcPartCandidate->Phi();
-               vectorMC[4] = fzMCVertex;
-               vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
-               vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level 
-               vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
-               break;
-       }
-
-       delete decay;
-       bGenValues = kTRUE;
-       return bGenValues;
+  AliDebug(3, Form("The candidate has pt = %f, y = %f", fmcPartCandidate->Pt(), fmcPartCandidate->Y()));
+
+  switch (fConfiguration){
+  case AliCFTaskVertexingHF::kSnail:
+    vectorMC[0] = fmcPartCandidate->Pt();
+    vectorMC[1] = fmcPartCandidate->Y() ;
+    vectorMC[2] = cosThetaStar ;
+    vectorMC[3] = vectorNeutrDaugh[0];
+    vectorMC[4] = vectorNeutrDaugh[1];
+    vectorMC[5] = cT*1.E4 ;  // in micron
+    vectorMC[6] = 0.;   // dummy value, meaningless in MC
+    vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
+    vectorMC[8] = 1.01;    // dummy value, meaningless in MC
+    vectorMC[9] = fmcPartCandidate->Phi(); 
+    vectorMC[10] = fzMCVertex;    // z of reconstructed of primary vertex
+    vectorMC[11] = fCentValue; // reconstructed centrality
+    vectorMC[12] = 1.;           // always filling with 1 at MC level 
+    vectorMC[13] = 1.01; // dummy value for cosPointingXY  multiplicity
+    vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
+    vectorMC[15] = fMultiplicity; // reconstructed multiplicity
+    break;
+  case AliCFTaskVertexingHF::kCheetah:
+    vectorMC[0] = fmcPartCandidate->Pt();
+    vectorMC[1] = fmcPartCandidate->Y() ;
+    vectorMC[2] = cT*1.E4; // in micron
+    vectorMC[3] = fmcPartCandidate->Phi();
+    vectorMC[4] = fzMCVertex;
+    vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
+    vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level 
+    vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
+    break;
+  }
+
+  delete decay;
+  bGenValues = kTRUE;
+  return bGenValues;
 }
 //____________________________________________
 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 { 
   // read the variables for the container
 
-  Bool_t bFillRecoValues=kFALSE;
-  
-  //Get the D* and the D0 from D*
-  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
-  AliAODRecoDecayHF2Prong* d0toKpi = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
+  Bool_t bFillRecoValues = kFALSE;
   
+  //Get the cascade and the neutral particle from it
+  AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
+  AliAODRecoDecay* neutrDaugh = NULL; 
+  if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
+  else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
+  else {
+    return kFALSE;
+  }
 
-  //if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi has primary vtx\n");
+  //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
   //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
 
-  Double_t pt =  dstarD0pi->Pt();
-  Double_t rapidity =  dstarD0pi->YDstar();
-  Double_t invMass=0.;
+  Double_t pt =  cascade->Pt();
+  Double_t rapidity =  cascade->Y(fPDGcascade);
+  Double_t invMass = 0.;
   Double_t cosThetaStar = 9999.;
-  Double_t pTpi = 0.;
-  Double_t pTK = 0.;
-  Double_t dca = d0toKpi->GetDCA();
-  Double_t d0pi = 0.;
-  Double_t d0K = 0.;
-  Double_t d0xd0 = d0toKpi->Prodd0d0();
-  Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
-  Double_t phi = dstarD0pi->Phi();
-  Double_t cosPointingAngleXY = d0toKpi->CosPointingAngleXY();
-  Double_t normDecayLengthXY = d0toKpi->NormalizedDecayLengthXY();
+  Double_t pTneutrDaughPos = 0.;
+  Double_t pTneutrDaughNeg = 0.;
+  Double_t dca = neutrDaugh->GetDCA();
+  Double_t d0neutrDaughPos = 0.;
+  Double_t d0neutrDaughNeg = 0.;
+  Double_t d0xd0 = neutrDaugh->Prodd0d0();
+  Double_t cosPointingAngle = neutrDaugh->CosPointingAngle(fPrimVtx);
+  Double_t phi = cascade->Phi();
+  Double_t cosPointingAngleXY = neutrDaugh->CosPointingAngleXY(fPrimVtx);
+  Double_t normDecayLengthXY = neutrDaugh->NormalizedDecayLengthXY(fPrimVtx);
 
   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
  
+  UInt_t pdgDaughCascade[2] = {fPDGbachelor, fPDGneutrDaugh};    // bachelor is first daughter of cascade
+  UInt_t pdgDaughBarCascade[2] = {fPDGneutrDaugh, fPDGbachelor}; // bachelor is second daughter in case of a cascade-bar 
+
   if (pdgCode > 0){
-    cosThetaStar = d0toKpi->CosThetaStarD0();
-    pTpi = d0toKpi->PtProng(0);
-    pTK = d0toKpi->PtProng(1);
-    d0pi = d0toKpi->Getd0Prong(0);
-    d0K = d0toKpi->Getd0Prong(1);
-    invMass=d0toKpi->InvMassD0();
+    cosThetaStar = neutrDaugh->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
+    pTneutrDaughPos = neutrDaugh->PtProng(0);
+    pTneutrDaughNeg = neutrDaugh->PtProng(1);
+    d0neutrDaughPos = neutrDaugh->Getd0Prong(0);
+    d0neutrDaughNeg = neutrDaugh->Getd0Prong(1);
+    invMass = neutrDaugh->InvMass(2, pdgDaughCascade);
   }
   else {
-    cosThetaStar = d0toKpi->CosThetaStarD0bar();
-    pTpi = d0toKpi->PtProng(1);
-    pTK = d0toKpi->PtProng(0);
-    d0pi = d0toKpi->Getd0Prong(1);
-    d0K = d0toKpi->Getd0Prong(0);
-    invMass= d0toKpi->InvMassD0bar();
+    cosThetaStar = neutrDaugh->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
+    pTneutrDaughPos = neutrDaugh->PtProng(1);
+    pTneutrDaughNeg = neutrDaugh->PtProng(0);
+    d0neutrDaughPos = neutrDaugh->Getd0Prong(1);
+    d0neutrDaughNeg = neutrDaugh->Getd0Prong(0);
+    invMass = neutrDaugh->InvMass(2, pdgDaughBarCascade);
   }
   
-  Double_t cT = d0toKpi->CtD0();
+  Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
   
-       switch (fConfiguration){
-       case AliCFTaskVertexingHF::kSnail:
-               vectorReco[0] = pt;
-               vectorReco[1] = rapidity;
-               vectorReco[2] = cosThetaStar;
-               vectorReco[3] = pTpi;
-               vectorReco[4] = pTK;
-               vectorReco[5] = cT*1.E4;  // in micron
-               vectorReco[6] = dca*1.E4;  // in micron
-               vectorReco[7] = d0xd0*1.E8;  // in micron^2
-               vectorReco[8] = cosPointingAngle;  // in micron
-               vectorReco[9] = phi;  
-               vectorReco[10] = fzPrimVertex;    // z of reconstructed of primary vertex
-               vectorReco[11] = fCentValue;
-               vectorReco[12] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
-               vectorReco[13] = cosPointingAngleXY; 
-               vectorReco[14] = normDecayLengthXY; // in cm
-               vectorReco[15] = fMultiplicity; // reconstructed multiplicity
-               break;
-       case AliCFTaskVertexingHF::kCheetah:
-               vectorReco[0] = pt;
-               vectorReco[1] = rapidity ;
-               vectorReco[2] = cT*1.E4; // in micron
-               vectorReco[3] = phi; 
-               vectorReco[4] = fzPrimVertex;
-               vectorReco[5] = fCentValue;   
-               vectorReco[6] = fFake ; 
-               vectorReco[7] = fMultiplicity;  
-               break;
-       }
-
-       bFillRecoValues = kTRUE;
-               
+  switch (fConfiguration){
+  case AliCFTaskVertexingHF::kSnail:
+    vectorReco[0] = pt;
+    vectorReco[1] = rapidity;
+    vectorReco[2] = cosThetaStar;
+    vectorReco[3] = pTneutrDaughPos;
+    vectorReco[4] = pTneutrDaughNeg;
+    vectorReco[5] = cT*1.E4;  // in micron
+    vectorReco[6] = dca*1.E4;  // in micron
+    vectorReco[7] = d0xd0*1.E8;  // in micron^2
+    vectorReco[8] = cosPointingAngle;  // in micron
+    vectorReco[9] = phi;
+    vectorReco[10] = fzPrimVertex;    // z of reconstructed of primary vertex
+    vectorReco[11] = fCentValue;
+    vectorReco[12] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
+    vectorReco[13] = cosPointingAngleXY;
+    vectorReco[14] = normDecayLengthXY; // in cm
+    vectorReco[15] = fMultiplicity; // reconstructed multiplicity
+    break;
+  case AliCFTaskVertexingHF::kCheetah:
+    vectorReco[0] = pt;
+    vectorReco[1] = rapidity;
+    vectorReco[2] = cT*1.E4; // in micron
+    vectorReco[3] = phi;
+    vectorReco[4] = fzPrimVertex;
+    vectorReco[5] = fCentValue;
+    vectorReco[6] = fFake;
+    vectorReco[7] = fMultiplicity; 
+    break;
+  }
+
+  bFillRecoValues = kTRUE;
+
   return bFillRecoValues;
 }
 
@@ -354,7 +389,6 @@ Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
   // check the required decay channel
 
   Bool_t checkCD = kFALSE;
-
   
   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
@@ -366,19 +400,32 @@ Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
     return checkCD;
   }
 
-  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
-       TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
-      !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
-       TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
-    AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
+  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGneutrDaughForMC &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGbachelor) && 
+      !(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGbachelor &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGneutrDaughForMC)) {
+    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()));
     return checkCD;  
   }
   
-  Double_t vectorD0[2] ={0.,0.};
+  // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
+  AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
+
+  // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the 
+  // charge of teh daughters to decide which is which
+  AliDebug(3, Form("Charge0 = %d, Charge1 = %d", mcPartDaughter0->Charge()/3, mcPartDaughter1->Charge()/3));
+  if (mcPartDaughter0->Charge()/3 != 0){
+    mcPartDaughterNeutrDaugh = mcPartDaughter1;
+  }
+  else {
+    mcPartDaughterNeutrDaugh = mcPartDaughter0;
+  }
 
-  // D* is a cascade ...evaluate the correct cascade
-  if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
-    AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+  Double_t vectorNeutrDaugh[2] ={0., 0.};
+
+  // We are looking at a cascade ...evaluate the correct cascade
+  if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
+    AliDebug(2, "Error! the NeutrDaugh MC doesn't have correct daughters!!");
     return checkCD;  
   }
    
@@ -388,49 +435,67 @@ Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
 }
 
 //__________________________________________
-Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
+Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
 {  
   //
-  // chack wether D0 is decaing into kpi
+  // check wether D0 is decaing into kpi
   //
   
   Bool_t isHadronic = kFALSE;
+  AliDebug(2, Form("neutralDaugh = %p, pdg = %d", neutralDaugh, neutralDaugh->GetPdgCode()));
+
+  if (fPDGcascade == 4122) {
+    Int_t labelresonanceDaugh = neutralDaugh->GetDaughter(0);
+    AliAODMCParticle* resonanceDaugh = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelresonanceDaugh));
+    if (!resonanceDaugh){
+      return kFALSE;
+    }
+    else {
+      if (TMath::Abs(resonanceDaugh->GetPdgCode()) != fPDGneutrDaugh){
+       return kFALSE;
+      }
+      else {
+       neutralDaugh = resonanceDaugh;
+      }
+    }
+  }
+
+  Int_t daughterNeutrDaugh0 = neutralDaugh->GetDaughter(0);
+  Int_t daughterNeutrDaugh1 = neutralDaugh->GetDaughter(1);
   
-  Int_t daughterD00 = neutralDaugh->GetDaughter(0);
-  Int_t daughterD01 = neutralDaugh->GetDaughter(1);
-  
-  AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
-  if (daughterD00 == 0 || daughterD01 == 0) {
+  AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
+  if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
     return isHadronic;  
   }
   
-  if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
+  Int_t numberOfExpectedDaughters = 2;
+  if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
     return isHadronic;  
   }
   
-  AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD00));
-  AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD01));
-  if (!mcPartDaughterD00 || !mcPartDaughterD01) {
+  AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
+  AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
+  if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping"); 
     return isHadronic;  
   }
   
-  if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
-       TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) && 
-      !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
-       TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
-    AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
+  if (!(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughPositive &&
+       TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughNegative) && 
+      !(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughNegative &&
+       TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughPositive)) {
+    AliDebug(2, "The neutral particle (MC) doesn't come from the required decay, skipping!!");
     return isHadronic;  
   }
   
-  Double_t sumPxDau=mcPartDaughterD00->Px()+mcPartDaughterD01->Px();
-  Double_t sumPyDau=mcPartDaughterD00->Py()+mcPartDaughterD01->Py();
-  Double_t sumPzDau=mcPartDaughterD00->Pz()+mcPartDaughterD01->Pz();
-  Double_t pxMother=neutralDaugh->Px();
-  Double_t pyMother=neutralDaugh->Py();
-  Double_t pzMother=neutralDaugh->Pz();
+  Double_t sumPxDau = mcPartDaughterNeutrDaugh0->Px()+mcPartDaughterNeutrDaugh1->Px();
+  Double_t sumPyDau = mcPartDaughterNeutrDaugh0->Py()+mcPartDaughterNeutrDaugh1->Py();
+  Double_t sumPzDau = mcPartDaughterNeutrDaugh0->Pz()+mcPartDaughterNeutrDaugh1->Pz();
+  Double_t pxMother = neutralDaugh->Px();
+  Double_t pyMother = neutralDaugh->Py();
+  Double_t pzMother = neutralDaugh->Pz();
   if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>0.00001 ||
      TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>0.00001 ||
      TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>0.00001){
@@ -438,23 +503,22 @@ Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh
     return isHadronic;  
   }
 
-  Double_t pTD0pi = 0;
-  Double_t pTD0K = 0;
-  
+  Double_t pTNeutrDaughPositive = 0;
+  Double_t pTNeutrDaughNegative = 0;
   
-  if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
-    pTD0pi = mcPartDaughterD00->Pt();
-    pTD0K = mcPartDaughterD01->Pt();
+  if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
+    pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
+    pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
   }
   else {
-    pTD0pi = mcPartDaughterD01->Pt();
-    pTD0K  = mcPartDaughterD00->Pt();
+    pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
+    pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
   }
   
   isHadronic = kTRUE;
   
-  vectorD0[0] = pTD0pi;
-  vectorD0[1] = pTD0K;
+  vectorNeutrDaugh[0] = pTNeutrDaughPositive;
+  vectorNeutrDaugh[1] = pTNeutrDaughNegative;
  
   return isHadronic;
 
@@ -464,17 +528,17 @@ Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh
 
 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
 {
-       //
-       // setting the pt cut to be used in the Acceptance steps (MC+Reco)
-       //
-
-       AliInfo("The 3rd element of the pt cut array will correspond to the cut applied to the soft pion - please check that it is correct");
-       if (fProngs>0){
-               for (Int_t iP=0; iP<fProngs; iP++){
-                       fPtAccCut[iP]=ptAccCut[iP];
-               }
-       }
-       return;
+  //
+  // setting the pt cut to be used in the Acceptance steps (MC+Reco)
+  //
+
+  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");
+  if (fProngs>0){
+    for (Int_t iP=0; iP<fProngs; iP++){
+      fPtAccCut[iP]=ptAccCut[iP];
+    }
+  }
+  return;
 }              
 
 
@@ -483,51 +547,53 @@ void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
 
 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
 {
-       //
-       // setting the eta cut to be used in the Acceptance steps (MC+Reco)
-       //
-
-       AliInfo("The 3rd element of the eta cut array will correspond to the cut applied to the soft pion - please check that it is correct");
-       if (fProngs>0){
-               for (Int_t iP=0; iP<fProngs; iP++){
-                       fEtaAccCut[iP]=etaAccCut[iP];
-               }
-       }
-       return;
-}      
+  //
+  // setting the eta cut to be used in the Acceptance steps (MC+Reco)
+  //
+
+  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");
+  if (fProngs>0){
+    for (Int_t iP=0; iP<fProngs; iP++){
+      fEtaAccCut[iP] = etaAccCut[iP];
+    }
+  }
+  return;
+}
 //___________________________________________________________
 
 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 {
-       //
-       // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
-       //
-
-       AliInfo("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");
-       if (fProngs>0){
-               for (Int_t iP=0; iP<fProngs; iP++){
-                       fPtAccCut[iP]=ptAccCut[iP];
-                       fEtaAccCut[iP]=etaAccCut[iP];
-               }
-       }
-       return;
-}              
+  //
+  // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
+  //
+
+  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");
+  if (fProngs>0){
+    for (Int_t iP=0; iP<fProngs; iP++){
+      fPtAccCut[iP]=ptAccCut[iP];
+      fEtaAccCut[iP]=etaAccCut[iP];
+    }
+  }
+  return;
+}
 
 //___________________________________________________________
 
 void AliCFVertexingHFCascade::SetAccCut()
 {
-       //
-       // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
-       //
-  
-  AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[2]));  // should be the soft pion...  
+  //
+  // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
+  //
+
+  Int_t bachelorPosition = 2;
+  if (fPDGcascade == 4122) bachelorPosition = 0;
+  AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[bachelorPosition]));  // should be the soft pion...  
   if(!mcPartDaughter) return;
   Int_t mother =  mcPartDaughter->GetMother();
   AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); 
   if(!mcMother) return;
 
-  if (TMath::Abs(mcPartDaughter->GetPdgCode())!= 211 || TMath::Abs(mcMother->GetPdgCode())!=413){
+  if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
     AliFatal("Apparently the soft pion is not in the third position, causing a crash!!");
   }             
   if (fProngs>0){
@@ -539,23 +605,26 @@ void AliCFVertexingHFCascade::SetAccCut()
     fEtaAccCut[2]=0.9;  // soft pion
   }
   return;
-}              
+}
 
 //_____________________________________________________________
 Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const 
 {
-       //
-       // getting eta of the prong - overload the mother class method
-       //
+  //
+  // getting eta of the prong - overload the mother class method
+  //
 
- if (fRecoCandidate){
+  if (fRecoCandidate){
    
-   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
-
-   Double_t etaProng =-9999;
-    if(iProng==0) etaProng =dstarD0pi->Get2Prong()->EtaProng(0);
-    if(iProng==1) etaProng =dstarD0pi->Get2Prong()->EtaProng(1);
-    if(iProng==2) etaProng =dstarD0pi->EtaProng(1);
+    AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
+
+    Double_t etaProng =-9999;
+    AliAODRecoDecay* neutrDaugh; 
+    if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
+    else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
+    if (iProng==0) etaProng = neutrDaugh->EtaProng(0);
+    if (iProng==1) etaProng = neutrDaugh->EtaProng(1);
+    if (iProng==2) etaProng = cascade->EtaProng(1);
     
     return etaProng;
     
@@ -565,22 +634,70 @@ Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const
 //_____________________________________________________________
 Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const 
 {
-       //
-       // getting pt of the prong
-       //
+  //
+  // getting pt of the prong
+  //
   
   if (fRecoCandidate){
 
-    AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
+    AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
     Double_t ptProng= -9999;
-    if(iProng==0) ptProng =dstarD0pi->Get2Prong()->PtProng(0);
-    if(iProng==1) ptProng =dstarD0pi->Get2Prong()->PtProng(1);
-    if(iProng==2) ptProng =dstarD0pi->PtProng(1);
+    AliAODRecoDecay* neutrDaugh; 
+    if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
+    else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
+    if (iProng == 0) ptProng = neutrDaugh->PtProng(0);
+    if (iProng == 1) ptProng = neutrDaugh->PtProng(1);
+    if (iProng == 2) ptProng = cascade->PtProng(1);
     
     // Double_t ptProng = fRecoCandidate->PtProng(iProng);  
     return ptProng;
     
   }
   return 999999;  
-  
+
+}
+//_____________________________________________________________
+Bool_t AliCFVertexingHFCascade::CheckAdditionalCuts(AliPIDResponse* pidResponse) const {
+
+  // function to check whether the candidate passes the additional cuts defined in the task to get the
+  // invariant mass spectra; these cuts are NOT pt-dependent
+
+  if (fPDGcascade == 4122){
+    // the method is implemented only in this case so far
+    AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
+    AliAODv0 * v0part = cascade->Getv0();
+    AliAODTrack *bachelor = (AliAODTrack*)cascade->GetBachelor();
+
+    Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
+    Double_t nSigmaTPCpr=-999.;
+    Double_t nSigmaTOFpr=-999.;
+    nSigmaTPCpr = pidResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
+    nSigmaTOFpr = pidResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
+    Double_t ptArm = v0part->PtArmV0();
+    Double_t invmassK0s = v0part->MassK0Short();
+    Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
+
+    Bool_t cutsForTMVA = ((nSigmaTOFpr > -800) && (nSigmaTOFpr < 3)) &&  
+      ((ptArm > 0.07) && (ptArm < 0.105)) &&
+      ((TMath::Abs(invmassK0s - mK0SPDG)) < 0.01);
+    
+    if (!fUseCutsForTMVA) cutsForTMVA = kTRUE;
+
+    Bool_t cutsForInvMassTask = !(onFlyV0) && 
+      (cascade->CosV0PointingAngle()>0.99) && 
+      (TMath::Abs(nSigmaTPCpr) <= 3) && 
+      (v0part->Getd0Prong(0) < 20) && 
+      (v0part->Getd0Prong(1) < 20);
+
+    if (cutsForTMVA && cutsForInvMassTask) {   
+      // K0Smass cut
+      // eta cut
+      // TOF PID cut
+      // Arm cut
+      return kTRUE;
+    }
+  }
+
+  return kFALSE;
+
 }
index 7162592..f0bac2b 100644 (file)
@@ -34,6 +34,7 @@ class TClonesArray;
 class AliCFVertexingHF;
 class AliESDtrack;
 class TDatabasePDG;
+class AliPIDResponse;
 
 class AliCFVertexingHFCascade : public AliCFVertexingHF{
  public:
@@ -49,7 +50,8 @@ class AliCFVertexingHFCascade : public AliCFVertexingHF{
   Bool_t CheckMCChannelDecay()const;
   
   Bool_t SetRecoCandidateParam(AliAODRecoDecayHF *recoCand);
-  Bool_t EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* VectorD0)const;
+  //Bool_t EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* VectorD0)const;
+  Bool_t EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* VectorD0)const;
 
   void SetPtAccCut(Float_t* ptAccCut);
   void SetEtaAccCut(Float_t* etaAccCut);
@@ -59,14 +61,45 @@ class AliCFVertexingHFCascade : public AliCFVertexingHF{
   Double_t GetEtaProng(Int_t iProng)const;
   Double_t GetPtProng(Int_t iProng) const;
 
+  void SetPDGcascade(Int_t pdg)    {fPDGcascade = pdg;}
+  void SetPDGbachelor(Int_t pdg)   {fPDGbachelor = pdg;}
+  void SetPDGneutrDaugh(Int_t pdg)         {fPDGneutrDaugh = pdg;}
+  void SetPDGneutrDaughForMC(Int_t pdg)         {fPDGneutrDaughForMC = pdg;}
+  void SetPDGneutrDaughPositive(Int_t pdg) {fPDGneutrDaughPositive = pdg;}
+  void SetPDGneutrDaughNegative(Int_t pdg) {fPDGneutrDaughNegative = pdg;}
+  void SetPrimaryVertex(AliAODVertex* vtx) {fPrimVtx = vtx;}
+
+  Int_t GetPDGcascade()    const {return fPDGcascade;}
+  Int_t GetPDGbachelor()   const {return fPDGbachelor;}
+  Int_t GetPDGneutrDaugh()         const {return fPDGneutrDaugh;}
+  Int_t GetPDGneutrDaughForMC()         const {return fPDGneutrDaughForMC;}
+  Int_t GetPDGneutrDaughPositive() const {return fPDGneutrDaughPositive;}
+  Int_t GetPDGneutrDaughNegative() const {return fPDGneutrDaughNegative;}
+  AliAODVertex* GetPrimaryVertex() const {return fPrimVtx;}
+
+  Bool_t CheckAdditionalCuts(AliPIDResponse* pidResponse) const;
+
+  void SetUseCutsForTMVA(Bool_t useCutsForTMVA) {fUseCutsForTMVA = useCutsForTMVA;}
+  Bool_t GetUseCutsForTMVA() const {return fUseCutsForTMVA;}
+
  protected:
   
   
  private:      
   AliCFVertexingHFCascade(const AliCFVertexingHFCascade& c);
   AliCFVertexingHFCascade& operator= (const AliCFVertexingHFCascade& other);
-  
-  ClassDef(AliCFVertexingHFCascade, 2); // CF class for D* and other cascades
+
+  Int_t fPDGcascade;   // pdg code of the cascade
+  Int_t fPDGbachelor;  // pdg code of the bachelor
+  Int_t fPDGneutrDaugh;        // pdg code of the V0
+  Int_t fPDGneutrDaughForMC;        // pdg code of the V0
+  Int_t fPDGneutrDaughPositive;  // pdg code of the positive daughter of the V0
+  Int_t fPDGneutrDaughNegative;  // pdg code of the negative daughter of the V0
+  AliAODVertex* fPrimVtx;        // primaryVertex
+  Bool_t fUseCutsForTMVA;        // flag to decide whether to use or not the preselection
+                                 // cuts of the TMVA when filling the CF
+
+  ClassDef(AliCFVertexingHFCascade, 3); // CF class for D* and other cascades
   
 };