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

/**************************************************************************
 * Copyright(c) 2007-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.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////
//                                                               //
// Class to compute variables for correction framework           //  
// for 3-body decays of D mesons (D+, Ds, Lc)                    //
// in bins of cut variables                                      //
// Origin:       Francesco Prino (prino@to.infn.it)              //
//               Renu Bala       (bala@to.infn.it)               //
//               Davide Caffarri (cafarri@pd.infn.it)            //
///////////////////////////////////////////////////////////////////

#include "AliAODMCParticle.h"
#include "AliAODEvent.h"
#include "TClonesArray.h"
#include "AliCFVertexingHF.h"
#include "AliESDtrack.h"
#include "TDatabasePDG.h"

#include "AliCFVertexingHF3Prong.h"
#include "AliCFContainer.h"

ClassImp(AliCFVertexingHF3Prong)


//_________________________________________
AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay):
AliCFVertexingHF(),
  fDecay(decay)
 {
  // 
  SetNProngs(3);

  fPtAccCut=new Float_t[fProngs];
  fEtaAccCut=new Float_t[fProngs];
  for(Int_t iP=0; iP<fProngs; iP++){
	  fPtAccCut[iP]=0.1;
	  fEtaAccCut[iP]=0.9;
  }

}
//_________________________________________
AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
  AliCFVertexingHF(mcArray, originDselection),
  fDecay(decay)
 
{
  //
  SetNProngs(3);
  fPtAccCut=new Float_t[fProngs];
  fEtaAccCut=new Float_t[fProngs];
  for(Int_t iP=0; iP<fProngs; iP++){
	  fPtAccCut[iP]=0.1;
	  fEtaAccCut[iP]=0.9;
  }
}


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

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

//__________________________________________
Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
  // Checks if candidate is signal and D meson is present in MC array
  
  Bool_t bSignAssoc = kFALSE;
  fRecoCandidate = recoCand;

  if (!fRecoCandidate) {
    AliError("fRecoCandidate not found, problem in assignement\n");
    return bSignAssoc;
  }
  
  Int_t pdgCand = -1;
  Int_t pdgDaughter[3]={-1,-1,-1};
  if(fDecay==kDplustoKpipi){
    pdgCand=411;
    pdgDaughter[0]=321;
    pdgDaughter[1]=211;
    pdgDaughter[2]=211;
  }else if(fDecay==kDstoKKpi){
    pdgCand=431;
    pdgDaughter[0]=321;
    pdgDaughter[1]=321;
    pdgDaughter[2]=211;
  }else if(fDecay==kLctopKpi){
	  pdgCand=4122;
	  pdgDaughter[0]=2212;
	  pdgDaughter[1]=321;
	  pdgDaughter[2]=211;	  
  }else{
    AliError("WRONG DECAY SETTING");
    return bSignAssoc;    
  }

  Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);  
  if (mcLabel == -1) return bSignAssoc;

  if (fRecoCandidate->NumberOfFakeDaughters()>0){
	  fFake = 0;    // fake candidate
	  if (fFakeSelection==1) return bSignAssoc;
  }
  if (fRecoCandidate->NumberOfFakeDaughters()==0){
	  fFake = 2;    // non-fake candidate
	  if (fFakeSelection==2) 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 AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
  // 
  // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
  // pt_D
  // y_D
  // phi_D
  // ctau
  // cos point
  // pt_1
  // pt_2
  // pt_3
  // d0_1
  // d0_2
  // d0_3
  // zPrimVert
  // centrality

  Bool_t bGenValues = kFALSE;
		
  Int_t pdgCand = -1;
  if(fDecay==kDplustoKpipi){
    pdgCand=411;
  }else if(fDecay==kDstoKKpi){
    pdgCand=431;
  }else if(fDecay==kLctopKpi){
	pdgCand=4122;
  }else{
    AliError("WRONG DECAY SETTING");
    return bGenValues;
  }

  Double_t vertD[3] = {0,0,0};   // D origin
  fmcPartCandidate->XvYvZv(vertD);  // cm

  Int_t nprongs = 3;
  Int_t daughter[3];
  Short_t charge = fmcPartCandidate->Charge();

  // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
  // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
  Int_t index=0;
  Int_t nDauLS=0;
  Int_t nDauOS=0;


  Int_t nDau=fmcPartCandidate->GetNDaughters();
  Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0); 
  if(nDau==3){
    for(Int_t iDau=0; iDau<3; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return bGenValues;
      }
      Short_t signDau=part->Charge();
      if(signDau==charge){
	nDauLS++;
	daughter[index] = ind;
	index=2;
      }else{
	daughter[1] = ind;
	nDauOS++;
      }
    }
  }else if(nDau==2){
    for(Int_t iDau=0; iDau<2; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return bGenValues;
      }
      Int_t pdgCode=TMath::Abs(part->GetPdgCode());
      if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
	Short_t signDau=part->Charge();
	if(signDau==charge){
	  nDauLS++;
	  daughter[index] = ind;
	  index=2;
	}else{
	  daughter[1] = ind;
	  nDauOS++;
	}
      }else{
	Int_t nDauRes=part->GetNDaughters();
	if(nDauRes!=2){
	  AliError("Wrong resonant decay");
	  return bGenValues;
	}
	Int_t labelFirstDauRes = part->GetDaughter(0); 	
	for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	  Int_t indDR = labelFirstDauRes+iDauRes;
	  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
	  if(!partDR){
	    AliError("Daughter particle not found in MC array");
	    return bGenValues;
	  }
	  Short_t signDau=partDR->Charge();
	  if(signDau==charge){
	    nDauLS++;
	    daughter[index] = ind;
	    index=2;
	  }else{
	    daughter[1] = ind;
	    nDauOS++;
	  }
	}
      }
    }
  }else{
    AliError(Form("Wrong number of daughters %d",nDau));
    return bGenValues;
  }

  if(nDauLS!=2 || nDauOS!=1){
    AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
    return bGenValues;
  }
  if(daughter[0]>daughter[2]){
    Int_t tmp=daughter[0];
    daughter[0]=daughter[2];
    daughter[2]=tmp;
  }
  
  // getting the momentum from the daughters and decay vertex
  Double_t px[3],py[3],pz[3],pt[3];
  Double_t vertDec[3] = {0,0,0};   // decay vertex		
  for(Int_t iDau=0; iDau<3; iDau++){
    AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
    if(!part){
      AliError("Daughter particle not found in MC array");
      return bGenValues;
    }
    px[iDau]=part->Px();
    py[iDau]=part->Py();
    pz[iDau]=part->Pz();
    pt[iDau]=part->Pt();
    if(iDau==0) part->XvYvZv(vertDec);
  }

  Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
  
  AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
  Double_t cT = decay->Ct(pdgCand);
		
  vectorMC[0] = fmcPartCandidate->Pt();
  vectorMC[1] = fmcPartCandidate->Y() ;
  vectorMC[2] = fmcPartCandidate->Phi();
  vectorMC[3] = cT*1.E4 ;  // in micron
  vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
  vectorMC[5] = pt[0];
  vectorMC[6] = pt[1];
  vectorMC[7] = pt[2];
  vectorMC[8] = 0.;   // imppar0, dummy value, meaningless in MC
  vectorMC[9] = 0.;   // imppar1, dummy value, meaningless in MC, in micron
  vectorMC[10] = 0.;   // imppar2, dummy value, meaningless in MC, in micron
  vectorMC[11] = fzMCVertex;    // z of reconstructed of primary vertex
  vectorMC[12] = fCentValue; // reconstructed centrality value 
  vectorMC[13] = 1.;           // always filling with 1 at MC level 

	if (fDecay==kLctopKpi){
		vectorMC[14] = 0.;
		vectorMC[15] = 0.;
		vectorMC[16] = 0.;
		vectorMC[17] = 0.;
	}
	

  bGenValues = kTRUE;
  return bGenValues;
}


//____________________________________________
Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
{ 
  // Fill vector (see above) with reconstructed quantities
  Bool_t bFillRecoValues=kFALSE;
  
  Int_t pdgCand = -1;
  if(fDecay==kDplustoKpipi){
    pdgCand=411;
  }else if(fDecay==kDstoKKpi){
    pdgCand=431;
  }else if(fDecay==kLctopKpi){
	pdgCand=4122;
   // AliError("LambdaC not yet implemented");
   // return bFillRecoValues;
  }else{
    AliError("WRONG DECAY SETTING");
    return bFillRecoValues;
  }

  AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
  Short_t charge=decay3->Charge();
  Double_t rapidity=decay3->Y(pdgCand);
  Double_t cT=decay3->Ct(pdgCand); 
  Double_t pt = decay3->Pt();
  Double_t cosPointingAngle = decay3->CosPointingAngle();
  Double_t phi = decay3->Phi();
  Double_t dist12= decay3->GetDist12toPrim();
  Double_t dist23 = decay3->GetDist23toPrim();
  Double_t sigmVert = decay3->GetSigmaVert();
	

  Int_t daughtSorted[3];
  Int_t tmpIndex=0;
  Int_t nDauLS=0;
  Int_t nDauOS=0;
  for(Int_t iDau=0; iDau<3; iDau++){
    AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
    Int_t label = TMath::Abs(trk->GetLabel());
    Short_t chargedau=trk->Charge();
    if(chargedau==charge){
      daughtSorted[tmpIndex]=label;
      tmpIndex=2;
      nDauLS++;
    }else{
      daughtSorted[1]=label;
      nDauOS++;
    }
  }

  if(nDauLS!=2 || nDauOS!=1){
    AliError("Wrong decay channel: number of OS and LS tracks not OK");
    return bFillRecoValues;
  }
  
  if(daughtSorted[0]>daughtSorted[2]){
    Int_t tmp=daughtSorted[0];
    daughtSorted[0]=daughtSorted[2];
    daughtSorted[2]=tmp;
  }
  
	Double_t d0prong0 = decay3->Getd0Prong(daughtSorted[0]);
	Double_t d0prong1 = decay3->Getd0Prong(daughtSorted[1]);
	Double_t d0prong2 = decay3->Getd0Prong(daughtSorted[2]);
	
  vectorReco[0] = pt;
  vectorReco[1] = rapidity;
  vectorReco[2] = phi;
  vectorReco[3] = cT*1.E4;  // in micron
  vectorReco[4] = cosPointingAngle;  // in micron
  vectorReco[5] = decay3->PtProng(daughtSorted[0]);
  vectorReco[6] = decay3->PtProng(daughtSorted[1]);
  vectorReco[7] = decay3->PtProng(daughtSorted[2]);
  vectorReco[8] = decay3->Getd0Prong(daughtSorted[0]);
  vectorReco[9] = decay3->Getd0Prong(daughtSorted[1]);
  vectorReco[10] = decay3->Getd0Prong(daughtSorted[2]);
  vectorReco[11] = fzPrimVertex;    // z of reconstructed of primary vertex
  vectorReco[12] = fCentValue; //reconstructed centrality value
  vectorReco[13] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 

	if(fDecay==kLctopKpi){  
		Double_t sumd02 =(d0prong0*d0prong0 + d0prong1*d0prong1 + d0prong2*d0prong2); 
		vectorReco[14] = dist12*1.E4;
		vectorReco[15] = dist23*1.E4;
		vectorReco[16] = sigmVert*1.E4;
		vectorReco[17] = sumd02*1.E8;
	}
	

  bFillRecoValues = kTRUE;
  return bFillRecoValues;
}


//_____________________________________________________________
Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
{ 
  // Check the pdg codes of the daughters
  Bool_t checkCD = kFALSE;

  Int_t pdgCand = -1;
  Int_t pdgDaughter[3]={-1,-1,-1};
  if(fDecay==kDplustoKpipi){
    pdgCand=411;
    pdgDaughter[0]=321;
    pdgDaughter[1]=211;
    pdgDaughter[2]=211;
  }else if(fDecay==kDstoKKpi){
    pdgCand=431;
    pdgDaughter[0]=321;
    pdgDaughter[1]=321;
    pdgDaughter[2]=211;
  }else if(fDecay==kLctopKpi){
	  pdgCand=4122;
	  pdgDaughter[0]=2212;
	  pdgDaughter[1]=321;
	  pdgDaughter[2]=211;
	  
  //  AliError("LambdaC not yet implemented");
  //  return checkCD;
  }else{
    AliError("WRONG DECAY SETTING");
    return checkCD;    
  }


  Int_t daughter[3];

  Int_t nDau=fmcPartCandidate->GetNDaughters();
  Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0); 
  if(nDau==3){
    for(Int_t iDau=0; iDau<3; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return checkCD;
      }
      daughter[iDau]=TMath::Abs(part->GetPdgCode());
    }
  }else if(nDau==2){
    Int_t nDauFound=0;
    for(Int_t iDau=0; iDau<2; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return checkCD;
      }
      Int_t pdgCode=TMath::Abs(part->GetPdgCode());
      if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
	if(nDauFound>=3) return checkCD;
	daughter[nDauFound]=pdgCode;
	nDauFound++;
      }else{
	Int_t nDauRes=part->GetNDaughters();
	if(nDauRes!=2) return checkCD;
	Int_t labelFirstDauRes = part->GetDaughter(0); 	
	for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	  Int_t indDR = labelFirstDauRes+iDauRes;
	  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
	  if(!partDR){
	    AliError("Daughter particle not found in MC array");
	    return checkCD;
	  }
	  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
	  if(nDauFound>=3) return checkCD;
	  daughter[nDauFound]=pdgCodeDR;
	  nDauFound++;
	}
      }
    }
  }else{
    return checkCD;
  }
  for(Int_t iDau1=0; iDau1<3; iDau1++){
    for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
      if(daughter[iDau1]<daughter[iDau2]){
	Int_t tmp=daughter[iDau1];
	daughter[iDau1]=daughter[iDau2];
	daughter[iDau2]=tmp;
      }
    }
  }
  for(Int_t iDau=0; iDau<3; iDau++){
    if(daughter[iDau]!=pdgDaughter[iDau]){
      AliDebug(2, "Wrong decay channel from MC, skipping!!");
      return checkCD;  
    }
  }
  
  checkCD = kTRUE;
  return checkCD;
  
}
Commit	Line	Data
043062fe	1	/**************************************************************************
	2	* Copyright(c) 2007-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
7b45817b	16	/* $Id$ */
043062fe	17
	18	///////////////////////////////////////////////////////////////////
	19	// //
	20	// Class to compute variables for correction framework //
	21	// for 3-body decays of D mesons (D+, Ds, Lc) //
	22	// in bins of cut variables //
	23	// Origin: Francesco Prino (prino@to.infn.it) //
	24	// Renu Bala (bala@to.infn.it) //
