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

/**************************************************************************
 * Copyright(c) 1998-2008, 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.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// AliAnalysisTaskSE for the comparison of heavy flavor
// decay candidates with the MC truth.
//
// Author: A.Dainese, andrea.dainese@lnl.infn.it
/////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TNtuple.h>
#include <TList.h>
#include <TH1F.h>

#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSECompareHF.h"

ClassImp(AliAnalysisTaskSECompareHF)


//________________________________________________________________________
AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
AliAnalysisTaskSE(),
fOutput(0), 
fNtupleD0Cmp(0),
fHistMass(0)
{
  // Default constructor
  SetD0toKpiCuts();

  // Output slot #1 writes into a TList container
  DefineOutput(1,TList::Class());  //My private output
}

//________________________________________________________________________
AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
AliAnalysisTaskSE(name),
fOutput(0), 
fNtupleD0Cmp(0),
fHistMass(0)
{
  // Default constructor
  SetD0toKpiCuts();

  // Output slot #1 writes into a TList container
  DefineOutput(1,TList::Class());  //My private output
}

//________________________________________________________________________
AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
{
  // Destructor
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}  

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::Init()
{
  // Initialization

  if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
{
  // Create the output container
  //
  if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList();
  fOutput->SetOwner();

  fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
  fHistMass->Sumw2();
  fHistMass->SetMinimum(0);
  fOutput->Add(fHistMass);

  fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
  fOutput->Add(fNtupleD0Cmp);

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event:
  // heavy flavor candidates association to MC truth
  
  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());

  // load D0->Kpi candidates                                                   
  TClonesArray *inputArrayD0toKpi =
    (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
  if(!inputArrayD0toKpi) {
    printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
    return;
  }

  // AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  //vtx1->Print();

  // load MC particles
  TClonesArray *mcArray = 
    (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
  if(!mcArray) {
    printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
    return;
  }

  // load MC header
  AliAODMCHeader *mcHeader = 
    (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
  if(!mcHeader) {
    printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
    return;
  }
  
    
  // loop over D0->Kpi candidates
  Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
  printf("Number of D0->Kpi: %d\n",nInD0toKpi);

  Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0;
  
  for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
    AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
    labD0daugh0=-1; 
    labD0daugh1=-1;
    Bool_t unsetvtx=kFALSE;
    if(!d->GetOwnPrimaryVtx()) {
      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
      unsetvtx=kTRUE;
    }
    Int_t okD0=0,okD0bar=0; 
    if(d->SelectD0(fD0toKpiCuts,okD0,okD0bar)) {
      // get daughter AOD tracks
      AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
      AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);

      lab0 = trk0->GetLabel();
      lab1 = trk1->GetLabel();


      AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(lab0);
      if(!part0) { 
	printf("no MC particle\n");
	continue;
      }
      while(part0->GetMother()>=0) {
	labMother=part0->GetMother();
	part0 = (AliAODMCParticle*)mcArray->At(labMother);
	if(!part0) {
	  printf("no MC mother particle\n");
	  break;
	}
	pdgMother = TMath::Abs(part0->GetPdgCode());
	if(pdgMother==421) {
	  labD0daugh0=labMother;
	  break;
	}
      }

      AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(lab1);
      if(!part1) {
	printf("no MC particle\n");
	continue;
      }
      while(part1->GetMother()>=0) {
	labMother=part1->GetMother();
	part1 = (AliAODMCParticle*)mcArray->At(labMother);
	if(!part1) {
	  printf("no MC mother particle\n");
	  break;
	}
	pdgMother = TMath::Abs(part1->GetPdgCode());
	if(pdgMother==421) {
	  labD0daugh1=labMother;
	  break;
	}
      }

      // check if the candidate is signal
      if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) {

	AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0daugh0);
	pdgD0 = partD0->GetPdgCode();
	Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
	fHistMass->Fill(invmass);

	// Post the data already here
	PostData(1,fOutput);

	fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
      }


    }
    if(unsetvtx) d->UnsetOwnPrimaryVtx();
  } // end loop on D0->Kpi


  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");

  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    printf("ERROR: fOutput not available\n");
    return;
  }

  fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
  fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
				   Double_t cut2,Double_t cut3,Double_t cut4,
				   Double_t cut5,Double_t cut6,
				   Double_t cut7,Double_t cut8) 
{
  // Set the cuts for D0 selection
  // cuts[0] = inv. mass half width [GeV]   
  // cuts[1] = dca [cm]
  // cuts[2] = cosThetaStar 
  // cuts[3] = pTK [GeV/c]
  // cuts[4] = pTPi [GeV/c]
  // cuts[5] = d0K [cm]   upper limit!
  // cuts[6] = d0Pi [cm]  upper limit!
  // cuts[7] = d0d0 [cm^2]
  // cuts[8] = cosThetaPoint

  fD0toKpiCuts[0] = cut0;
  fD0toKpiCuts[1] = cut1;
  fD0toKpiCuts[2] = cut2;
  fD0toKpiCuts[3] = cut3;
  fD0toKpiCuts[4] = cut4;
  fD0toKpiCuts[5] = cut5;
  fD0toKpiCuts[6] = cut6;
  fD0toKpiCuts[7] = cut7;
  fD0toKpiCuts[8] = cut8;

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(const Double_t cuts[9]) 
{
  // Set the cuts for D0 selection
  // cuts[0] = inv. mass half width [GeV]   
  // cuts[1] = dca [cm]
  // cuts[2] = cosThetaStar 
  // cuts[3] = pTK [GeV/c]
  // cuts[4] = pTPi [GeV/c]
  // cuts[5] = d0K [cm]   upper limit!
  // cuts[6] = d0Pi [cm]  upper limit!
  // cuts[7] = d0d0 [cm^2]
  // cuts[8] = cosThetaPoint

  for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];

  return;
}
Commit	Line	Data
3a219f60	1	/**************************************************************************
	2	* Copyright(c) 1998-2008, 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	//
	18	// AliAnalysisTaskSE for the comparison of heavy flavor
	19	// decay candidates with the MC truth.
	20	//
	21	// Author: A.Dainese, andrea.dainese@lnl.infn.it
	22	/////////////////////////////////////////////////////////////
	23
	24	#include <TClonesArray.h>
	25	#include <TNtuple.h>
	26	#include <TList.h>
	27	#include <TH1F.h>
	28
	29	#include "AliAODEvent.h"
	30	#include "AliAODVertex.h"
	31	#include "AliAODTrack.h"
	32	#include "AliAODMCHeader.h"
	33	#include "AliAODMCParticle.h"
	34	#include "AliAODRecoDecayHF2Prong.h"
	35	#include "AliAnalysisTaskSE.h"
	36	#include "AliAnalysisTaskSECompareHF.h"
	37
	38	ClassImp(AliAnalysisTaskSECompareHF)
	39
	40
	41	//________________________________________________________________________
	42	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
	43	AliAnalysisTaskSE(),
	44	fOutput(0),
	45	fNtupleD0Cmp(0),
	46	fHistMass(0)
	47	{
	48	// Default constructor
	49	SetD0toKpiCuts();
	50
	51	// Output slot #1 writes into a TList container
	52	DefineOutput(1,TList::Class()); //My private output
	53	}
	54
	55	//________________________________________________________________________
	56	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
	57	AliAnalysisTaskSE(name),
	58	fOutput(0),
	59	fNtupleD0Cmp(0),
	60	fHistMass(0)
	61	{
	62	// Default constructor
	63	SetD0toKpiCuts();
	64
65	// Output slot #1 writes into a TList container
66	DefineOutput(1,TList::Class()); //My private output
67	}
68
69	//________________________________________________________________________
70	AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
71	{
72	// Destructor
73	if (fOutput) {
74	delete fOutput;
75	fOutput = 0;
76	}
77	}
78
79	//________________________________________________________________________
80	void AliAnalysisTaskSECompareHF::Init()
81	{
82	// Initialization
83
84	if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
85
86	return;
87	}
88
89	//________________________________________________________________________
90	void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
91	{
92	// Create the output container
93	//
94	if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
95
96	// Several histograms are more conveniently managed in a TList
97	fOutput = new TList();
98	fOutput->SetOwner();
99
100	fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
101	fHistMass->Sumw2();
102	fHistMass->SetMinimum(0);
103	fOutput->Add(fHistMass);
104
105	fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
106	fOutput->Add(fNtupleD0Cmp);
107
108	return;
109	}
110
111	//________________________________________________________________________
112	void AliAnalysisTaskSECompareHF::UserExec(Option_t /option*/)
113	{
114	// Execute analysis for current event:
115	// heavy flavor candidates association to MC truth
116
117	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
118
119	// load D0->Kpi candidates
120	TClonesArray *inputArrayD0toKpi =
121	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
122	if(!inputArrayD0toKpi) {
123	printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
124	return;
125	}
126
127	// AOD primary vertex
128	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
129	//vtx1->Print();
130
131	// load MC particles
132	TClonesArray *mcArray =
133	(TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
134	if(!mcArray) {
135	printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
136	return;
137	}
138
139	// load MC header
140	AliAODMCHeader *mcHeader =
141	(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
142	if(!mcHeader) {
143	printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
144	return;
145	}
146
147
148
149	// loop over D0->Kpi candidates
150	Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
151	printf("Number of D0->Kpi: %d\n",nInD0toKpi);
152
153	Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0;
154
155	for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
156	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
157	labD0daugh0=-1;
158	labD0daugh1=-1;
159	Bool_t unsetvtx=kFALSE;
160	if(!d->GetOwnPrimaryVtx()) {
161	d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
162	unsetvtx=kTRUE;
163	}
164	Int_t okD0=0,okD0bar=0;
165	if(d->SelectD0(fD0toKpiCuts,okD0,okD0bar)) {
166	// get daughter AOD tracks
167	AliAODTrack trk0 = (AliAODTrack)d->GetDaughter(0);
168	AliAODTrack trk1 = (AliAODTrack)d->GetDaughter(1);
169
170	lab0 = trk0->GetLabel();
171	lab1 = trk1->GetLabel();
172
173
174	AliAODMCParticle part0 = (AliAODMCParticle)mcArray->At(lab0);
175	if(!part0) {
176	printf("no MC particle\n");
177	continue;
178	}
179	while(part0->GetMother()>=0) {
180	labMother=part0->GetMother();
181	part0 = (AliAODMCParticle*)mcArray->At(labMother);
182	if(!part0) {
183	printf("no MC mother particle\n");
184	break;
185	}
186	pdgMother = TMath::Abs(part0->GetPdgCode());
187	if(pdgMother==421) {
188	labD0daugh0=labMother;
189	break;
190	}
191	}
192
193	AliAODMCParticle part1 = (AliAODMCParticle)mcArray->At(lab1);
194	if(!part1) {
195	printf("no MC particle\n");
196	continue;
197	}
198	while(part1->GetMother()>=0) {
199	labMother=part1->GetMother();
200	part1 = (AliAODMCParticle*)mcArray->At(labMother);
201	if(!part1) {
202	printf("no MC mother particle\n");
203	break;
204	}
205	pdgMother = TMath::Abs(part1->GetPdgCode());
206	if(pdgMother==421) {
207	labD0daugh1=labMother;
208	break;
209	}
210	}
211
212	// check if the candidate is signal
213	if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) {
214
215	AliAODMCParticle partD0 = (AliAODMCParticle)mcArray->At(labD0daugh0);
216	pdgD0 = partD0->GetPdgCode();
217	Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
218	fHistMass->Fill(invmass);
219
220	// Post the data already here
221	PostData(1,fOutput);
222
223	fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
224	}
225
226
227	}
228	if(unsetvtx) d->UnsetOwnPrimaryVtx();
229	} // end loop on D0->Kpi
230
231
232	return;
233	}
234
235	//________________________________________________________________________
236	void AliAnalysisTaskSECompareHF::Terminate(Option_t /option*/)
237	{
238	// Terminate analysis
239	//
240	if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
241
242	fOutput = dynamic_cast<TList*> (GetOutputData(1));
243	if (!fOutput) {
244	printf("ERROR: fOutput not available\n");
245	return;
246	}
247
248	fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
249	fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
250
251	return;
252	}
253
254	//________________________________________________________________________
255	void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
256	Double_t cut2,Double_t cut3,Double_t cut4,
257	Double_t cut5,Double_t cut6,
258	Double_t cut7,Double_t cut8)
259	{
260	// Set the cuts for D0 selection
261	// cuts[0] = inv. mass half width [GeV]
262	// cuts[1] = dca [cm]
263	// cuts[2] = cosThetaStar
264	// cuts[3] = pTK [GeV/c]
265	// cuts[4] = pTPi [GeV/c]
266	// cuts[5] = d0K [cm] upper limit!
267	// cuts[6] = d0Pi [cm] upper limit!
268	// cuts[7] = d0d0 [cm^2]
269	// cuts[8] = cosThetaPoint
270
271	fD0toKpiCuts[0] = cut0;
272	fD0toKpiCuts[1] = cut1;
273	fD0toKpiCuts[2] = cut2;
274	fD0toKpiCuts[3] = cut3;
275	fD0toKpiCuts[4] = cut4;
276	fD0toKpiCuts[5] = cut5;
277	fD0toKpiCuts[6] = cut6;
278	fD0toKpiCuts[7] = cut7;
279	fD0toKpiCuts[8] = cut8;
280
281	return;
282	}
283
284	//________________________________________________________________________
285	void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(const Double_t cuts[9])
286	{
287	// Set the cuts for D0 selection
288	// cuts[0] = inv. mass half width [GeV]
289	// cuts[1] = dca [cm]
290	// cuts[2] = cosThetaStar
291	// cuts[3] = pTK [GeV/c]
292	// cuts[4] = pTPi [GeV/c]
293	// cuts[5] = d0K [cm] upper limit!
294	// cuts[6] = d0Pi [cm] upper limit!
295	// cuts[7] = d0d0 [cm^2]
296	// cuts[8] = cosThetaPoint
297
298	for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
299
300	return;
301	}