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

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

//-----------------------------------------------------------------------
// Class for HF corrections as a function of many variables and step 
// Author : C. Zampolli, CERN
// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
// Base class for HF Unfolding - agrelli@uu.nl
//-----------------------------------------------------------------------

#include "TParticle.h"
#include "TClonesArray.h"
#include "AliAODMCParticle.h"
#include "AliAODRecoDecayHF.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoDecayHF3Prong.h"
#include "AliAODRecoDecayHF4Prong.h"
#include "AliAODMCHeader.h"
#include "AliAODEvent.h"
#include "AliLog.h"
#include "AliESDtrackCuts.h"
#include "AliESDtrack.h"

#include "AliCFVertexingHF.h"

//___________________________________________________________
AliCFVertexingHF::AliCFVertexingHF() :
	fmcArray(0x0),
	fRecoCandidate(0),
	fmcPartCandidate(0x0),
	fNDaughters(0),
	fNVar(0),
	fzPrimVertex(0),
	fzMCVertex(0),
	fFillFromGenerated(0),
	fOriginDselection(0),
	fKeepDfromB(kFALSE),
	fKeepDfromBOnly(kFALSE),
	fmcLabel(0),
	fProngs(-1)
{
	return;
}


//_____________________________________________________
AliCFVertexingHF::AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection) :
	fmcArray(mcArray),
	fRecoCandidate(0),
	fmcPartCandidate(0x0),
	fNDaughters(0),
	fNVar(0),
	fzPrimVertex(0),
	fzMCVertex(0),
	fFillFromGenerated(0),
	fOriginDselection(0),
	fKeepDfromB(kFALSE),
	fKeepDfromBOnly(kFALSE),
	fmcLabel(0),
	fProngs(-1)
{

  
  SetDselection(originDselection);


  return;
}


//_______________________________________________________
AliCFVertexingHF::~AliCFVertexingHF()
{
	if (fmcArray) delete fmcArray;
	if (fRecoCandidate) delete fRecoCandidate;
	if (fmcPartCandidate) delete fmcPartCandidate;
}


//_____________________________________________________
AliCFVertexingHF& AliCFVertexingHF::operator=(const AliCFVertexingHF& c){
	
	if (this!= &c){
		TObject::operator=(c);
		fmcArray = c.fmcArray;
		fRecoCandidate = c.fRecoCandidate;
		fmcPartCandidate = c.fmcPartCandidate;
		fNDaughters = c.fNDaughters;
		fNVar = c.fNVar;
		fzPrimVertex = c.fzPrimVertex;
		fzMCVertex = c.fzMCVertex;
		fFillFromGenerated = c.fFillFromGenerated;
		fOriginDselection = c.fOriginDselection;
		fKeepDfromB = c.fKeepDfromB;
		fKeepDfromBOnly = c.fKeepDfromBOnly;
		fmcLabel = c.fmcLabel;
		
	}
	
	return *this;
}

//____________________________________________________
AliCFVertexingHF::AliCFVertexingHF(const AliCFVertexingHF &c) :
        TObject(c),
	fmcArray(c.fmcArray),
	fRecoCandidate(c.fRecoCandidate),
	fmcPartCandidate(c.fmcPartCandidate),
	fNDaughters(c.fNDaughters),
	fNVar(c.fNVar),
	fzPrimVertex(c.fzPrimVertex),
	fzMCVertex(c.fzMCVertex),
	fFillFromGenerated(c.fFillFromGenerated),
	fOriginDselection (c.fOriginDselection),
	fKeepDfromB (c.fKeepDfromB),
	fKeepDfromBOnly (c.fKeepDfromBOnly),
	fmcLabel(c.fmcLabel),
	fProngs(c.fProngs)
{
   

}

//___________________________________________________________
void AliCFVertexingHF::SetDselection(UShort_t originDselection){

 fOriginDselection = originDselection;
 
 if (fOriginDselection == 0) {
   fKeepDfromB = kFALSE;
   fKeepDfromBOnly = kFALSE;
 }
 
 if (fOriginDselection == 1) {
   fKeepDfromB = kTRUE;
   fKeepDfromBOnly = kTRUE;
 }
 
 if (fOriginDselection == 2) {
   fKeepDfromB = kTRUE;
   fKeepDfromBOnly = kFALSE;
 }
 
 return;

}