1c44b663	25	// Davide Caffarri (cafarri@pd.infn.it) //
043062fe	26	///////////////////////////////////////////////////////////////////
	27
	28	#include "AliAODMCParticle.h"
	29	#include "AliAODEvent.h"
	30	#include "TClonesArray.h"
	31	#include "AliCFVertexingHF.h"
	32	#include "AliESDtrack.h"
	33	#include "TDatabasePDG.h"
	34
	35	#include "AliCFVertexingHF3Prong.h"
	36	#include "AliCFContainer.h"
	37
	38	ClassImp(AliCFVertexingHF3Prong)
	39
	40
	41	//_________________________________________
	42	AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay):
	43	AliCFVertexingHF(),
	44	fDecay(decay)
	45	{
	46	//
	47	SetNProngs(3);
2bf2e62b	48
	49	fPtAccCut=new Float_t[fProngs];
	50	fEtaAccCut=new Float_t[fProngs];
	51	for(Int_t iP=0; iP<fProngs; iP++){
	52	fPtAccCut[iP]=0.1;
	53	fEtaAccCut[iP]=0.9;
	54	}
	55
043062fe	56	}
	57	//_________________________________________
	58	AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
	59	AliCFVertexingHF(mcArray, originDselection),
	60	fDecay(decay)
	61
	62	{
	63	//
	64	SetNProngs(3);
2bf2e62b	65	fPtAccCut=new Float_t[fProngs];
	66	fEtaAccCut=new Float_t[fProngs];
	67	for(Int_t iP=0; iP<fProngs; iP++){
	68	fPtAccCut[iP]=0.1;
	69	fEtaAccCut[iP]=0.9;
	70	}
043062fe	71	}
	72
	73
	74	//_____________________________________
	75	AliCFVertexingHF3Prong& AliCFVertexingHF3Prong::operator=(const AliCFVertexingHF3Prong& c){
	76	//
	77	if (this != &c) {
	78
	79	AliCFVertexingHF::operator=(c);
	80
	81	}
	82	return *this;
	83	}
	84
	85	//__________________________________________
	86	Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
	87	// Checks if candidate is signal and D meson is present in MC array
	88
	89	Bool_t bSignAssoc = kFALSE;
	90	fRecoCandidate = recoCand;
	91
	92	if (!fRecoCandidate) {
	93	AliError("fRecoCandidate not found, problem in assignement\n");
	94	return bSignAssoc;
	95	}
	96
	97	Int_t pdgCand = -1;
	98	Int_t pdgDaughter[3]={-1,-1,-1};
	99	if(fDecay==kDplustoKpipi){
	100	pdgCand=411;
	101	pdgDaughter[0]=321;
	102	pdgDaughter[1]=211;
	103	pdgDaughter[2]=211;
	104	}else if(fDecay==kDstoKKpi){
	105	pdgCand=431;
	106	pdgDaughter[0]=321;
	107	pdgDaughter[1]=321;
	108	pdgDaughter[2]=211;
	109	}else if(fDecay==kLctopKpi){
1c44b663	110	pdgCand=4122;
	111	pdgDaughter[0]=2212;
	112	pdgDaughter[1]=321;
	113	pdgDaughter[2]=211;
043062fe	114	}else{
	115	AliError("WRONG DECAY SETTING");
	116	return bSignAssoc;
	117	}
	118
	119	Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);
	120	if (mcLabel == -1) return bSignAssoc;
