[u/mrichter/AliRoot.git] / PWGHF / 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),
  fGenDsOption(kCountResonant)
 {
  // 
  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),
  fGenDsOption(kCountResonant)
{
  //
  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;
  }

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