//______________________________________________________
void AliCFVertexingHF::SetMCCandidateParam(Int_t label){
	
  fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
  fNDaughters = fmcPartCandidate->GetNDaughters();
  return;
}


//____________________________________________________________
Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const{
	
 
  Int_t cquarks = 0;
  if (mcPart->GetPdgCode() == 4) cquarks++; 
  if (mcPart->GetPdgCode() == -4) cquarks++; 
  if (!mcPart) {
    AliWarning("Particle not found in tree, skipping\n"); 
    return cquarks;
  } 
  
  return cquarks;
}


//________________________________________________________
Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const {

  Int_t pdgGranma = CheckOrigin();
  if (pdgGranma == -9999){
    printf ("This particle come from a B decay channel but the we keep only the prompt charm particles\n");	
    return kFALSE;
  }	
  
 if (pdgGranma == -999){
    printf ("This particle come from a prompt charm particles but we want only the ones coming from B\n");	
    return kFALSE;
  }	

  if (!CheckMCDaughters()) return kFALSE;
  if (!CheckMCChannelDecay()) return kFALSE;
  return kTRUE;
}

//_________________________________________________________________________________________________
Int_t AliCFVertexingHF::CheckOrigin() const {		
  //
  // checking whether the mother of the particles come from a charm or a bottom quark
  //
  
  Int_t pdgGranma = 0;
  Int_t mother = 0;
  mother = fmcPartCandidate->GetMother();
  Int_t istep = 0;
  while (mother >0 ){
    istep++;
    AliDebug(2,Form("mother at step %d = %d", istep, mother));
    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
    pdgGranma = mcGranma->GetPdgCode();
    AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
    Int_t abspdgGranma = TMath::Abs(pdgGranma);
    if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
      if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
      
      else{
	return pdgGranma;
      }
    }
    else {
    if (fKeepDfromBOnly) return -999;
    }
    
    mother = mcGranma->GetMother();
  }
  
  return pdgGranma;
}


//___________________________________________
Bool_t AliCFVertexingHF::CheckMCDaughters()const {
	
  AliAODMCParticle *mcPartDaughter;
  Bool_t checkDaughters = kFALSE;

  Int_t label0 = fmcPartCandidate->GetDaughter(0);
  Int_t label1 = fmcPartCandidate->GetDaughter(1);
  if (label1==0 || label0 == 0){
    AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
	return checkDaughters;  
  }

  if (label1 - label0 != fProngs-1){
    AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
    return checkDaughters;  
  }

  for (Int_t iProng = 0; iProng<fProngs; iProng++){
    mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));    
    if (!mcPartDaughter) {
      AliWarning("At least one Daughter Particle not found in tree, skipping"); 
      return checkDaughters;  
    }
  }

  checkDaughters = kTRUE;
  return checkDaughters;
}

	
//______________________________________________________
Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
{
// fill the container for Generator level selection
  Bool_t mcContainerFilled = kFALSE;  

  Double_t* vectorMC = new Double_t[fNVar];
  for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
  
  if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
    for (Int_t iVar = 0; iVar<fNVar; iVar++){
      
      containerInputMC[iVar] = vectorMC[iVar];
    }
    
    mcContainerFilled = kTRUE;		
  }
  delete vectorMC;
  return mcContainerFilled;	
}

//______________________________________________________
Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput) {  

// fill the container for Reconstrucred level selection

  Bool_t recoContainerFilled = kFALSE;
  Double_t* vectorValues = new Double_t[fNVar];
  Double_t* vectorReco = new Double_t[fNVar];  
  
  for (Int_t iVar = 0; iVar<fNVar; iVar++) {
    vectorValues[iVar]= 9999.;
    vectorReco[iVar]=9999.;
  }

  if(fFillFromGenerated){
    //filled with MC values
    if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
      for (Int_t iVar = 0; iVar<fNVar; iVar++){
		  containerInput[iVar] = vectorValues[iVar];
		}
		recoContainerFilled = kTRUE;		
		}
	}
	else{
    //filled with Reco values
		
		if (GetRecoValuesFromCandidate(&vectorReco[0])){
			for (Int_t iVar = 0; iVar<fNVar; iVar++){
				containerInput[iVar] = vectorReco[iVar];
			}
			recoContainerFilled = kTRUE;		
		}
	}
  
  delete vectorValues;
  delete vectorReco;
  return recoContainerFilled;	
}

