[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"
#include "AliCFTaskVertexingHF.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);
	
	switch (fConfiguration){
	case AliCFTaskVertexingHF::kSnail:
		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] = fzMCVertex;    // z of reconstructed of primary vertex
		vectorMC[9] = fCentValue; // reconstructed centrality value 
		vectorMC[10] = 1.;           // fake: always filling with 1 at MC level 
		vectorMC[11] = 1.01; // dummy value for cosPointingXY  multiplicity
		vectorMC[12] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
		vectorMC[13] = fMultiplicity; // reconstructed multiplicity
		
		if (fDecay==kLctopKpi){
			vectorMC[11] = 0.; //dist12
			vectorMC[12] = 0.; //dist23
			vectorMC[13] = 0.; //sigmaVtx
			vectorMC[14] = 0.; //sumd02
			vectorMC[15] = 1.01; // dummy value for cosPointingXY  multiplicity
			vectorMC[16] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
			vectorMC[17] = 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;
	}
	
	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();
	Double_t cosPointingAngleXY = decay3->CosPointingAngleXY();
	Double_t normDecayLengthXY = decay3->NormalizedDecayLengthXY();
		
	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]);
	
	switch (fConfiguration){
	case AliCFTaskVertexingHF::kSnail:
		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] = fzPrimVertex;    // z of reconstructed of primary vertex
		vectorReco[9] = fCentValue; //reconstructed centrality value
		vectorReco[10] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
		vectorReco[11] = cosPointingAngleXY; 
		vectorReco[12] = normDecayLengthXY; // in cm
		vectorReco[13] = fMultiplicity; // reconstructed multiplicity
		
		if(fDecay==kLctopKpi){  
			Double_t sumd02 =(d0prong0*d0prong0 + d0prong1*d0prong1 + d0prong2*d0prong2); 
			vectorReco[11] = dist12*1.E4;
			vectorReco[12] = dist23*1.E4;
			vectorReco[13] = sigmVert*1.E4;
			vectorReco[14] = sumd02*1.E8;
			vectorReco[15] = cosPointingAngleXY; 
			vectorReco[16] = normDecayLengthXY; // in cm
			vectorReco[17] = 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;
}


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