Since it contains fixes of coding rule violations, all classes are involved. Further...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskHFE.cxx
CommitLineData
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:
259c3296 22 * Raphaelle Bailhache <R.Bailhache@gsi.de>
809a4336 23 * Markus Fasel <M.Fasel@gsi.de>
259c3296 24 * MinJung Kweon <minjung@physi.uni-heidelberg.de>
809a4336 25 */
26#include <TAxis.h>
27#include <TCanvas.h>
28#include <TChain.h>
259c3296 29#include <TDirectory.h>
809a4336 30#include <TH1F.h>
31#include <TH1I.h>
32#include <TH2F.h>
33#include <TIterator.h>
34#include <TList.h>
35#include <TLegend.h>
36#include <TMath.h>
37#include <TObjArray.h>
38#include <TParticle.h>
39#include <TProfile.h>
75d81601 40#include <TString.h>
809a4336 41#include <TTree.h>
42
43#include "AliCFContainer.h"
44#include "AliCFManager.h"
45#include "AliESDEvent.h"
46#include "AliESDInputHandler.h"
47#include "AliESDtrack.h"
48#include "AliESDtrackCuts.h"
49#include "AliLog.h"
50#include "AliAnalysisManager.h"
51#include "AliMCEvent.h"
52#include "AliMCEventHandler.h"
53#include "AliMCParticle.h"
54#include "AliPID.h"
259c3296 55#include "AliStack.h"
809a4336 56
57#include "AliHFEpid.h"
58#include "AliHFEcuts.h"
259c3296 59#include "AliHFEmcQA.h"
60#include "AliHFEsecVtx.h"
809a4336 61#include "AliAnalysisTaskHFE.h"
62
63//____________________________________________________________
64AliAnalysisTaskHFE::AliAnalysisTaskHFE():
259c3296 65 AliAnalysisTask("PID efficiency Analysis", "")
75d81601 66 , fQAlevel(0)
67 , fPIDdetectors("")
259c3296 68 , fESD(0x0)
69 , fMC(0x0)
70 , fCFM(0x0)
71 , fCorrelation(0x0)
dbe3abbe 72 , fFakeElectrons(0x0)
259c3296 73 , fPID(0x0)
809a4336 74 , fCuts(0x0)
259c3296 75 , fSecVtx(0x0)
dbe3abbe 76 , fMCQA(0x0)
259c3296 77 , fNEvents(0x0)
75d81601 78 , fNElectronTracksEvent(0x0)
259c3296 79 , fQA(0x0)
80 , fOutput(0x0)
81 , fHistMCQA(0x0)
82 , fHistSECVTX(0x0)
809a4336 83{
dbe3abbe 84 //
85 // Default constructor
86 //
87 DefineInput(0, TChain::Class());
88 DefineOutput(0, TH1I::Class());
89 DefineOutput(1, TList::Class());
259c3296 90 DefineOutput(2, TList::Class());
809a4336 91
92 // Initialize cuts
93 fCuts = new AliHFEcuts;
94 fPID = new AliHFEpid;
95}
96
dbe3abbe 97//____________________________________________________________
98AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
75d81601 99 AliAnalysisTask(ref)
100 , fQAlevel(ref.fQAlevel)
101 , fPIDdetectors(ref.fPIDdetectors)
102 , fESD(ref.fESD)
103 , fMC(ref.fMC)
104 , fCFM(ref.fCFM)
105 , fCorrelation(ref.fCorrelation)
106 , fFakeElectrons(ref.fFakeElectrons)
107 , fPID(ref.fPID)
108 , fCuts(ref.fCuts)
109 , fSecVtx(ref.fSecVtx)
110 , fMCQA(ref.fMCQA)
111 , fNEvents(ref.fNEvents)
112 , fNElectronTracksEvent(ref.fNElectronTracksEvent)
113 , fQA(ref.fQA)
114 , fOutput(ref.fOutput)
115 , fHistMCQA(ref.fHistMCQA)
116 , fHistSECVTX(ref.fHistSECVTX)
dbe3abbe 117{
118 //
119 // Copy Constructor
120 //
dbe3abbe 121}
122
123//____________________________________________________________
124AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
125 //
126 // Assignment operator
127 //
75d81601 128 if(this == &ref) return *this;
129 AliAnalysisTask::operator=(ref);
130 fQAlevel = ref.fQAlevel;
131 fPIDdetectors = ref.fPIDdetectors;
132 fESD = ref.fESD;
133 fMC = ref.fMC;
134 fCFM = ref.fCFM;
135 fCorrelation = ref.fCorrelation;
136 fFakeElectrons = ref.fFakeElectrons;
137 fPID = ref.fPID;
138 fCuts = ref.fCuts;
139 fSecVtx = ref.fSecVtx;
140 fMCQA = ref.fMCQA;
141 fNEvents = ref.fNEvents;
142 fNElectronTracksEvent = ref.fNElectronTracksEvent;
143 fQA = ref.fQA;
144 fOutput = ref.fOutput;
145 fHistMCQA = ref.fHistMCQA;
146 fHistSECVTX = ref.fHistSECVTX;
dbe3abbe 147 return *this;
148}
149
809a4336 150//____________________________________________________________
151AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
dbe3abbe 152 //
153 // Destructor
154 //
155 if(fESD) delete fESD;
156 if(fMC) delete fMC;
157 if(fPID) delete fPID;
158 if(fQA){
259c3296 159 fQA->Clear();
160 delete fQA;
161 }
162 if(fOutput){
163 fOutput->Clear();
164 delete fOutput;
165 }
166 if(fHistMCQA){
167 fHistMCQA->Clear();
168 delete fHistMCQA;
169 }
170 if(fHistSECVTX){
171 fHistSECVTX->Clear();
172 delete fHistSECVTX;
173 }
809a4336 174 if(fCuts) delete fCuts;
259c3296 175 if(fSecVtx) delete fSecVtx;
dbe3abbe 176 if(fMCQA) delete fMCQA;
177 if(fNEvents) delete fNEvents;
259c3296 178 if(fCorrelation) delete fCorrelation;
179 if(fFakeElectrons) delete fFakeElectrons;
809a4336 180}
181
182//____________________________________________________________
183void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
75d81601 184/* TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
dbe3abbe 185 if(!esdchain){
186 AliError("ESD chain empty");
187 return;
188 } else {
189 esdchain->SetBranchStatus("Tracks", 1);
190 }
75d81601 191*/
dbe3abbe 192 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
193 if(!esdH){
194 AliError("No ESD input handler");
195 return;
196 } else {
197 fESD = esdH->GetEvent();
198 }
199 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
200 if(!mcH){
201 AliError("No MC truth handler");
202 return;
203 } else {
204 fMC = mcH->MCEvent();
205 }
809a4336 206}
207
208//____________________________________________________________
209void AliAnalysisTaskHFE::CreateOutputObjects(){
dbe3abbe 210 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
75d81601 211 fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
dbe3abbe 212 // First Step: TRD alone
213 if(!fQA) fQA = new TList;
214 fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
215 fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
216 fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
217 fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
218 fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
219 fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
220 fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
809a4336 221
259c3296 222 if(!fOutput) fOutput = new TList;
809a4336 223 // Initialize correction Framework and Cuts
224 fCFM = new AliCFManager;
225 MakeParticleContainer();
226 // Temporary fix: Initialize particle cuts with 0x0
227 for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
228 fCFM->SetParticleCutsList(istep, 0x0);
75d81601 229 if(IsQAOn(kCUTqa)){
230 AliInfo("QA on for Cuts");
809a4336 231 fCuts->SetDebugMode();
259c3296 232 fQA->Add(fCuts->GetQAhistograms());
809a4336 233 }
234 fCuts->CreateStandardCuts();
75d81601 235 fCuts->SetMinNTrackletsTRD(0); // Minimal requirement to get a minimum biased electron sample (only TPC and ITS requirements allowed)
809a4336 236 fCuts->Initialize(fCFM);
259c3296 237 // add output objects to the List
238 fOutput->AddAt(fCFM->GetParticleContainer(), 0);
239 fOutput->AddAt(fCorrelation, 1);
240 fOutput->AddAt(fFakeElectrons, 2);
75d81601 241 fOutput->AddAt(fNElectronTracksEvent, 3);
809a4336 242
243 // Initialize PID
75d81601 244 if(IsQAOn(kPIDqa)){
245 AliInfo("PID QA switched on");
809a4336 246 fPID->SetQAOn();
259c3296 247 fQA->Add(fPID->GetQAhistograms());
809a4336 248 }
249 fPID->SetHasMCData(kTRUE);
75d81601 250 if(!fPIDdetectors.Length()) AddPIDdetector("TPC");
251 fPID->InitializePID(fPIDdetectors.Data()); // Only restrictions to TPC allowed
259c3296 252
dbe3abbe 253 // mcQA----------------------------------
75d81601 254 if (IsQAOn(kMCqa)) {
dbe3abbe 255 AliInfo("MC QA on");
256 if(!fMCQA) fMCQA = new AliHFEmcQA;
259c3296 257 if(!fHistMCQA) fHistMCQA = new TList();
75d81601 258 fHistMCQA->SetName("MCqa");
dbe3abbe 259 fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm
260 fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty
261 fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm
262 fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty
75d81601 263 TIter next(gDirectory->GetList());
264 TObject *obj;
265 int counter = 0;
dbe3abbe 266 TString objname;
75d81601 267 while ((obj = next.Next())) {
268 objname = obj->GetName();
dbe3abbe 269 TObjArray *toks = objname.Tokenize("_");
270 if (toks->GetEntriesFast()){
271 TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
75d81601 272 if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
dbe3abbe 273 }
274 }
259c3296 275 fQA->Add(fHistMCQA);
dbe3abbe 276 }
259c3296 277
dbe3abbe 278 // secvtx----------------------------------
259c3296 279 if (IsSecVtxOn()) {
dbe3abbe 280 AliInfo("Secondary Vertex Analysis on");
281 fSecVtx = new AliHFEsecVtx;
282
259c3296 283 if(!fHistSECVTX) fHistSECVTX = new TList();
75d81601 284 fHistSECVTX->SetName("SecVtx");
259c3296 285 fSecVtx->CreateHistograms("secvtx_");
75d81601 286 TIter next(gDirectory->GetList());
287 TObject *obj;
288 int counter = 0;
dbe3abbe 289 TString objname;
75d81601 290 while ((obj = next.Next())) {
291 objname = obj->GetName();
dbe3abbe 292 TObjArray *toks = objname.Tokenize("_");
293 if (toks->GetEntriesFast()){
294 TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
75d81601 295 if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
259c3296 296 }
297 }
75d81601 298 fOutput->Add(fHistSECVTX);
dbe3abbe 299 }
809a4336 300}
301
302//____________________________________________________________
303void AliAnalysisTaskHFE::Exec(Option_t *){
dbe3abbe 304 //
305 // Run the analysis
306 //
307 if(!fESD){
308 AliError("No ESD Event");
309 return;
310 }
311 if(!fMC){
312 AliError("No MC Event");
313 return;
314 }
315 fCFM->SetEventInfo(fMC);
809a4336 316 fPID->SetMCEvent(fMC);
317
dbe3abbe 318 //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
319 Double_t container[6];
809a4336 320
dbe3abbe 321 // Loop over the Monte Carlo tracks to see whether we have overlooked any track
322 AliMCParticle *mctrack = 0x0;
323 Int_t nElectrons = 0;
809a4336 324
dbe3abbe 325 if (IsSecVtxOn()) {
326 fSecVtx->SetEvent(fESD);
327 fSecVtx->SetStack(fMC->Stack());
328 }
259c3296 329
dbe3abbe 330 // run mc QA ------------------------------------------------
75d81601 331 if (IsQAOn(kMCqa)) {
dbe3abbe 332 AliDebug(2, "Running MC QA");
259c3296 333
dbe3abbe 334 fMCQA->SetStack(fMC->Stack());
335 fMCQA->Init();
259c3296 336
337 Int_t nPrims = fMC->Stack()->GetNprimary();
338 Int_t nMCTracks = fMC->Stack()->GetNtrack();
339
340 // loop over primary particles for quark and heavy hadrons
341 for (Int_t igen = 0; igen < nPrims; igen++){
dbe3abbe 342 fMCQA->GetQuarkKine(igen, AliHFEmcQA::kCharm);
343 fMCQA->GetQuarkKine(igen, AliHFEmcQA::kBeauty);
344 fMCQA->GetHadronKine(igen, AliHFEmcQA::kCharm);
345 fMCQA->GetHadronKine(igen, AliHFEmcQA::kBeauty);
259c3296 346 }
dbe3abbe 347 fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
348 fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
259c3296 349
350 // loop over all tracks for decayed electrons
351 for (Int_t igen = 0; igen < nMCTracks; igen++){
dbe3abbe 352 fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0);
353 fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0);
354 fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
355 fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
259c3296 356 }
357
358 } // end of MC QA loop
359 // -----------------------------------------------------------------
360
dbe3abbe 361 //
362 // Loop MC
363 //
364 for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
365 mctrack = fMC->GetTrack(imc);
366
367 container[0] = mctrack->Pt();
368 container[1] = mctrack->Eta();
809a4336 369 container[2] = mctrack->Phi();
370
dbe3abbe 371 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
372 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
373 (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
374 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
375 // find the label in the vector
376 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
377 nElectrons++;
378 }
379 (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
380
381 // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
382
383
384 //
385 // Loop ESD
386 //
387
388 Int_t nbtrackrec = fESD->GetNumberOfTracks();
75d81601 389 Int_t nElectronCandidates = 0;
dbe3abbe 390 AliESDtrack *track = 0x0, *htrack = 0x0;
259c3296 391 Int_t pid = 0;
dbe3abbe 392 // For double counted tracks
393 Int_t *indexRecKineTPC = new Int_t[nbtrackrec];
394 Int_t *indexRecKineITS = new Int_t[nbtrackrec];
395 Int_t *indexRecPrim = new Int_t[nbtrackrec];
396 Int_t *indexHFEcutsITS = new Int_t[nbtrackrec];
397 Int_t *indexHFEcutsTPC = new Int_t[nbtrackrec];
398 Int_t *indexHFEcutsTRD = new Int_t[nbtrackrec];
399 Int_t *indexPid = new Int_t[nbtrackrec];
400
401 for(Int_t nb = 0; nb < nbtrackrec; nb++){
402 indexRecKineTPC[nb] = -1;
403 indexRecKineITS[nb] = -1;
404 indexRecPrim[nb] = -1;
405 indexHFEcutsITS[nb] = -1;
406 indexHFEcutsTPC[nb] = -1;
407 indexHFEcutsTRD[nb] = -1;
408 indexPid[nb] = -1;
409 }
410
411 Bool_t alreadyseen = kFALSE;
412
413 for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
414
415 track = fESD->GetTrack(itrack);
416
417 container[0] = track->Pt();
418 container[1] = track->Eta();
809a4336 419 container[2] = track->Phi();
dbe3abbe 420
421 Bool_t signal = kTRUE;
422
423 // RecKine: TPC cuts
424 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineTPC, track)) continue;
425 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineTPC);
426
427
428 // Check if it is signal electrons
429 if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue;
430
431 container[3] = mctrack->Pt();
432 container[4] = mctrack->Eta();
433 container[5] = mctrack->Phi();
434
435 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
436
437 if(signal) {
438 alreadyseen = kFALSE;
439 for(Int_t nb = 0; nb < itrack; nb++){
440 if(indexRecKineTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
441 }
442 indexRecKineTPC[itrack] = TMath::Abs(track->GetLabel());
443
444 fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineTPC + AliHFEcuts::kNcutESDSteps + 1));
445 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)));
446 if(alreadyseen) {
447 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
448 }
449 else {
450 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
451 }
452 }
453
454
455 // RecKine: ITS cuts
456 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITS, track)) continue;
457 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineITS);
458 if(signal) {
459 alreadyseen = kFALSE;
460 for(Int_t nb = 0; nb < itrack; nb++){
461 if(indexRecKineITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
462 }
463 indexRecKineITS[itrack] = TMath::Abs(track->GetLabel());
464 fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineITS + AliHFEcuts::kNcutESDSteps + 1));
465 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 2*(AliHFEcuts::kNcutESDSteps + 1)));
466 if(alreadyseen) {
467 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
468 }
469 else {
470 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
471 }
472 }
473
474
475 // Check TRD criterions (outside the correction framework)
476 if(track->GetTRDncls()){
477 (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
478 (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts
479 (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
480 (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
481 }
482
483
484 // RecPrim
485 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
809a4336 486 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
dbe3abbe 487 if(signal) {
488 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
489 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
490 alreadyseen = kFALSE;
491 for(Int_t nb = 0; nb < itrack; nb++){
492 if(indexRecPrim[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
493 }
494 indexRecPrim[itrack] = TMath::Abs(track->GetLabel());
495 if(alreadyseen) {
496 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 4*(AliHFEcuts::kNcutESDSteps + 1));
497 }
498 else {
499 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 3*(AliHFEcuts::kNcutESDSteps + 1));
500 }
501 }
502
503
504 // HFEcuts: ITS layers cuts
505 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
506 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS);
507 if(signal) {
508 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
509 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
510 alreadyseen = kFALSE;
511 for(Int_t nb = 0; nb < itrack; nb++){
512 if(indexHFEcutsITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
513 }
514 indexHFEcutsITS[itrack] = TMath::Abs(track->GetLabel());
515 if(alreadyseen) {
516 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
517 }
518 else {
519 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
520 }
521 }
522
523 // HFEcuts: TPC ratio clusters
524 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
525 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC);
526 if(signal) {
527 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutESDSteps + 1);
528 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTPC + 2*(AliHFEcuts::kNcutESDSteps + 1));
529 alreadyseen = kFALSE;
530 for(Int_t nb = 0; nb < itrack; nb++){
531 if(indexHFEcutsTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
532 }
533 indexHFEcutsTPC[itrack] = TMath::Abs(track->GetLabel());
534 if(alreadyseen) {
535 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
536 }
537 else {
538 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
539 }
540 }
541
542 // HFEcuts: Nb of tracklets TRD
543 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
544 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD);
545 if(signal) {
546 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 1);
547 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
548 alreadyseen = kFALSE;
549 for(Int_t nb = 0; nb < itrack; nb++){
550 if(indexHFEcutsTRD[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
551 }
552 indexHFEcutsTRD[itrack] = TMath::Abs(track->GetLabel());
553 if(alreadyseen) {
554 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 4*(AliHFEcuts::kNcutESDSteps + 1)));
555 }
556 else {
557 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 3*(AliHFEcuts::kNcutESDSteps + 1)));
558 }
559 }
560
561
809a4336 562 // track accepted, do PID
dbe3abbe 563 if(!fPID->IsSelected(track)) continue;
75d81601 564 nElectronCandidates++;
dbe3abbe 565
566
567 // Fill Containers
568 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 1);
569 if(signal) {
570 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 2);
571 fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)+1);
572 alreadyseen = kFALSE;
573 for(Int_t nb = 0; nb < itrack; nb++){
574 if(indexPid[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
575 }
576 indexPid[itrack] = TMath::Abs(track->GetLabel());
577 if(alreadyseen) {
578 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)));
579 }
580 else {
581 fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 3*(AliHFEcuts::kNcutESDSteps + 1)));
582 }
259c3296 583 // dimensions 3&4&5 : pt,eta,phi (MC)
259c3296 584 fCorrelation->Fill(container);
dbe3abbe 585 }
586
587
588 // Track selected: distinguish between true and fake
589 if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
259c3296 590 // pair analysis [mj]
591 if (IsSecVtxOn()) {
dbe3abbe 592 AliDebug(2, "Running Secondary Vertex Analysis");
593 fSecVtx->InitAnaPair();
259c3296 594 for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
595 htrack = fESD->GetTrack(jtrack);
596 if ( itrack == jtrack ) continue;
597 //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
598 if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
dbe3abbe 599 fSecVtx->AnaPair(track, htrack, jtrack);
259c3296 600 }
dbe3abbe 601 // based on the partner of e info, you run secandary vertexing function
602 fSecVtx->RunSECVTX(track);
603 } // end of pair analysis
259c3296 604 } else {
605 // Fill THnSparse with the information for Fake Electrons
606 switch(pid){
dbe3abbe 607 case 13: container[3] = AliPID::kMuon;
608 case 211: container[3] = AliPID::kPion;
609 case 321: container[3] = AliPID::kKaon;
610 case 2212: container[3] = AliPID::kProton;
259c3296 611 }
612 fFakeElectrons->Fill(container);
613 }
dbe3abbe 614 }
615 fNEvents->Fill(1);
75d81601 616 fNElectronTracksEvent->Fill(nElectronCandidates);
617
dbe3abbe 618 // release memory
619 delete[] indexRecKineTPC;
620 delete[] indexRecKineITS;
621 delete[] indexRecPrim;
622 delete[] indexHFEcutsITS;
623 delete[] indexHFEcutsTPC;
624 delete[] indexHFEcutsTRD;
625 delete[] indexPid;
626
627 // Done!!!
628 PostData(0, fNEvents);
629 PostData(1, fOutput);
630 PostData(2, fQA);
809a4336 631}
632
633//____________________________________________________________
634void AliAnalysisTaskHFE::Terminate(Option_t *){
dbe3abbe 635 //
636 // Terminate not implemented at the moment
637 //
809a4336 638}
639
640//____________________________________________________________
641void AliAnalysisTaskHFE::MakeParticleContainer(){
642 //
643 // Create the particle container for the correction framework manager and
644 // link it
645 //
259c3296 646 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
dbe3abbe 647 const Double_t kPtmin = 0.1, kPtmax = 10.;
259c3296 648 const Double_t kEtamin = -0.9, kEtamax = 0.9;
649 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
809a4336 650
651 //arrays for the number of bins in each dimension
259c3296 652 Int_t iBin[kNvar];
dbe3abbe 653 iBin[0] = 40; //bins in pt
809a4336 654 iBin[1] = 8; //bins in eta
655 iBin[2] = 18; // bins in phi
656
657 //arrays for lower bounds :
259c3296 658 Double_t* binEdges[kNvar];
659 for(Int_t ivar = 0; ivar < kNvar; ivar++)
660 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
809a4336 661
662 //values for bin lower bounds
dbe3abbe 663 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
259c3296 664 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
665 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
809a4336 666
667 //one "container" for MC
dbe3abbe 668 AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin);
669
809a4336 670 //setting the bin limits
dbe3abbe 671 for(Int_t ivar = 0; ivar < kNvar; ivar++)
672 container -> SetBinLimits(ivar, binEdges[ivar]);
809a4336 673 fCFM->SetParticleContainer(container);
259c3296 674
675 //create correlation matrix for unfolding
676 Int_t thnDim[2*kNvar];
677 for (int k=0; k<kNvar; k++) {
678 //first half : reconstructed
679 //second half : MC
680 thnDim[k] = iBin[k];
681 thnDim[k+kNvar] = iBin[k];
682 }
683
684 fCorrelation = new THnSparseF("correlation","THnSparse with correlations",2*kNvar,thnDim);
685 for (int k=0; k<kNvar; k++) {
686 fCorrelation->SetBinEdges(k,binEdges[k]);
687 fCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
688 }
689 fCorrelation->Sumw2();
690
691 // Add a histogram for Fake electrons
692 thnDim[3] = AliPID::kSPECIES;
693 fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
694 for(Int_t idim = 0; idim < kNvar; idim++)
695 fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
809a4336 696}
dbe3abbe 697
75d81601 698void AliAnalysisTaskHFE::AddPIDdetector(Char_t *detector){
699 if(!fPIDdetectors.Length())
700 fPIDdetectors = detector;
701 else
702 fPIDdetectors += Form(":%s", detector);
703}
704
705//____________________________________________________________
706void AliAnalysisTaskHFE::PrintStatus(){
707 printf("\n\tAnalysis Settings\n\t========================================\n\n");
708 printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
709 printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
710 printf("\n");
711 printf("\tParticle Identification Detectors:\n");
712 TObjArray *detectors = fPIDdetectors.Tokenize(":");
713 for(Int_t idet = 0; idet < detectors->GetEntries(); idet++)
714 printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data());
715 printf("\n");
716 printf("\tQA: \n");
717 printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" : "NO");
718 printf("\t\tCUTS: %s\n", IsQAOn(kCUTqa) ? "YES" : "NO");
719 printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
720 printf("\n");
721}