//_____________________________________________________
Bool_t AliCFVertexingHF::MCAcceptanceStep() const
{

  Bool_t bMCAccStep = kFALSE;
  
  AliAODMCParticle *mcPartDaughter;
  Int_t label0 = fmcPartCandidate->GetDaughter(0);
  Int_t label1 = fmcPartCandidate->GetDaughter(1);
  if (label1==0 || label0 == 0){
    AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
    return bMCAccStep;  
  }

  if (label1 - label0 != fProngs-1){
    AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
	return bMCAccStep;  
  }

  for (Int_t iProng = 0; iProng<fProngs; iProng++){
    mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));    
    if (!mcPartDaughter) {
      AliWarning("At least one Daughter Particle not found in tree, skipping"); 
      return bMCAccStep;  
    }
    Double_t eta = mcPartDaughter->Eta();
    Double_t pt = mcPartDaughter->Pt();
	  
	  //set values of eta and pt in the constructor.
    if (TMath::Abs(eta) > 0.9 || pt < 0.1){
		AliDebug(3,"At least one daughter has eta>0.9 or pt < 0.1 \n"); 
      return bMCAccStep;
    }

  }

  
  bMCAccStep = kTRUE;
  return bMCAccStep; 
  
}

//_____________________________________________________
Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts *trackCuts) const
{		
  // check on the kTPCrefit and kITSrefit conditions of the daughters
	Bool_t bRefitStep = kFALSE;

	Int_t label0 = fmcPartCandidate->GetDaughter(0);
	Int_t label1 = fmcPartCandidate->GetDaughter(1);

	if (label1==0 || label0 == 0){
		AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
        return bRefitStep;  
	}

	if (label1 - label0 != fProngs-1){
		AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
		//AliInfo(Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
		return bRefitStep;  
	}

	Int_t foundDaughters = 0;

	if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
    
		for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
			AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
			if (track->GetLabel()>= label0 && track->GetLabel()<=label1){
				foundDaughters++;
				printf("daughter %d \n",foundDaughters);
				if(trackCuts->GetRequireTPCRefit()){
					if((track->GetStatus()&AliESDtrack::kTPCrefit)) {
						bRefitStep = kTRUE;
					}
					else {
						AliDebug(3, "Refit cut not passed , missing TPC refit\n");
						AliInfo( "Refit cut not passed , missing TPC refit\n");
						return kFALSE;
					}
				}
				
				if (trackCuts->GetRequireITSRefit()) {
					if((track->GetStatus()&AliESDtrack::kITSrefit)){
						bRefitStep = kTRUE;
					}
					else {
						AliDebug(3, "Refit cut not passed , missing ITS refit\n");
						AliInfo("Refit cut not passed , missing ITS refit\n");
						return kFALSE;
					}
				}
			}      	
			if (foundDaughters == fProngs){
	 
				break;
			}
      
		}
    
	}
  
	return bRefitStep;
}

//____________________________________________________________________________

Bool_t AliCFVertexingHF::RecoStep() 
//check also vertex and ITS Refit and TPC Refit
{ 
  Bool_t bRecoStep = kFALSE;
  Int_t mcLabel = GetMCLabel();
  
  if (mcLabel == -1) {
    AliDebug(2,"No MC particle found");
    return bRecoStep;
  }
  else{
    fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
    if (!fmcPartCandidate){
      AliWarning("Could not find associated MC in AOD MC tree");
      return bRecoStep;
    }
    
  }
  
  Int_t pdgGranma = CheckOrigin();
  
  if (pdgGranma == -9999){
    printf ("This particle come from a B decay channel but we keep only prompt charm particles\n");
    return bRecoStep;
  }

  if (pdgGranma == -999){
    printf ("This particle come from a  prompt charm particle but we want only the ones coming from B\n");
    return bRecoStep;
  }
  
   
  bRecoStep=kTRUE;
  return bRecoStep;
  
  
}	
//____________________________________________
Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const {

  if (fRecoCandidate){
    Double_t etaProng = fRecoCandidate->EtaProng(iProng);  
    return etaProng;
  }
  return 999999;  
}
//______________________________________________________
Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const {

  if (fRecoCandidate){
    Double_t ptProng = fRecoCandidate->PtProng(iProng);  
    return ptProng;
  }
  return 999999;  
  
}

