[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpriVtx.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/
//
// QA class of primary vertex study for Heavy Flavor electrons
// this has functionality to reject electrons from primary vertex
// and check primary vertex characteristics
//  
// Authors:
//   MinJung Kweon <minjung@physi.uni-heidelberg.de>
//


#include <TH2F.h>
#include <TParticle.h>

#include "AliESDEvent.h"
#include <AliESDtrack.h>
#include <AliStack.h>

#include <AliLog.h>
#include <AliKFParticle.h>
#include <AliKFVertex.h>

#include "AliHFEpriVtx.h"


ClassImp(AliHFEpriVtx)

//_______________________________________________________________________________________________
AliHFEpriVtx::AliHFEpriVtx():
        fESD1(0x0)
        ,fStack(0x0)
        ,fNtrackswoPid(0)
        ,fHNtrackswoPid(0x0)
        ,fNESDprimVtxContributor(0x0)
        ,fNESDprimVtxIndices(0x0)
        ,fDiffDCAvsPt(0x0)
        ,fDiffDCAvsNt(0x0)
{ 
        //
        // Default constructor
        //

        Init();

}

//_______________________________________________________________________________________________
AliHFEpriVtx::AliHFEpriVtx(const AliHFEpriVtx &p):
         TObject(p)
        ,fESD1(0x0)
        ,fStack(0x0)
        ,fNtrackswoPid(p.fNtrackswoPid)
        ,fHNtrackswoPid(0x0)
        ,fNESDprimVtxContributor(0x0)
        ,fNESDprimVtxIndices(0x0)
        ,fDiffDCAvsPt(0x0)
        ,fDiffDCAvsNt(0x0)
{
        //
        // Copy constructor
        //
}

//_______________________________________________________________________________________________
AliHFEpriVtx&
AliHFEpriVtx::operator=(const AliHFEpriVtx &)
{
        //
        // Assignment operator
        //

        AliInfo("Not yet implemented.");
        return *this;
}

//_______________________________________________________________________________________________
AliHFEpriVtx::~AliHFEpriVtx()
{
        //
        // Destructor
        //

        AliInfo("Analysis Done.");
}

//__________________________________________
void AliHFEpriVtx::Init()
{
        //
        // initialize counters
        //

        fNtrackswoPid = 0;
        for (int i=0; i<10; i++){       
                fPrimVtx[i].fNtrackCount = 0 ;
                fPrimVtx[i].fNprimVtxContributorCount = 0 ;
        }
}

//_______________________________________________________________________________________________
void AliHFEpriVtx::CreateHistograms(TString hnopt)
{ 
        //
        // create histograms
        //

        fkSourceLabel[kAll]="all";
        fkSourceLabel[kDirectCharm]="directCharm";
        fkSourceLabel[kDirectBeauty]="directBeauty";
        fkSourceLabel[kBeautyCharm]="beauty2charm";
        fkSourceLabel[kGamma]="gamma";
        fkSourceLabel[kPi0]="pi0";
        fkSourceLabel[kElse]="others";
        fkSourceLabel[kBeautyGamma]="beauty22gamma";
        fkSourceLabel[kBeautyPi0]="beauty22pi0";
        fkSourceLabel[kBeautyElse]="beauty22others";


        TString hname;
        for (Int_t isource = 0; isource < 10; isource++ ){

           hname=hnopt+"ntracks_"+fkSourceLabel[isource];
           fPrimVtx[isource].fNtracks = new TH1F(hname,hname,50,0,50);
           hname=hnopt+"nPrimVtxContributor_"+fkSourceLabel[isource];
           fPrimVtx[isource].fNprimVtxContributor = new TH1F(hname,hname,100,0,100);
           hname=hnopt+"PtElec_"+fkSourceLabel[isource];
           fPrimVtx[isource].fPtElec = new TH1F(hname,hname,250,0,50);
           hname=hnopt+"PtElecContributor_"+fkSourceLabel[isource];
           fPrimVtx[isource].fPtElecContributor = new TH1F(hname,hname,250,0,50);

        }

        hname=hnopt+"ntrackswopid";
        fHNtrackswoPid = new TH1F(hname,hname,50,0,50);
        hname=hnopt+"nESDprimVtxContributor";
        fNESDprimVtxContributor = new TH1I(hname,hname,100,0,100);
        hname=hnopt+"nESDprimVtxIndices";
        fNESDprimVtxIndices= new TH1I(hname,hname,100,0,100);
        hname=hnopt+"diffDCAvsPt";
        fDiffDCAvsPt = new TH2F(hname,hname,250,0,50,500,0,1);
        hname=hnopt+"diffDCAvsNt";
        fDiffDCAvsNt = new TH2F(hname,hname,100,0,100,500,0,1);

}

//_______________________________________________________________________________________________
void AliHFEpriVtx::CountNtracks(Int_t sourcePart, Int_t recpid, Double_t recprob)
{
        //
        // count number of tracks passed certain cuts
        //

        fNtrackswoPid++;

        if (!recpid && recprob>0.5)     
        fPrimVtx[kAll].fNtrackCount++;
        if(sourcePart<0) return;
        fPrimVtx[sourcePart].fNtrackCount++;

}

//_______________________________________________________________________________________________
void AliHFEpriVtx::FillNtracks()
{
        //
        // count number of tracks passed certain cuts
        //

        fHNtrackswoPid->Fill(fNtrackswoPid);
        for (int i=0; i<10; i++){
          fPrimVtx[i].fNtracks->Fill(fPrimVtx[i].fNtrackCount);
        }

}

//_______________________________________________________________________________________________
Int_t AliHFEpriVtx::GetMCPID(AliESDtrack *track) 
{
        //
        // get MC pid
        //

        Int_t label = TMath::Abs(track->GetLabel());
        TParticle* mcpart = fStack->Particle(label);
        if ( !mcpart ) return 0;
        Int_t pdgCode = mcpart->GetPdgCode();

        return pdgCode;
}

//_______________________________________________________________________________________________
void AliHFEpriVtx::GetNPriVxtContributor() 
{
        //
        // count number of primary vertex contributor
        //

        const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
        fNESDprimVtxContributor->Fill(primvtx->GetNContributors());
        fNESDprimVtxIndices->Fill(primvtx->GetNIndices());
}

//_______________________________________________________________________________________________
void AliHFEpriVtx::CountPriVxtElecContributor(AliESDtrack *ESDelectron, Int_t sourcePart, Int_t recpid, Double_t recprob) 
{
        //
        // count number of electrons contributing to the primary vertex
        //


        if (recpid || recprob<0.5) return;

        // get track id of our selected electron
        Int_t elecTrkID = ESDelectron->GetID();

        Int_t label = TMath::Abs(ESDelectron->GetLabel());
        TParticle* mcpart = fStack->Particle(label);

        AliKFParticle::SetField(fESD1->GetMagneticField());
        AliKFParticle kfElectron(*ESDelectron,11);
 
        // prepare kfprimary vertex
        AliKFVertex kfESDprimary;

        // Reconstructed Primary Vertex (with ESD tracks)
        const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
//         Int_t nt = primvtx->GetNContributors();
        Int_t n=primvtx->GetNIndices();

        if (n>0 && primvtx->GetStatus()){

                kfESDprimary = AliKFVertex(*primvtx);
                Double_t dcaBFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);

                UShort_t *priIndex = primvtx->GetIndices();

                fPrimVtx[kAll].fPtElec->Fill(mcpart->Pt());
                if(sourcePart>=0) fPrimVtx[sourcePart].fPtElec->Fill(mcpart->Pt());

                for (Int_t i=0;i<n;i++){

                        Int_t idx = Int_t(priIndex[i]);
                        if (idx == elecTrkID){
                                fPrimVtx[kAll].fNprimVtxContributorCount++;
                                fPrimVtx[kAll].fPtElecContributor->Fill(mcpart->Pt());
                                if(sourcePart<0) continue;
                                fPrimVtx[sourcePart].fNprimVtxContributorCount++;
                                fPrimVtx[sourcePart].fPtElecContributor->Fill(mcpart->Pt());

                                kfESDprimary -= kfElectron;
                                Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
                                fDiffDCAvsPt->Fill(mcpart->Pt(),dcaBFrmElec-dcaAFrmElec);
                                fDiffDCAvsNt->Fill(n,dcaBFrmElec-dcaAFrmElec);
                        }
                } 
        }  

}

