[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFVertexingHFCascade.cxx

/**************************************************************************
 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------
// 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 "AliAODRecoDecayHF2Prong.h"
#include "AliAODMCParticle.h"
#include "AliAODEvent.h"
#include "TClonesArray.h"
#include "AliCFVertexingHF.h"
#include "AliESDtrack.h"
#include "TDatabasePDG.h"
#include "AliAODRecoCascadeHF.h"
#include "AliCFVertexingHFCascade.h"
#include "AliCFContainer.h"

ClassImp(AliCFVertexingHFCascade)


//_________________________________________
  AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
  AliCFVertexingHF(mcArray, originDselection)
{	
  // standard constructor

  SetNProngs(3);
}


//_____________________________________
AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
{
  // operator =
 
  if  (this != &c) {

    AliCFVertexingHF::operator=(c);
   
  }
  return *this;
}

//__________________________________________
Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
{
  // set the AliAODRecoDecay candidate
  
  Bool_t bSignAssoc = kFALSE;
 
  fRecoCandidate = recoCand;
  AliAODRecoCascadeHF* dstarD0pi = (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 nentries = fmcArray->GetEntriesFast();

  AliDebug(3,Form("nentries = %d\n", nentries));
 
  Int_t mcLabel =  dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray); 
  
  if (mcLabel < 0) return bSignAssoc;
  SetMCLabel(mcLabel);
  fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel)); 

  if (!fmcPartCandidate){
    AliDebug(3,"No part candidate");
    return bSignAssoc;
  }
 
  bSignAssoc = kTRUE;
  return bSignAssoc;
}

//______________________________________________
Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) 
{
  // 
  // 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));
	if (!mcPartDaughterD0 || !mcPartDaughterPis){
		if (!mcPartDaughterD0){
			AliError("Error while casting the MC D0, returning FALSE");
		} 
		if (!mcPartDaughterD0){
			AliError("Error while casting the MC soft pion, returning FALSE");
		} 
		return bGenValues;
	}

	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

	//D0 daughters
	Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
	Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);

	AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
	AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis

	if (!mcPartDaughter0 || !mcPartDaughter1){
		if (!mcPartDaughter0){
			AliError("Error while casting the MC daughter 0, returning FALSE");
		} 
		if (!mcPartDaughter1){
			AliError("Error while casting the MC daughter 1, returning FALSE");
		} 
		return bGenValues;
	}

	// 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);
	
	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...");
	  return vectorMC;
	}
	if (cosThetaStar < -1 || cosThetaStar > 1) {
	  AliWarning("Invalid value for cosine Theta star, returning");
	  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!!");
	  return bGenValues;  
	}	

	//ct
	Double_t cT = decay->Ct(421);
	// get the pT of the daughters
	Double_t pTD0 = 0.;
	Double_t pTpi = 0.;
	
	if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
		pTD0 = mcPartDaughterD0->Pt();
		pTpi = mcPartDaughterPis->Pt();
	}

	
	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] = 0.;   // dummy value, meaningless in MC, in micron
	vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
	vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
	vectorMC[10] = 1.01;    // dummy value, meaningless in MC
	vectorMC[11] = fmcPartCandidate->Phi(); 
	vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
       

	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();
  

  if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi 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 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();

  Int_t pdgCode = fmcPartCandidate->GetPdgCode();
 
  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();
  }
  else {
    cosThetaStar = d0toKpi->CosThetaStarD0bar();
    pTpi = d0toKpi->PtProng(1);
    pTK = d0toKpi->PtProng(0);
    d0pi = d0toKpi->Getd0Prong(1);
    d0K = d0toKpi->Getd0Prong(0);
    invMass= d0toKpi->InvMassD0bar();
  }
  
  Double_t cT = d0toKpi->CtD0();
  
  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] = d0pi*1.E4;  // in micron
  vectorReco[8] = d0K*1.E4;  // in micron
  vectorReco[9] = d0xd0*1.E8;  // in micron^2
  vectorReco[10] = cosPointingAngle;  // in micron
  vectorReco[11] = phi;  
  vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex

  bFillRecoValues = kTRUE;

  return bFillRecoValues;
}


//_____________________________________________________________
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);
  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));

  if (!mcPartDaughter0 || !mcPartDaughter1) {
    AliDebug (2,"Problems in the MC Daughters\n");
    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!!");
    return checkCD;  
  }
  
  Double_t vectorD0[2] ={0.,0.};

  // D* is a cascade ...evaluate the correct cascade
  if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
    AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
    return checkCD;  
  }
   
  checkCD = kTRUE;
  return checkCD;
  
}

//__________________________________________
Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
{  
  //
  // chack wether D0 is decaing into kpi
  //
  
  Bool_t isHadronic = kFALSE;
  
  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, "Error! the D0 MC doesn't have correct daughters!!");
    return isHadronic;  
  }
  
  if (TMath::Abs(daughterD01 - daughterD00) != 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) {
    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!!");
    return isHadronic;  
  }
  
  Double_t pTD0pi = 0;
  Double_t pTD0K = 0;
  
  
  if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
    pTD0pi = mcPartDaughterD00->Pt();
    pTD0K = mcPartDaughterD01->Pt();
  }
  else {
    pTD0pi = mcPartDaughterD01->Pt();
    pTD0K  = mcPartDaughterD00->Pt();
  }
  
  isHadronic = kTRUE;
  
  vectorD0[0] = pTD0pi;
  vectorD0[1] = pTD0K;
 
  return isHadronic;

}
Commit	Line	Data
a9f238cf	1	/**************************************************************************
	2	* Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------------
	17	// Class for HF corrections as a function of many variables and steps
	18	// For D* and other cascades
	19	//
	20	// Author : A.Grelli a.grelli@uu.nl UTECHT
	21	//-----------------------------------------------------------------------
	22
	23	#include "AliAODRecoDecayHF2Prong.h"
	24	#include "AliAODMCParticle.h"
	25	#include "AliAODEvent.h"
	26	#include "TClonesArray.h"
	27	#include "AliCFVertexingHF.h"
	28	#include "AliESDtrack.h"
	29	#include "TDatabasePDG.h"
	30	#include "AliAODRecoCascadeHF.h"
	31	#include "AliCFVertexingHFCascade.h"
	32	#include "AliCFContainer.h"
	33
	34	ClassImp(AliCFVertexingHFCascade)
	35
	36
	37	//_________________________________________
	38	AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
	39	AliCFVertexingHF(mcArray, originDselection)
	40	{
	41	// standard constructor
	42
	43	SetNProngs(3);
	44	}
	45
	46
	47	//_____________________________________
	48	AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
	49	{
	50	// operator =
	51
	52	if (this != &c) {
	53
	54	AliCFVertexingHF::operator=(c);
	55
	56	}
	57	return *this;
	58	}
	59
	60	//__________________________________________
	61	Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
	62	{
	63	// set the AliAODRecoDecay candidate
	64
65	Bool_t bSignAssoc = kFALSE;
66
67	fRecoCandidate = recoCand;
68	AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)recoCand;
69
70	if (!fRecoCandidate) {
71	AliError("fRecoCandidate not found, problem in assignement\n");
72	return bSignAssoc;
73	}
74
75	if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
76	//if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
77
78	//Int_t pdgCand = 413;
79
80	Int_t pdgDgDStartoD0pi[2]={421,211};
81	Int_t pdgDgD0toKpi[2]={321,211};
82
83	Int_t nentries = fmcArray->GetEntriesFast();
84
85	AliDebug(3,Form("nentries = %d\n", nentries));
86
87	Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
88
89	if (mcLabel < 0) return bSignAssoc;
90	SetMCLabel(mcLabel);
91	fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
92
93	if (!fmcPartCandidate){
94	AliDebug(3,"No part candidate");
95	return bSignAssoc;
96	}
97
98	bSignAssoc = kTRUE;
99	return bSignAssoc;
100	}
101
102	//______________________________________________
103	Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
104	{
105	//
106	// collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
107	//
108
109	Bool_t bGenValues = kFALSE;
110
111
112	Int_t daughter0ds = fmcPartCandidate->GetDaughter(0);
113	Int_t daughter1ds = fmcPartCandidate->GetDaughter(1);
114
115	//the D0
116	AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0ds));
117	AliAODMCParticle* mcPartDaughterPis = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1ds));
2c16db30	118	if (!mcPartDaughterD0 \|\| !mcPartDaughterPis){
	119	if (!mcPartDaughterD0){
	120	AliError("Error while casting the MC D0, returning FALSE");
	121	}
	122	if (!mcPartDaughterD0){
	123	AliError("Error while casting the MC soft pion, returning FALSE");
	124	}
	125	return bGenValues;
	126	}
a9f238cf	127
	128	Double_t vtx1[3] = {0,0,0}; // primary vertex
	129	Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
	130	Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
	131	fmcPartCandidate->XvYvZv(vtx1); // cm
	132
	133	//D0 daughters
	134	Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
	135	Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);
	136
	137	AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
	138	AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis
	139
2c16db30	140	if (!mcPartDaughter0 \|\| !mcPartDaughter1){
	141	if (!mcPartDaughter0){
	142	AliError("Error while casting the MC daughter 0, returning FALSE");
	143	}
	144	if (!mcPartDaughter1){
	145	AliError("Error while casting the MC daughter 1, returning FALSE");
	146	}
	147	return bGenValues;
	148	}
	149
a9f238cf	150	// getting vertex from daughters
	151	mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
	152	mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
	153	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) {
	154	AliError("Daughters have different secondary vertex, skipping the track");
	155	return bGenValues;
	156	}
	157
	158	Int_t nprongs = 2;
	159	Short_t charge = 0;
	160	// always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
	161	AliAODMCParticle* positiveDaugh = mcPartDaughter0;
	162	AliAODMCParticle* negativeDaugh = mcPartDaughter1;
	163	if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
	164	// inverting in case the positive daughter is the second one
	165	positiveDaugh = mcPartDaughter1;
	166	negativeDaugh = mcPartDaughter0;
	167	}
	168	// getting the momentum from the daughters
	169	Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
	170	Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
	171	Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
	172
	173	Double_t d0[2] = {0.,0.};
	174
	175	AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
	176
	177	Double_t cosThetaStar = 0.;
	178	Double_t cosThetaStarD0 = 0.;
	179	Double_t cosThetaStarD0bar = 0.;
	180	cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
	181	cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
	182	if (mcPartDaughterD0->GetPdgCode() == 421){ // D0
	183	AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
	184	cosThetaStar = cosThetaStarD0;
	185	}
	186	else if (mcPartDaughterD0->GetPdgCode() == -421){ // D0bar{
	187	AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
	188	cosThetaStar = cosThetaStarD0bar;
	189	}
	190	else{
	191	AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
	192	return vectorMC;
	193	}
	194	if (cosThetaStar < -1 \|\| cosThetaStar > 1) {
	195	AliWarning("Invalid value for cosine Theta star, returning");
	196	return bGenValues;
	197	}
	198
	199	Double_t vectorD0[2] ={0.,0.};
	200
	201	// evaluate the correct cascade
	202	if (!EvaluateIfD0toKpi(mcPartDaughterD0,vectorD0)) {
	203	AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
	204	return bGenValues;
	205	}
	206
	207	//ct
	208	Double_t cT = decay->Ct(421);
	209	// get the pT of the daughters
	210	Double_t pTD0 = 0.;
	211	Double_t pTpi = 0.;
	212
	213	if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
214	pTD0 = mcPartDaughterD0->Pt();
215	pTpi = mcPartDaughterPis->Pt();
216	}
217
218
219	vectorMC[0] = fmcPartCandidate->Pt();
220	vectorMC[1] = fmcPartCandidate->Y() ;
221	vectorMC[2] = cosThetaStar ;
222	vectorMC[3] = vectorD0[0];
223	vectorMC[4] = vectorD0[1];
224	vectorMC[5] = cT*1.E4 ; // in micron
225	vectorMC[6] = 0.; // dummy value, meaningless in MC
226	vectorMC[7] = 0.; // dummy value, meaningless in MC, in micron
227	vectorMC[8] = 0.; // dummy value, meaningless in MC, in micron
228	vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
229	vectorMC[10] = 1.01; // dummy value, meaningless in MC
230	vectorMC[11] = fmcPartCandidate->Phi();
231	vectorMC[12] = fzMCVertex; // z of reconstructed of primary vertex
232
233
234	delete decay;
235	bGenValues = kTRUE;
236	return bGenValues;
237	}
238	//____________________________________________
239	Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
240	{
241	// read the variables for the container
242
243	Bool_t bFillRecoValues=kFALSE;
244
245	//Get the D* and the D0 from D*
246	AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
247	AliAODRecoDecayHF2Prong* d0toKpi = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
248
249
250	if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi has primary vtx\n");
251	if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
252
253	Double_t pt = dstarD0pi->Pt();
254	Double_t rapidity = dstarD0pi->YDstar();
255	Double_t invMass=0.;
256	Double_t cosThetaStar = 9999.;
257	Double_t pTpi = 0.;
258	Double_t pTK = 0.;
259	Double_t dca = d0toKpi->GetDCA();
260	Double_t d0pi = 0.;
261	Double_t d0K = 0.;
262	Double_t d0xd0 = d0toKpi->Prodd0d0();
263	Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
264	Double_t phi = dstarD0pi->Phi();
265
266	Int_t pdgCode = fmcPartCandidate->GetPdgCode();
267
268	if (pdgCode > 0){
269	cosThetaStar = d0toKpi->CosThetaStarD0();
270	pTpi = d0toKpi->PtProng(0);
271	pTK = d0toKpi->PtProng(1);
272	d0pi = d0toKpi->Getd0Prong(0);
273	d0K = d0toKpi->Getd0Prong(1);
274	invMass=d0toKpi->InvMassD0();
275	}
276	else {
277	cosThetaStar = d0toKpi->CosThetaStarD0bar();
278	pTpi = d0toKpi->PtProng(1);
279	pTK = d0toKpi->PtProng(0);
280	d0pi = d0toKpi->Getd0Prong(1);
281	d0K = d0toKpi->Getd0Prong(0);
282	invMass= d0toKpi->InvMassD0bar();
283	}
284
285	Double_t cT = d0toKpi->CtD0();
286
287	vectorReco[0] = pt;
288	vectorReco[1] = rapidity;
289	vectorReco[2] = cosThetaStar;
290	vectorReco[3] = pTpi;
291	vectorReco[4] = pTK;
292	vectorReco[5] = cT*1.E4; // in micron
293	vectorReco[6] = dca*1.E4; // in micron
294	vectorReco[7] = d0pi*1.E4; // in micron
295	vectorReco[8] = d0K*1.E4; // in micron
296	vectorReco[9] = d0xd0*1.E8; // in micron^2
297	vectorReco[10] = cosPointingAngle; // in micron
298	vectorReco[11] = phi;
299	vectorReco[12] = fzPrimVertex; // z of reconstructed of primary vertex
300
301	bFillRecoValues = kTRUE;
302
303	return bFillRecoValues;
304	}
305
306
307	//_____________________________________________________________
308	Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
309	{
310	// check the required decay channel
311
312	Bool_t checkCD = kFALSE;
313
314
315	Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
316	Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
317	AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
318	AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
319
320	if (!mcPartDaughter0 \|\| !mcPartDaughter1) {
321	AliDebug (2,"Problems in the MC Daughters\n");
322	return checkCD;
323	}
324
325	if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
326	TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
327	!(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
328	TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
329	AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
330	return checkCD;
331	}
332
333	Double_t vectorD0[2] ={0.,0.};
334
335	// D* is a cascade ...evaluate the correct cascade
336	if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
337	AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
338	return checkCD;
339	}
340
341	checkCD = kTRUE;
342	return checkCD;
343
344	}
345
346	//__________________________________________
347	Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
348	{
349	//
350	// chack wether D0 is decaing into kpi
351	//
352
353	Bool_t isHadronic = kFALSE;
354
355	Int_t daughterD00 = neutralDaugh->GetDaughter(0);
356	Int_t daughterD01 = neutralDaugh->GetDaughter(1);
357
358	AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
359	if (daughterD00 == 0 \|\| daughterD01 == 0) {
360	AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
361	return isHadronic;
362	}
363
364	if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
365	AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
366	return isHadronic;
367	}
368
369	AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD00));
370	AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD01));
371	if (!mcPartDaughterD00 \|\| !mcPartDaughterD01) {
372	AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
373	return isHadronic;
374	}
375
376	if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
377	TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
378	!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
379	TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
380	AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
381	return isHadronic;
382	}
383
384	Double_t pTD0pi = 0;
385	Double_t pTD0K = 0;
386
387
388	if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
389	pTD0pi = mcPartDaughterD00->Pt();
390	pTD0K = mcPartDaughterD01->Pt();
391	}
392	else {
393	pTD0pi = mcPartDaughterD01->Pt();
394	pTD0K = mcPartDaughterD00->Pt();
395	}
396
397	isHadronic = kTRUE;
398
399	vectorD0[0] = pTD0pi;
400	vectorD0[1] = pTD0K;
401
402	return isHadronic;
403
404	}
405
406
407
408