//____________________________________________________________________

Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts *trackCuts) const
{
  Bool_t bRecoAccStep = kFALSE;
  
  Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
  trackCuts->GetEtaRange(etaCutMin, etaCutMax);
  trackCuts->GetPtRange(ptCutMin, ptCutMax);
  
  Float_t etaProng=0., ptProng=0.; 
  
  for (Int_t iProng =0; iProng<fProngs; iProng++){
    
    etaProng = GetEtaProng(iProng);
    ptProng = GetPtProng(iProng);
    
    Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
    if (!acceptanceProng) {
      printf ("At least one reconstructed prong isn't in the acceptance\n");
      return bRecoAccStep;
    }
  }
  
  bRecoAccStep=kTRUE;
  return bRecoAccStep;
}
//___________________________________________________________

Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t *fill) const{

  fill = new Double_t[4];
  
  if(fmcPartCandidate){
    
    fill[0] = GetPtCand();
    fill[1] = GetYCand();
    
    fill[2] =  fmcPartCandidate->Pt(); 
    fill[3] =  fmcPartCandidate->Y(); 
    
    return kTRUE;
  }

  delete fill;
  return kFALSE;
}


//___________________________________________________________
Commit	Line	Data
379592af	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------------
	17	// Class for HF corrections as a function of many variables and step
	18	// Author : C. Zampolli, CERN
	19	// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
	20	// Base class for HF Unfolding - agrelli@uu.nl
	21	//-----------------------------------------------------------------------
	22
	23	#include "TParticle.h"
	24	#include "TClonesArray.h"
	25	#include "AliAODMCParticle.h"
	26	#include "AliAODRecoDecayHF.h"
	27	#include "AliAODRecoDecayHF2Prong.h"
	28	#include "AliAODRecoDecayHF3Prong.h"
	29	#include "AliAODRecoDecayHF4Prong.h"
	30	#include "AliAODMCHeader.h"
	31	#include "AliAODEvent.h"
	32	#include "AliLog.h"
	33	#include "AliESDtrackCuts.h"
	34	#include "AliESDtrack.h"
	35
	36	#include "AliCFVertexingHF.h"
	37
	38	//___________________________________________________________
	39	AliCFVertexingHF::AliCFVertexingHF() :
	40	fmcArray(0x0),
	41	fRecoCandidate(0),
	42	fmcPartCandidate(0x0),
	43	fNDaughters(0),
	44	fNVar(0),
	45	fzPrimVertex(0),
	46	fzMCVertex(0),
	47	fFillFromGenerated(0),
	48	fOriginDselection(0),
	49	fKeepDfromB(kFALSE),
	50	fKeepDfromBOnly(kFALSE),
	51	fmcLabel(0),
	52	fProngs(-1)
	53	{
	54	return;
	55	}
	56
	57
	58
	59	//_____________________________________________________
	60	AliCFVertexingHF::AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection) :
	61	fmcArray(mcArray),
	62	fRecoCandidate(0),
	63	fmcPartCandidate(0x0),
	64	fNDaughters(0),