fbec9fa9	121
1f780958	122	if (fRecoCandidate->NumberOfFakeDaughters()>0){
	123	fFake = 0; // fake candidate
	124	if (fFakeSelection==1) return bSignAssoc;
	125	}
	126	if (fRecoCandidate->NumberOfFakeDaughters()==0){
	127	fFake = 2; // non-fake candidate
	128	if (fFakeSelection==2) return bSignAssoc;
	129	}
fbec9fa9	130
043062fe	131	SetMCLabel(mcLabel);
	132	fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
	133
	134	if (!fmcPartCandidate){
	135	AliDebug(3,"No part candidate");
	136	return bSignAssoc;
	137	}
	138
	139	bSignAssoc = kTRUE;
	140	return bSignAssoc;
	141	}
	142
	143	//______________________________________________
	144	Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
	145	//
	146	// collecting all the necessary info from MC particle and fill vectorMC: 12 variables
	147	// pt_D
	148	// y_D
	149	// phi_D
	150	// ctau
	151	// cos point
	152	// pt_1
	153	// pt_2
	154	// pt_3
	155	// d0_1
	156	// d0_2
	157	// d0_3
	158	// zPrimVert
b7af2639	159	// centrality
043062fe	160
	161	Bool_t bGenValues = kFALSE;
	162
	163	Int_t pdgCand = -1;
	164	if(fDecay==kDplustoKpipi){
	165	pdgCand=411;
	166	}else if(fDecay==kDstoKKpi){
	167	pdgCand=431;
	168	}else if(fDecay==kLctopKpi){
1c44b663	169	pdgCand=4122;
043062fe	170	}else{
	171	AliError("WRONG DECAY SETTING");
	172	return bGenValues;
	173	}
	174
	175	Double_t vertD[3] = {0,0,0}; // D origin
	176	fmcPartCandidate->XvYvZv(vertD); // cm
	177
	178	Int_t nprongs = 3;
	179	Int_t daughter[3];
	180	Short_t charge = fmcPartCandidate->Charge();
	181
	182	// order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
	183	// the 2 LS are ordered so that in pos. 0 there is the one with lower label value
	184	Int_t index=0;
	185	Int_t nDauLS=0;
	186	Int_t nDauOS=0;
	187
	188
	189	Int_t nDau=fmcPartCandidate->GetNDaughters();
	190	Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
	191	if(nDau==3){
	192	for(Int_t iDau=0; iDau<3; iDau++){
	193	Int_t ind = labelFirstDau+iDau;
	194	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
968b84b9	195	if(!part){
	196	AliError("Daughter particle not found in MC array");
	197	return bGenValues;
	198	}
043062fe	199	Short_t signDau=part->Charge();
	200	if(signDau==charge){
	201	nDauLS++;
	202	daughter[index] = ind;
	203	index=2;
	204	}else{
	205	daughter[1] = ind;
	206	nDauOS++;
	207	}
	208	}
	209	}else if(nDau==2){
	210	for(Int_t iDau=0; iDau<2; iDau++){
	211	Int_t ind = labelFirstDau+iDau;
	212	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
968b84b9	213	if(!part){
	214	AliError("Daughter particle not found in MC array");
	215	return bGenValues;
	216	}
043062fe	217	Int_t pdgCode=TMath::Abs(part->GetPdgCode());
	218	if(pdgCode==211 \|\| pdgCode==321 \|\| pdgCode==2212){
	219	Short_t signDau=part->Charge();
	220	if(signDau==charge){
	221	nDauLS++;
	222	daughter[index] = ind;
	223	index=2;
	224	}else{
	225	daughter[1] = ind;
	226	nDauOS++;
	227	}
	228	}else{
	229	Int_t nDauRes=part->GetNDaughters();
	230	if(nDauRes!=2){
	231	AliError("Wrong resonant decay");
	232	return bGenValues;
	233	}
	234	Int_t labelFirstDauRes = part->GetDaughter(0);
	235	for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	236	Int_t indDR = labelFirstDauRes+iDauRes;
	237	AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
968b84b9	238	if(!partDR){
	239	AliError("Daughter particle not found in MC array");
	240	return bGenValues;
	241	}
043062fe	242	Short_t signDau=partDR->Charge();
	243	if(signDau==charge){
	244	nDauLS++;
	245	daughter[index] = ind;
	246	index=2;
	247	}else{
	248	daughter[1] = ind;
	249	nDauOS++;
	250	}
	251	}
	252	}
	253	}
	254	}else{
	255	AliError(Form("Wrong number of daughters %d",nDau));
	256	return bGenValues;
	257	}
	258
	259	if(nDauLS!=2 \|\| nDauOS!=1){
	260	AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
	261	return bGenValues;
	262	}
	263	if(daughter[0]>daughter[2]){
	264	Int_t tmp=daughter[0];
	265	daughter[0]=daughter[2];
	266	daughter[2]=tmp;
	267	}
	268
	269	// getting the momentum from the daughters and decay vertex
	270	Double_t px[3],py[3],pz[3],pt[3];
	271	Double_t vertDec[3] = {0,0,0}; // decay vertex
	272	for(Int_t iDau=0; iDau<3; iDau++){
	273	AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
968b84b9	274	if(!part){
	275	AliError("Daughter particle not found in MC array");
	276	return bGenValues;
	277	}
043062fe	278	px[iDau]=part->Px();
	279	py[iDau]=part->Py();
	280	pz[iDau]=part->Pz();
	281	pt[iDau]=part->Pt();
	282	if(iDau==0) part->XvYvZv(vertDec);
	283	}
	284
	285	Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
	286
	287	AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
	288	Double_t cT = decay->Ct(pdgCand);
	289
	290	vectorMC[0] = fmcPartCandidate->Pt();
	291	vectorMC[1] = fmcPartCandidate->Y() ;
	292	vectorMC[2] = fmcPartCandidate->Phi();
	293	vectorMC[3] = cT*1.E4 ; // in micron
	294	vectorMC[4] = 1.01; // cos pointing angle, dummy value, meaningless in MC
	295	vectorMC[5] = pt[0];
	296	vectorMC[6] = pt[1];
	297	vectorMC[7] = pt[2];
	298	vectorMC[8] = 0.; // imppar0, dummy value, meaningless in MC
	299	vectorMC[9] = 0.; // imppar1, dummy value, meaningless in MC, in micron
	300	vectorMC[10] = 0.; // imppar2, dummy value, meaningless in MC, in micron
	301	vectorMC[11] = fzMCVertex; // z of reconstructed of primary vertex