//_______________________________________________________________________________________________
void AliHFEpriVtx::FillNprimVtxContributor() const
{
        //
        // Fill histogram with number of electrons contributing to the primary vertex
        //

        for (int i=0; i<10; i++){
          fPrimVtx[i].fNprimVtxContributor->Fill(fPrimVtx[i].fNprimVtxContributorCount);
        }

}

//_______________________________________________________________________________________________
Double_t AliHFEpriVtx::GetDistanceFromRecalVertexXY(AliESDtrack *ESDelectron) 
{
        //
        // return recalculated DCA after removing input track from the primary vertex
        //

        // get track id of our selected electron
        Int_t elecTrkID = ESDelectron->GetID();

        AliKFParticle::SetField(fESD1->GetMagneticField());
        AliKFParticle kfElectron(*ESDelectron,11);
 
        // prepare kfprimary vertex
        AliKFVertex kfESDprimary;

        // Reconstructed Primary Vertex (with ESD tracks)
        const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
        Int_t n=primvtx->GetNIndices();

        if (n>0 && primvtx->GetStatus()){

                kfESDprimary = AliKFVertex(*primvtx);
                UShort_t *priIndex = primvtx->GetIndices();

                for (Int_t i=0;i<n;i++){

                        Int_t idx = Int_t(priIndex[i]);
                        if (idx == elecTrkID){
                                kfESDprimary -= kfElectron;
                                Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);

                                return dcaAFrmElec;
                        }
                } 
        }  
	return -1;

}
Commit	Line	Data
259c3296	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	**************************************************************************/
50685501	15	//
9bcfd1ab	16	// QA class of primary vertex study for Heavy Flavor electrons
	17	// this has functionality to reject electrons from primary vertex
	18	// and check primary vertex characteristics
	19	//