65	fNVar(0),
66	fzPrimVertex(0),
67	fzMCVertex(0),
68	fFillFromGenerated(0),
69	fOriginDselection(0),
70	fKeepDfromB(kFALSE),
71	fKeepDfromBOnly(kFALSE),
72	fmcLabel(0),
73	fProngs(-1)
74	{
75
76
77	SetDselection(originDselection);
78
79
80	return;
81	}
82
83
84
85	//_______________________________________________________
86	AliCFVertexingHF::~AliCFVertexingHF()
87	{
88	if (fmcArray) delete fmcArray;
89	if (fRecoCandidate) delete fRecoCandidate;
90	if (fmcPartCandidate) delete fmcPartCandidate;
91	}
92
93
94	//_____________________________________________________
95	AliCFVertexingHF& AliCFVertexingHF::operator=(const AliCFVertexingHF& c){
96
97	if (this!= &c){
98	TObject::operator=(c);
99	fmcArray = c.fmcArray;
100	fRecoCandidate = c.fRecoCandidate;
101	fmcPartCandidate = c.fmcPartCandidate;
102	fNDaughters = c.fNDaughters;
103	fNVar = c.fNVar;
104	fzPrimVertex = c.fzPrimVertex;
105	fzMCVertex = c.fzMCVertex;
106	fFillFromGenerated = c.fFillFromGenerated;
107	fOriginDselection = c.fOriginDselection;
108	fKeepDfromB = c.fKeepDfromB;
109	fKeepDfromBOnly = c.fKeepDfromBOnly;
110	fmcLabel = c.fmcLabel;
111
112	}
113
114	return *this;
115	}
116
117	//____________________________________________________
118	AliCFVertexingHF::AliCFVertexingHF(const AliCFVertexingHF &c) :
119	TObject(c),
120	fmcArray(c.fmcArray),
121	fRecoCandidate(c.fRecoCandidate),
122	fmcPartCandidate(c.fmcPartCandidate),
123	fNDaughters(c.fNDaughters),
124	fNVar(c.fNVar),
125	fzPrimVertex(c.fzPrimVertex),
126	fzMCVertex(c.fzMCVertex),
127	fFillFromGenerated(c.fFillFromGenerated),
128	fOriginDselection (c.fOriginDselection),
129	fKeepDfromB (c.fKeepDfromB),
130	fKeepDfromBOnly (c.fKeepDfromBOnly),
131	fmcLabel(c.fmcLabel),
132	fProngs(c.fProngs)
133	{
134
135
136	}
137
138	//___________________________________________________________
139	void AliCFVertexingHF::SetDselection(UShort_t originDselection){
140
141	fOriginDselection = originDselection;
142
143	if (fOriginDselection == 0) {
144	fKeepDfromB = kFALSE;
145	fKeepDfromBOnly = kFALSE;
146	}
147
148	if (fOriginDselection == 1) {
149	fKeepDfromB = kTRUE;
150	fKeepDfromBOnly = kTRUE;
151	}
152
153	if (fOriginDselection == 2) {
154	fKeepDfromB = kTRUE;
155	fKeepDfromBOnly = kFALSE;
156	}
157
158	return;
159
160	}
161
162	//______________________________________________________
163	void AliCFVertexingHF::SetMCCandidateParam(Int_t label){
164
165	fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
166	fNDaughters = fmcPartCandidate->GetNDaughters();
167	return;
168	}
169
170
171	//____________________________________________________________
172	Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const{
173
174
175	Int_t cquarks = 0;
176	if (mcPart->GetPdgCode() == 4) cquarks++;
177	if (mcPart->GetPdgCode() == -4) cquarks++;
178	if (!mcPart) {
179	AliWarning("Particle not found in tree, skipping\n");
180	return cquarks;
181	}
182
183	return cquarks;
184	}
185
186
187	//________________________________________________________
188	Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle /mcPart/, TClonesArray /mcArray/) const {
189
190	Int_t pdgGranma = CheckOrigin();
191	if (pdgGranma == -9999){
192	printf ("This particle come from a B decay channel but the we keep only the prompt charm particles\n");
193	return kFALSE;
194	}
195
196	if (pdgGranma == -999){
197	printf ("This particle come from a prompt charm particles but we want only the ones coming from B\n");
198	return kFALSE;
199	}
200
201	if (!CheckMCDaughters()) return kFALSE;
202	if (!CheckMCChannelDecay()) return kFALSE;
203	return kTRUE;
204	}
205
206	//_________________________________________________________________________________________________
207	Int_t AliCFVertexingHF::CheckOrigin() const {
208	//
209	// checking whether the mother of the particles come from a charm or a bottom quark
210	//
211
212	Int_t pdgGranma = 0;
213	Int_t mother = 0;
214	mother = fmcPartCandidate->GetMother();
215	Int_t istep = 0;
216	while (mother >0 ){
217	istep++;
218	AliDebug(2,Form("mother at step %d = %d", istep, mother));
219	AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
220	pdgGranma = mcGranma->GetPdgCode();
221	AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
222	Int_t abspdgGranma = TMath::Abs(pdgGranma);
223	if ((abspdgGranma > 500 && abspdgGranma < 600) \|\| (abspdgGranma > 5000 && abspdgGranma < 6000)) {
224	if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
225
226	else{
227	return pdgGranma;
228	}
229	}
230	else {
231	if (fKeepDfromBOnly) return -999;
232	}
233
234	mother = mcGranma->GetMother();
235	}
236
237	return pdgGranma;
238	}
239
240
241	//___________________________________________
242	Bool_t AliCFVertexingHF::CheckMCDaughters()const {
243
244	AliAODMCParticle *mcPartDaughter;
245	Bool_t checkDaughters = kFALSE;
246
247	Int_t label0 = fmcPartCandidate->GetDaughter(0);
248	Int_t label1 = fmcPartCandidate->GetDaughter(1);
249	if (label1==0 \|\| label0 == 0){
250	AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
251	return checkDaughters;
252	}
253
254	if (label1 - label0 != fProngs-1){
255	AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
256	return checkDaughters;
257	}
258
259	for (Int_t iProng = 0; iProng<fProngs; iProng++){
260	mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
261	if (!mcPartDaughter) {
262	AliWarning("At least one Daughter Particle not found in tree, skipping");
263	return checkDaughters;
264	}
265	}
266
267	checkDaughters = kTRUE;
268	return checkDaughters;
269	}
270
271
272
273	//______________________________________________________
274	Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
275	{
276	// fill the container for Generator level selection
277	Bool_t mcContainerFilled = kFALSE;
278
279	Double_t* vectorMC = new Double_t[fNVar];
280	for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
281
282	if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
283	for (Int_t iVar = 0; iVar<fNVar; iVar++){
284
285	containerInputMC[iVar] = vectorMC[iVar];
286	}
287
288	mcContainerFilled = kTRUE;
289	}
290	delete vectorMC;
291	return mcContainerFilled;
292	}
293
294	//______________________________________________________
295	Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput) {
296
297	// fill the container for Reconstrucred level selection
298
299	Bool_t recoContainerFilled = kFALSE;
300	Double_t* vectorValues = new Double_t[fNVar];
301	Double_t* vectorReco = new Double_t[fNVar];
302
303	for (Int_t iVar = 0; iVar<fNVar; iVar++) {
304	vectorValues[iVar]= 9999.;
305	vectorReco[iVar]=9999.;
306	}
307
308	if(fFillFromGenerated){
309	//filled with MC values
310	if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
311	for (Int_t iVar = 0; iVar<fNVar; iVar++){
312	containerInput[iVar] = vectorValues[iVar];
313	}
314	recoContainerFilled = kTRUE;
315	}
316	}
317	else{
318	//filled with Reco values
319
320	if (GetRecoValuesFromCandidate(&vectorReco[0])){
321	for (Int_t iVar = 0; iVar<fNVar; iVar++){
322	containerInput[iVar] = vectorReco[iVar];
323	}
324	recoContainerFilled = kTRUE;
325	}
326	}
327
328	delete vectorValues;
329	delete vectorReco;
330	return recoContainerFilled;
331	}
332
333	//_____________________________________________________
334	Bool_t AliCFVertexingHF::MCAcceptanceStep() const
335	{
336
337	Bool_t bMCAccStep = kFALSE;
338
339	AliAODMCParticle *mcPartDaughter;
340	Int_t label0 = fmcPartCandidate->GetDaughter(0);
341	Int_t label1 = fmcPartCandidate->GetDaughter(1);
342	if (label1==0 \|\| label0 == 0){
343	AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
344	return bMCAccStep;
345	}
346
347	if (label1 - label0 != fProngs-1){
348	AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
349	return bMCAccStep;
350	}
351
352	for (Int_t iProng = 0; iProng<fProngs; iProng++){
353	mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
354	if (!mcPartDaughter) {
355	AliWarning("At least one Daughter Particle not found in tree, skipping");
356	return bMCAccStep;
357	}
358	Double_t eta = mcPartDaughter->Eta();
359	Double_t pt = mcPartDaughter->Pt();
360
361	//set values of eta and pt in the constructor.
362	if (TMath::Abs(eta) > 0.9 \|\| pt < 0.1){
363	AliDebug(3,"At least one daughter has eta>0.9 or pt < 0.1 \n");
364	return bMCAccStep;
365	}
366
367	}
368
369
370	bMCAccStep = kTRUE;
371	return bMCAccStep;
372
373	}
374
375	//_____________________________________________________
376	Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent aodEvent, AliESDtrackCuts trackCuts) const
377	{
378	// check on the kTPCrefit and kITSrefit conditions of the daughters
379	Bool_t bRefitStep = kFALSE;
380
381	Int_t label0 = fmcPartCandidate->GetDaughter(0);
382	Int_t label1 = fmcPartCandidate->GetDaughter(1);
383
384	if (label1==0 \|\| label0 == 0){
385	AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
386	return bRefitStep;
387	}
388
389	if (label1 - label0 != fProngs-1){
390	AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
391	//AliInfo(Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
392	return bRefitStep;
393	}
394
395	Int_t foundDaughters = 0;
396
397	if (trackCuts->GetRequireTPCRefit() \|\| trackCuts->GetRequireITSRefit()){
398
399	for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
400	AliAODTrack track = (AliAODTrack)aodEvent->GetTrack(iaod);
401	if (track->GetLabel()>= label0 && track->GetLabel()<=label1){
402	foundDaughters++;
403	printf("daughter %d \n",foundDaughters);
404	if(trackCuts->GetRequireTPCRefit()){
405	if((track->GetStatus()&AliESDtrack::kTPCrefit)) {
406	bRefitStep = kTRUE;
407	}
408	else {
409	AliDebug(3, "Refit cut not passed , missing TPC refit\n");
410	AliInfo( "Refit cut not passed , missing TPC refit\n");
411	return kFALSE;
412	}
413	}
414
415	if (trackCuts->GetRequireITSRefit()) {
416	if((track->GetStatus()&AliESDtrack::kITSrefit)){
417	bRefitStep = kTRUE;
418	}
419	else {
420	AliDebug(3, "Refit cut not passed , missing ITS refit\n");
421	AliInfo("Refit cut not passed , missing ITS refit\n");
422	return kFALSE;
423	}
424	}
425	}
426	if (foundDaughters == fProngs){
427
428	break;
429	}
430
431	}
432
433	}
434
435	return bRefitStep;
436	}
437
438	//____________________________________________________________________________
439
440	Bool_t AliCFVertexingHF::RecoStep()
441	//check also vertex and ITS Refit and TPC Refit
442	{
443	Bool_t bRecoStep = kFALSE;
444	Int_t mcLabel = GetMCLabel();
445
446	if (mcLabel == -1) {
447	AliDebug(2,"No MC particle found");
448	return bRecoStep;
449	}
450	else{
451	fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
452	if (!fmcPartCandidate){
453	AliWarning("Could not find associated MC in AOD MC tree");
454	return bRecoStep;
455	}
456
457	}
458
459	Int_t pdgGranma = CheckOrigin();
460
461	if (pdgGranma == -9999){
462	printf ("This particle come from a B decay channel but we keep only prompt charm particles\n");
463	return bRecoStep;
464	}
465
466	if (pdgGranma == -999){
467	printf ("This particle come from a prompt charm particle but we want only the ones coming from B\n");
468	return bRecoStep;
469	}
470
471
472	bRecoStep=kTRUE;
473	return bRecoStep;
474
475
476	}
477	//____________________________________________
478	Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const {
479
480	if (fRecoCandidate){
481	Double_t etaProng = fRecoCandidate->EtaProng(iProng);
482	return etaProng;
483	}
484	return 999999;
485	}
486	//______________________________________________________
487	Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const {
488
489	if (fRecoCandidate){
490	Double_t ptProng = fRecoCandidate->PtProng(iProng);
491	return ptProng;
492	}
493	return 999999;
494
495	}
496
497	//____________________________________________________________________
498
499	Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts *trackCuts) const
500	{
501	Bool_t bRecoAccStep = kFALSE;
502
503	Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
504	trackCuts->GetEtaRange(etaCutMin, etaCutMax);
505	trackCuts->GetPtRange(ptCutMin, ptCutMax);
506
507	Float_t etaProng=0., ptProng=0.;
508
509	for (Int_t iProng =0; iProng<fProngs; iProng++){
510
511	etaProng = GetEtaProng(iProng);
512	ptProng = GetPtProng(iProng);
513
514	Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
515	if (!acceptanceProng) {
516	printf ("At least one reconstructed prong isn't in the acceptance\n");
517	return bRecoAccStep;
518	}
519	}
520
521	bRecoAccStep=kTRUE;
522	return bRecoAccStep;
523	}
524	//___________________________________________________________
525
526	Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t *fill) const{
527
528	fill = new Double_t[4];
529
530	if(fmcPartCandidate){
531
532	fill[0] = GetPtCand();
533	fill[1] = GetYCand();
534
535	fill[2] = fmcPartCandidate->Pt();
536	fill[3] = fmcPartCandidate->Y();
537
538	return kTRUE;
539	}
540
541	delete fill;
542	return kFALSE;
543	}
544
545
546	//___________________________________________________________