[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisElectronTask.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.                  *
**************************************************************************/
/*
 * The analysis task:
 * Filling an AliCFContainer with the quantities pt, eta and phi
 * for tracks which survivied the particle cuts (MC resp. ESD tracks)
 * Track selection is done using the AliHFE package
 * 
 * Author:
 *  Markus Fasel <M.Fasel@gsi.de>
 */
#include <TAxis.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TH1F.h>
#include <TH1I.h>
#include <TH2F.h>
#include <TIterator.h>
#include <TList.h>
#include <TLegend.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TParticle.h>
#include <TProfile.h>
#include <TTree.h>

#include "AliCFContainer.h"
#include "AliCFManager.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliLog.h"
#include "AliAnalysisManager.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMCParticle.h"
#include "AliPID.h"

#include "AliHFEpid.h"
#include "AliHFEcuts.h"
#include "AliAnalysisElectronTask.h"

//____________________________________________________________
AliAnalysisElectronTask::AliAnalysisElectronTask():
	AliAnalysisTask("PID efficiency Analysis", "")
	, fESD(0x0)
	, fMC(0x0)
	, fCFM(0x0)
	, fPID(0x0)
  , fCuts(0x0)
	, fNEvents(0x0)
	, fQA(0x0)
{
	DefineInput(0, TChain::Class());
	DefineOutput(0, TH1I::Class());
	DefineOutput(1, AliCFContainer::Class());
	DefineOutput(2, TList::Class());

  // Initialize cuts
  fCuts = new AliHFEcuts;
  fPID = new AliHFEpid;
}

//____________________________________________________________
AliAnalysisElectronTask::~AliAnalysisElectronTask(){
	if(fESD) delete fESD;
	if(fMC) delete fMC;
	if(fPID) delete fPID;
	if(fQA) delete fQA;
  if(fCuts) delete fCuts;
	if(fNEvents) delete fNEvents;
	fQA = 0x0; fMC = 0x0; fESD = 0x0; fCFM = 0x0; fPID = 0x0; fCuts = 0x0;
}

//____________________________________________________________
void AliAnalysisElectronTask::ConnectInputData(Option_t *){
	TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
	if(!esdchain){
		AliError("ESD chain empty");
		return;
	} else {
		esdchain->SetBranchStatus("Tracks", 1);
	}
	AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	if(!esdH){	
		AliError("No ESD input handler");
		return;
	} else {
		fESD = esdH->GetEvent();
	}
	AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	if(!mcH){	
		AliError("No MC truth handler");
		return;
	} else {
		fMC = mcH->MCEvent();
	}
}

//____________________________________________________________
void AliAnalysisElectronTask::CreateOutputObjects(){
	fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
	// First Step: TRD alone
	if(!fQA) fQA = new TList;
	fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
	fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
	fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
	fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
	fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
	fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
	fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);

  // Initialize correction Framework and Cuts
  fCFM = new AliCFManager;
  MakeParticleContainer();
  // Temporary fix: Initialize particle cuts with 0x0
  for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
    fCFM->SetParticleCutsList(istep, 0x0);
  if(IsQAOn()){
    fCuts->SetDebugMode();
    fQA->AddAt(fCuts->GetQAhistograms(), 7);
  }
  fCuts->CreateStandardCuts();
  fCuts->Initialize(fCFM);

  // Initialize PID
  if(IsQAOn()){
    fPID->SetQAOn();
    fQA->AddAt(fPID->GetQAhistograms(), 8);
  }
  fPID->SetHasMCData(kTRUE);
  fPID->InitializePID("TRD");
}