50685501	20	// Authors:
	21	// MinJung Kweon <minjung@physi.uni-heidelberg.de>
	22	//
259c3296	23
259c3296	24
259c3296	25	#include <TH2F.h>
259c3296	26	#include <TParticle.h>
	27
	28	#include "AliESDEvent.h"
	29	#include <AliESDtrack.h>
	30	#include <AliStack.h>
	31
	32	#include <AliLog.h>
	33	#include <AliKFParticle.h>
	34	#include <AliKFVertex.h>
	35
	36	#include "AliHFEpriVtx.h"
	37
	38
	39	ClassImp(AliHFEpriVtx)
	40
	41	//_______________________________________________________________________________________________
	42	AliHFEpriVtx::AliHFEpriVtx():
	43	fESD1(0x0)
	44	,fStack(0x0)
	45	,fNtrackswoPid(0)
	46	,fHNtrackswoPid(0x0)
	47	,fNESDprimVtxContributor(0x0)
	48	,fNESDprimVtxIndices(0x0)
	49	,fDiffDCAvsPt(0x0)
	50	,fDiffDCAvsNt(0x0)
	51	{
	52	//
	53	// Default constructor
	54	//
	55
	56	Init();
	57
	58	}
	59
	60	//_______________________________________________________________________________________________
	61	AliHFEpriVtx::AliHFEpriVtx(const AliHFEpriVtx &p):
	62	TObject(p)
	63	,fESD1(0x0)
	64	,fStack(0x0)
	65	,fNtrackswoPid(p.fNtrackswoPid)
	66	,fHNtrackswoPid(0x0)
	67	,fNESDprimVtxContributor(0x0)
	68	,fNESDprimVtxIndices(0x0)
	69	,fDiffDCAvsPt(0x0)
	70	,fDiffDCAvsNt(0x0)
	71	{
	72	//
	73	// Copy constructor
	74	//
	75	}
	76
	77	//_______________________________________________________________________________________________
	78	AliHFEpriVtx&
	79	AliHFEpriVtx::operator=(const AliHFEpriVtx &)
	80	{
	81	//
	82	// Assignment operator
	83	//
	84
	85	AliInfo("Not yet implemented.");
	86	return *this;
	87	}
	88
	89	//_______________________________________________________________________________________________