b7af2639	302	vectorMC[12] = fCentValue; // reconstructed centrality value
1f780958	303	vectorMC[13] = 1.; // always filling with 1 at MC level
b7af2639	304
1c44b663	305	if (fDecay==kLctopKpi){
	306	vectorMC[14] = 0.;
	307	vectorMC[15] = 0.;
	308	vectorMC[16] = 0.;
	309	vectorMC[17] = 0.;
	310	}
	311
043062fe	312
	313	bGenValues = kTRUE;
	314	return bGenValues;
	315	}
	316
	317
	318	//____________________________________________
	319	Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
	320	{
	321	// Fill vector (see above) with reconstructed quantities
	322	Bool_t bFillRecoValues=kFALSE;
	323
	324	Int_t pdgCand = -1;
	325	if(fDecay==kDplustoKpipi){
	326	pdgCand=411;
	327	}else if(fDecay==kDstoKKpi){
	328	pdgCand=431;
	329	}else if(fDecay==kLctopKpi){
1c44b663	330	pdgCand=4122;
	331	// AliError("LambdaC not yet implemented");
	332	// return bFillRecoValues;
043062fe	333	}else{
	334	AliError("WRONG DECAY SETTING");
	335	return bFillRecoValues;
	336	}
	337
	338	AliAODRecoDecayHF3Prong decay3 = (AliAODRecoDecayHF3Prong)fRecoCandidate;
	339	Short_t charge=decay3->Charge();
	340	Double_t rapidity=decay3->Y(pdgCand);
	341	Double_t cT=decay3->Ct(pdgCand);
	342	Double_t pt = decay3->Pt();
	343	Double_t cosPointingAngle = decay3->CosPointingAngle();
	344	Double_t phi = decay3->Phi();
1c44b663	345	Double_t dist12= decay3->GetDist12toPrim();
	346	Double_t dist23 = decay3->GetDist23toPrim();
	347	Double_t sigmVert = decay3->GetSigmaVert();
	348
043062fe	349
	350	Int_t daughtSorted[3];
	351	Int_t tmpIndex=0;
	352	Int_t nDauLS=0;
	353	Int_t nDauOS=0;
	354	for(Int_t iDau=0; iDau<3; iDau++){
	355	AliAODTrack trk = (AliAODTrack)decay3->GetDaughter(iDau);
b7af2639	356	Int_t label = TMath::Abs(trk->GetLabel());
043062fe	357	Short_t chargedau=trk->Charge();
	358	if(chargedau==charge){
	359	daughtSorted[tmpIndex]=label;
	360	tmpIndex=2;
	361	nDauLS++;
	362	}else{
	363	daughtSorted[1]=label;
	364	nDauOS++;
	365	}
	366	}
	367
	368	if(nDauLS!=2 \|\| nDauOS!=1){
	369	AliError("Wrong decay channel: number of OS and LS tracks not OK");
	370	return bFillRecoValues;
	371	}
	372
	373	if(daughtSorted[0]>daughtSorted[2]){
	374	Int_t tmp=daughtSorted[0];
	375	daughtSorted[0]=daughtSorted[2];
	376	daughtSorted[2]=tmp;
	377	}
	378
1c44b663	379	Double_t d0prong0 = decay3->Getd0Prong(daughtSorted[0]);
	380	Double_t d0prong1 = decay3->Getd0Prong(daughtSorted[1]);
	381	Double_t d0prong2 = decay3->Getd0Prong(daughtSorted[2]);
	382
043062fe	383	vectorReco[0] = pt;
	384	vectorReco[1] = rapidity;
	385	vectorReco[2] = phi;
	386	vectorReco[3] = cT*1.E4; // in micron
	387	vectorReco[4] = cosPointingAngle; // in micron
	388	vectorReco[5] = decay3->PtProng(daughtSorted[0]);
	389	vectorReco[6] = decay3->PtProng(daughtSorted[1]);
	390	vectorReco[7] = decay3->PtProng(daughtSorted[2]);
	391	vectorReco[8] = decay3->Getd0Prong(daughtSorted[0]);
	392	vectorReco[9] = decay3->Getd0Prong(daughtSorted[1]);
	393	vectorReco[10] = decay3->Getd0Prong(daughtSorted[2]);
	394	vectorReco[11] = fzPrimVertex; // z of reconstructed of primary vertex
b7af2639	395	vectorReco[12] = fCentValue; //reconstructed centrality value
1f780958	396	vectorReco[13] = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
b7af2639	397
1c44b663	398	if(fDecay==kLctopKpi){
	399	Double_t sumd02 =(d0prong0d0prong0 + d0prong1d0prong1 + d0prong2*d0prong2);
	400	vectorReco[14] = dist12*1.E4;
	401	vectorReco[15] = dist23*1.E4;
	402	vectorReco[16] = sigmVert*1.E4;
	403	vectorReco[17] = sumd02*1.E8;
	404	}
	405
b7af2639	406
043062fe	407	bFillRecoValues = kTRUE;
	408	return bFillRecoValues;
	409	}
	410
	411
	412	//_____________________________________________________________
	413	Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
	414	{
	415	// Check the pdg codes of the daughters
	416	Bool_t checkCD = kFALSE;
	417
	418	Int_t pdgCand = -1;
	419	Int_t pdgDaughter[3]={-1,-1,-1};
	420	if(fDecay==kDplustoKpipi){
	421	pdgCand=411;
	422	pdgDaughter[0]=321;
	423	pdgDaughter[1]=211;
	424	pdgDaughter[2]=211;
	425	}else if(fDecay==kDstoKKpi){
	426	pdgCand=431;
	427	pdgDaughter[0]=321;
	428	pdgDaughter[1]=321;
	429	pdgDaughter[2]=211;
	430	}else if(fDecay==kLctopKpi){
1c44b663	431	pdgCand=4122;
	432	pdgDaughter[0]=2212;
	433	pdgDaughter[1]=321;
	434	pdgDaughter[2]=211;
	435
	436	// AliError("LambdaC not yet implemented");
	437	// return checkCD;
043062fe	438	}else{
	439	AliError("WRONG DECAY SETTING");
	440	return checkCD;
	441	}
	442
	443
	444	Int_t daughter[3];
	445
	446	Int_t nDau=fmcPartCandidate->GetNDaughters();
	447	Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
	448	if(nDau==3){
	449	for(Int_t iDau=0; iDau<3; iDau++){
	450	Int_t ind = labelFirstDau+iDau;
	451	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
968b84b9	452	if(!part){
	453	AliError("Daughter particle not found in MC array");
	454	return checkCD;
	455	}
043062fe	456	daughter[iDau]=TMath::Abs(part->GetPdgCode());
	457	}
	458	}else if(nDau==2){
	459	Int_t nDauFound=0;
	460	for(Int_t iDau=0; iDau<2; iDau++){
	461	Int_t ind = labelFirstDau+iDau;
	462	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
968b84b9	463	if(!part){
	464	AliError("Daughter particle not found in MC array");
	465	return checkCD;
	466	}
043062fe	467	Int_t pdgCode=TMath::Abs(part->GetPdgCode());
	468	if(pdgCode==211 \|\| pdgCode==321 \|\| pdgCode==2212){
	469	if(nDauFound>=3) return checkCD;
	470	daughter[nDauFound]=pdgCode;
	471	nDauFound++;
	472	}else{
	473	Int_t nDauRes=part->GetNDaughters();
	474	if(nDauRes!=2) return checkCD;
	475	Int_t labelFirstDauRes = part->GetDaughter(0);
	476	for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	477	Int_t indDR = labelFirstDauRes+iDauRes;
	478	AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
968b84b9	479	if(!partDR){
	480	AliError("Daughter particle not found in MC array");
	481	return checkCD;
	482	}
043062fe	483	Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
	484	if(nDauFound>=3) return checkCD;
	485	daughter[nDauFound]=pdgCodeDR;
	486	nDauFound++;
	487	}
	488	}
	489	}
	490	}else{
	491	return checkCD;
	492	}
	493	for(Int_t iDau1=0; iDau1<3; iDau1++){
	494	for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
	495	if(daughter[iDau1]<daughter[iDau2]){
	496	Int_t tmp=daughter[iDau1];
	497	daughter[iDau1]=daughter[iDau2];
	498	daughter[iDau2]=tmp;
	499	}
	500	}
	501	}
	502	for(Int_t iDau=0; iDau<3; iDau++){
	503	if(daughter[iDau]!=pdgDaughter[iDau]){
	504	AliDebug(2, "Wrong decay channel from MC, skipping!!");
	505	return checkCD;
	506	}
	507	}
	508
	509	checkCD = kTRUE;
	510	return checkCD;
	511
	512	}