//____________________________________________________________
void AliAnalysisElectronTask::Exec(Option_t *){
	//
	// Run the analysis
	//
	if(!fESD){
		AliError("No ESD Event");
		return;
	}
	if(!fMC){
		AliError("No MC Event");
		return;
	}
	fCFM->SetEventInfo(fMC);
  fPID->SetMCEvent(fMC);

	//fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);

	Double_t pt = 0;
	Double_t container[3];

	// Loop over the Monte Carlo tracks to see whether we have overlooked any track
	AliMCParticle *mctrack = 0x0;
	Int_t nElectrons = 0;
	for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
		mctrack = fMC->GetTrack(imc);
		container[0] = mctrack->Pt();
		container[1] = mctrack->Eta();
    container[2] = mctrack->Phi();

		if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
		fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
		(dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
		if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
		// find the label in the vector
		fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
		nElectrons++;
	}
	(dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);

	// fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
	AliESDtrack *track = 0x0;
	for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
		track = fESD->GetTrack(itrack);
		container[0] = track->Pt();
		container[1] = track->Eta();
    container[2] = track->Phi();
		if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue;
		fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine);

		// Check TRD criterions (outside the correction framework)
		if(track->GetTRDncls()){
			(dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
			(dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());	// Check the acceptance without tight cuts
			(dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
			(dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
		}
		if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
    if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue;
    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts);
    // track accepted, do PID
		if(!fPID->IsSelected(track)) continue;
  	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1);
	}
	fNEvents->Fill(1);

	// Done!!!
	PostData(0, fNEvents);
	PostData(1, fCFM->GetParticleContainer());
	PostData(2, fQA);
}

//____________________________________________________________
void AliAnalysisElectronTask::Terminate(Option_t *){
	//
	// Terminate not implemented at the moment
	//
}

//____________________________________________________________
void AliAnalysisElectronTask::MakeParticleContainer(){
  //
  // Create the particle container for the correction framework manager and 
  // link it
  //
  const Int_t nvar   = 3 ; //number of variables on the grid:pt,eta, phi
  const Double_t ptmin = 0., ptmax = 10.;
  const Double_t etamin = -0.9, etamax = 0.9;
  const Double_t phimin = 0., phimax = 2. * TMath::Pi();
  

  //arrays for the number of bins in each dimension
  Int_t iBin[nvar];
  iBin[0] = 20; //bins in pt
  iBin[1] =  8; //bins in eta 
  iBin[2] = 18; // bins in phi

  //arrays for lower bounds :
  Double_t *binLim1 = new Double_t[iBin[0] + 1];
  Double_t *binLim2 = new Double_t[iBin[1] + 1];
  Double_t *binLim3 = new Double_t[iBin[2] + 1];

  //values for bin lower bounds
  for(Int_t i=0; i<=iBin[0]; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/iBin[0]*(Double_t)i; 
  for(Int_t i=0; i<=iBin[1]; i++) binLim2[i]=(Double_t)etamin  + (etamax-etamin)/iBin[1]*(Double_t)i;
  for(Int_t i=0; i<=iBin[2]; i++) binLim3[i]=(Double_t)phimin  + (phimax-phimin)/iBin[2]*(Double_t)i;

  //one "container" for MC
  AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, nvar, iBin);
  //setting the bin limits
  container -> SetBinLimits(0,binLim1);
  container -> SetBinLimits(1,binLim2);
  container -> SetBinLimits(2,binLim3);
  fCFM->SetParticleContainer(container);
}
Commit	Line	Data
809a4336	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	**************************************************************************/
	15	/*
	16	* The analysis task:
	17	* Filling an AliCFContainer with the quantities pt, eta and phi
	18	* for tracks which survivied the particle cuts (MC resp. ESD tracks)
	19	* Track selection is done using the AliHFE package
	20	*
	21	* Author:
	22	* Markus Fasel <M.Fasel@gsi.de>
	23	*/
	24	#include <TAxis.h>
	25	#include <TCanvas.h>
	26	#include <TChain.h>
	27	#include <TH1F.h>
	28	#include <TH1I.h>
	29	#include <TH2F.h>
	30	#include <TIterator.h>
	31	#include <TList.h>
	32	#include <TLegend.h>
	33	#include <TMath.h>
	34	#include <TObjArray.h>
	35	#include <TParticle.h>
	36	#include <TProfile.h>
	37	#include <TTree.h>
	38
	39	#include "AliCFContainer.h"
	40	#include "AliCFManager.h"
	41	#include "AliESDEvent.h"
	42	#include "AliESDInputHandler.h"
	43	#include "AliESDtrack.h"
	44	#include "AliESDtrackCuts.h"
	45	#include "AliLog.h"
	46	#include "AliAnalysisManager.h"
	47	#include "AliMCEvent.h"
	48	#include "AliMCEventHandler.h"
	49	#include "AliMCParticle.h"
	50	#include "AliPID.h"
	51
	52	#include "AliHFEpid.h"
	53	#include "AliHFEcuts.h"
	54	#include "AliAnalysisElectronTask.h"
	55
	56	//____________________________________________________________
	57	AliAnalysisElectronTask::AliAnalysisElectronTask():
	58	AliAnalysisTask("PID efficiency Analysis", "")
	59	, fESD(0x0)
	60	, fMC(0x0)
	61	, fCFM(0x0)
	62	, fPID(0x0)
	63	, fCuts(0x0)
	64	, fNEvents(0x0)