90	AliHFEpriVtx::~AliHFEpriVtx()
91	{
92	//
93	// Destructor
94	//
95
96	AliInfo("Analysis Done.");
97	}
98
99	//__________________________________________
100	void AliHFEpriVtx::Init()
101	{
102	//
103	// initialize counters
104	//
105
106	fNtrackswoPid = 0;
107	for (int i=0; i<10; i++){
108	fPrimVtx[i].fNtrackCount = 0 ;
109	fPrimVtx[i].fNprimVtxContributorCount = 0 ;
110	}
111	}
112
113	//_______________________________________________________________________________________________
114	void AliHFEpriVtx::CreateHistograms(TString hnopt)
115	{
116	//
117	// create histograms
118	//
119
120	fkSourceLabel[kAll]="all";
121	fkSourceLabel[kDirectCharm]="directCharm";
122	fkSourceLabel[kDirectBeauty]="directBeauty";
123	fkSourceLabel[kBeautyCharm]="beauty2charm";
124	fkSourceLabel[kGamma]="gamma";
125	fkSourceLabel[kPi0]="pi0";
126	fkSourceLabel[kElse]="others";
127	fkSourceLabel[kBeautyGamma]="beauty22gamma";
128	fkSourceLabel[kBeautyPi0]="beauty22pi0";
129	fkSourceLabel[kBeautyElse]="beauty22others";
130
131
132	TString hname;
133	for (Int_t isource = 0; isource < 10; isource++ ){
134
135	hname=hnopt+"ntracks_"+fkSourceLabel[isource];
136	fPrimVtx[isource].fNtracks = new TH1F(hname,hname,50,0,50);
137	hname=hnopt+"nPrimVtxContributor_"+fkSourceLabel[isource];
138	fPrimVtx[isource].fNprimVtxContributor = new TH1F(hname,hname,100,0,100);
139	hname=hnopt+"PtElec_"+fkSourceLabel[isource];
78ea5ef4	140	fPrimVtx[isource].fPtElec = new TH1F(hname,hname,250,0,50);
259c3296	141	hname=hnopt+"PtElecContributor_"+fkSourceLabel[isource];
78ea5ef4	142	fPrimVtx[isource].fPtElecContributor = new TH1F(hname,hname,250,0,50);
259c3296	143
	144	}
	145
	146	hname=hnopt+"ntrackswopid";
	147	fHNtrackswoPid = new TH1F(hname,hname,50,0,50);
	148	hname=hnopt+"nESDprimVtxContributor";
	149	fNESDprimVtxContributor = new TH1I(hname,hname,100,0,100);
	150	hname=hnopt+"nESDprimVtxIndices";
	151	fNESDprimVtxIndices= new TH1I(hname,hname,100,0,100);
	152	hname=hnopt+"diffDCAvsPt";
78ea5ef4	153	fDiffDCAvsPt = new TH2F(hname,hname,250,0,50,500,0,1);
259c3296	154	hname=hnopt+"diffDCAvsNt";
	155	fDiffDCAvsNt = new TH2F(hname,hname,100,0,100,500,0,1);
	156
	157	}
	158
	159	//_______________________________________________________________________________________________
	160	void AliHFEpriVtx::CountNtracks(Int_t sourcePart, Int_t recpid, Double_t recprob)
	161	{
	162	//
	163	// count number of tracks passed certain cuts
	164	//
	165
	166	fNtrackswoPid++;
	167
	168	if (!recpid && recprob>0.5)
	169	fPrimVtx[kAll].fNtrackCount++;
	170	if(sourcePart<0) return;
	171	fPrimVtx[sourcePart].fNtrackCount++;
	172
	173	}
	174
	175	//_______________________________________________________________________________________________
	176	void AliHFEpriVtx::FillNtracks()
	177	{
	178	//
	179	// count number of tracks passed certain cuts
	180	//
	181
	182	fHNtrackswoPid->Fill(fNtrackswoPid);
	183	for (int i=0; i<10; i++){
	184	fPrimVtx[i].fNtracks->Fill(fPrimVtx[i].fNtrackCount);
	185	}
	186
	187	}
	188
	189	//_______________________________________________________________________________________________
	190	Int_t AliHFEpriVtx::GetMCPID(AliESDtrack *track)
	191	{
	192	//
	193	// get MC pid
	194	//
	195
	196	Int_t label = TMath::Abs(track->GetLabel());
	197	TParticle* mcpart = fStack->Particle(label);
	198	if ( !mcpart ) return 0;
	199	Int_t pdgCode = mcpart->GetPdgCode();
	200
	201	return pdgCode;
	202	}
	203
	204	//_______________________________________________________________________________________________
	205	void AliHFEpriVtx::GetNPriVxtContributor()
	206	{
	207	//
	208	// count number of primary vertex contributor
	209	//
	210
	211	const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
	212	fNESDprimVtxContributor->Fill(primvtx->GetNContributors());
	213	fNESDprimVtxIndices->Fill(primvtx->GetNIndices());
	214	}
	215
	216	//_______________________________________________________________________________________________
	217	void AliHFEpriVtx::CountPriVxtElecContributor(AliESDtrack *ESDelectron, Int_t sourcePart, Int_t recpid, Double_t recprob)
218	{
219	//
220	// count number of electrons contributing to the primary vertex
221	//
222
223
224	if (recpid \|\| recprob<0.5) return;
225
226	// get track id of our selected electron
227	Int_t elecTrkID = ESDelectron->GetID();
228
229	Int_t label = TMath::Abs(ESDelectron->GetLabel());
230	TParticle* mcpart = fStack->Particle(label);
231
232	AliKFParticle::SetField(fESD1->GetMagneticField());
233	AliKFParticle kfElectron(*ESDelectron,11);
234
235	// prepare kfprimary vertex
236	AliKFVertex kfESDprimary;
237
238	// Reconstructed Primary Vertex (with ESD tracks)
239	const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
240	// Int_t nt = primvtx->GetNContributors();
241	Int_t n=primvtx->GetNIndices();
242
243	if (n>0 && primvtx->GetStatus()){
244
245	kfESDprimary = AliKFVertex(*primvtx);
246	Double_t dcaBFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
247
248	UShort_t *priIndex = primvtx->GetIndices();
249
250	fPrimVtx[kAll].fPtElec->Fill(mcpart->Pt());
251	if(sourcePart>=0) fPrimVtx[sourcePart].fPtElec->Fill(mcpart->Pt());
252
253	for (Int_t i=0;i<n;i++){
254
255	Int_t idx = Int_t(priIndex[i]);
256	if (idx == elecTrkID){
257	fPrimVtx[kAll].fNprimVtxContributorCount++;
258	fPrimVtx[kAll].fPtElecContributor->Fill(mcpart->Pt());
259	if(sourcePart<0) continue;
260	fPrimVtx[sourcePart].fNprimVtxContributorCount++;
261	fPrimVtx[sourcePart].fPtElecContributor->Fill(mcpart->Pt());
262
263	kfESDprimary -= kfElectron;
264	Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
265	fDiffDCAvsPt->Fill(mcpart->Pt(),dcaBFrmElec-dcaAFrmElec);
266	fDiffDCAvsNt->Fill(n,dcaBFrmElec-dcaAFrmElec);
267	}
268	}
269	}
270
271	}
272
273	//_______________________________________________________________________________________________
75d81601	274	void AliHFEpriVtx::FillNprimVtxContributor() const
259c3296	275	{
	276	//
	277	// Fill histogram with number of electrons contributing to the primary vertex
	278	//
	279
	280	for (int i=0; i<10; i++){
	281	fPrimVtx[i].fNprimVtxContributor->Fill(fPrimVtx[i].fNprimVtxContributorCount);
	282	}
	283
	284	}
70da6c5a	285
	286	//_______________________________________________________________________________________________
	287	Double_t AliHFEpriVtx::GetDistanceFromRecalVertexXY(AliESDtrack *ESDelectron)
	288	{
	289	//
	290	// return recalculated DCA after removing input track from the primary vertex
	291	//
	292
	293	// get track id of our selected electron
	294	Int_t elecTrkID = ESDelectron->GetID();
	295
	296	AliKFParticle::SetField(fESD1->GetMagneticField());
	297	AliKFParticle kfElectron(*ESDelectron,11);
	298
	299	// prepare kfprimary vertex
	300	AliKFVertex kfESDprimary;
	301
	302	// Reconstructed Primary Vertex (with ESD tracks)
	303	const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
	304	Int_t n=primvtx->GetNIndices();
	305
	306	if (n>0 && primvtx->GetStatus()){
	307
	308	kfESDprimary = AliKFVertex(*primvtx);
	309	UShort_t *priIndex = primvtx->GetIndices();
	310
	311	for (Int_t i=0;i<n;i++){
	312
	313	Int_t idx = Int_t(priIndex[i]);
	314	if (idx == elecTrkID){
	315	kfESDprimary -= kfElectron;
	316	Double_t dcaAFrmElec = kfElectron.GetDistanceFromVertexXY(kfESDprimary);
	317
	318	return dcaAFrmElec;
	319	}
	320	}
	321	}
	322	return -1;
	323
	324	}