65	, fQA(0x0)
66	{
67	DefineInput(0, TChain::Class());
68	DefineOutput(0, TH1I::Class());
69	DefineOutput(1, AliCFContainer::Class());
70	DefineOutput(2, TList::Class());
71
72	// Initialize cuts
73	fCuts = new AliHFEcuts;
74	fPID = new AliHFEpid;
75	}
76
77	//____________________________________________________________
78	AliAnalysisElectronTask::~AliAnalysisElectronTask(){
79	if(fESD) delete fESD;
80	if(fMC) delete fMC;
81	if(fPID) delete fPID;
82	if(fQA) delete fQA;
83	if(fCuts) delete fCuts;
84	if(fNEvents) delete fNEvents;
85	fQA = 0x0; fMC = 0x0; fESD = 0x0; fCFM = 0x0; fPID = 0x0; fCuts = 0x0;
86	}
87
88	//____________________________________________________________
89	void AliAnalysisElectronTask::ConnectInputData(Option_t *){
90	TTree esdchain = dynamic_cast<TChain >(GetInputData(0));
91	if(!esdchain){
92	AliError("ESD chain empty");
93	return;
94	} else {
95	esdchain->SetBranchStatus("Tracks", 1);
96	}
97	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler >(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
98	if(!esdH){
99	AliError("No ESD input handler");
100	return;
101	} else {
102	fESD = esdH->GetEvent();
103	}
104	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler >(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
105	if(!mcH){
106	AliError("No MC truth handler");
107	return;
108	} else {
109	fMC = mcH->MCEvent();
110	}
111	}
112
113	//____________________________________________________________
114	void AliAnalysisElectronTask::CreateOutputObjects(){
115	fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
116	// First Step: TRD alone
117	if(!fQA) fQA = new TList;
118	fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
119	fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
120	fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
121	fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
122	fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
123	fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
124	fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
125
126	// Initialize correction Framework and Cuts
127	fCFM = new AliCFManager;
128	MakeParticleContainer();
129	// Temporary fix: Initialize particle cuts with 0x0
130	for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
131	fCFM->SetParticleCutsList(istep, 0x0);
132	if(IsQAOn()){
133	fCuts->SetDebugMode();
134	fQA->AddAt(fCuts->GetQAhistograms(), 7);
135	}
136	fCuts->CreateStandardCuts();
137	fCuts->Initialize(fCFM);
138
139	// Initialize PID
140	if(IsQAOn()){
141	fPID->SetQAOn();
142	fQA->AddAt(fPID->GetQAhistograms(), 8);
143	}
144	fPID->SetHasMCData(kTRUE);
145	fPID->InitializePID("TRD");
146	}
147
148	//____________________________________________________________
149	void AliAnalysisElectronTask::Exec(Option_t *){
150	//
151	// Run the analysis
152	//
153	if(!fESD){
154	AliError("No ESD Event");
155	return;
156	}
157	if(!fMC){
158	AliError("No MC Event");
159	return;
160	}
161	fCFM->SetEventInfo(fMC);
162	fPID->SetMCEvent(fMC);
163
164	//fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
165
166	Double_t pt = 0;
167	Double_t container[3];
168
169	// Loop over the Monte Carlo tracks to see whether we have overlooked any track
170	AliMCParticle *mctrack = 0x0;
171	Int_t nElectrons = 0;
172	for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
173	mctrack = fMC->GetTrack(imc);
174	container[0] = mctrack->Pt();
175	container[1] = mctrack->Eta();
176	container[2] = mctrack->Phi();
177
178	if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
179	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
180	(dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
181	if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
182	// find the label in the vector
183	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
184	nElectrons++;
185	}
186	(dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
187
188	// fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
189	AliESDtrack *track = 0x0;
190	for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
191	track = fESD->GetTrack(itrack);
192	container[0] = track->Pt();
193	container[1] = track->Eta();
194	container[2] = track->Phi();
195	if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue;
196	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine);
197
198	// Check TRD criterions (outside the correction framework)
199	if(track->GetTRDncls()){
200	(dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
201	(dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts
202	(dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
203	(dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
204	}
205	if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
206	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
207	if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue;
208	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts);
209	// track accepted, do PID
210	if(!fPID->IsSelected(track)) continue;
211	fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1);
212	}
213	fNEvents->Fill(1);
214
215	// Done!!!
216	PostData(0, fNEvents);
217	PostData(1, fCFM->GetParticleContainer());
218	PostData(2, fQA);
219	}
220
221	//____________________________________________________________
222	void AliAnalysisElectronTask::Terminate(Option_t *){
223	//
224	// Terminate not implemented at the moment
225	//
226	}
227
228	//____________________________________________________________
229	void AliAnalysisElectronTask::MakeParticleContainer(){
230	//
231	// Create the particle container for the correction framework manager and
232	// link it
233	//
234	const Int_t nvar = 3 ; //number of variables on the grid:pt,eta, phi
235	const Double_t ptmin = 0., ptmax = 10.;
236	const Double_t etamin = -0.9, etamax = 0.9;
237	const Double_t phimin = 0., phimax = 2. * TMath::Pi();
238
239
240	//arrays for the number of bins in each dimension
241	Int_t iBin[nvar];
242	iBin[0] = 20; //bins in pt
243	iBin[1] = 8; //bins in eta
244	iBin[2] = 18; // bins in phi
245
246	//arrays for lower bounds :
247	Double_t *binLim1 = new Double_t[iBin[0] + 1];
248	Double_t *binLim2 = new Double_t[iBin[1] + 1];
249	Double_t *binLim3 = new Double_t[iBin[2] + 1];
250
251	//values for bin lower bounds
252	for(Int_t i=0; i<=iBin[0]; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/iBin[0]*(Double_t)i;
253	for(Int_t i=0; i<=iBin[1]; i++) binLim2[i]=(Double_t)etamin + (etamax-etamin)/iBin[1]*(Double_t)i;
254	for(Int_t i=0; i<=iBin[2]; i++) binLim3[i]=(Double_t)phimin + (phimax-phimin)/iBin[2]*(Double_t)i;
255
256	//one "container" for MC
257	AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, nvar, iBin);
258	//setting the bin limits
259	container -> SetBinLimits(0,binLim1);
260	container -> SetBinLimits(1,binLim2);
261	container -> SetBinLimits(2,binLim3);
262	fCFM->SetParticleContainer(container);
263	}