]>
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 | **************************************************************************/ | |
50685501 | 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 | // Raphaelle Bailhache <R.Bailhache@gsi.de> | |
23 | // Markus Fasel <M.Fasel@gsi.de> | |
24 | // Matus Kalisky <matus.kalisky@cern.ch> | |
25 | // MinJung Kweon <minjung@physi.uni-heidelberg.de> | |
26 | // | |
809a4336 | 27 | #include <TAxis.h> |
28 | #include <TCanvas.h> | |
29 | #include <TChain.h> | |
259c3296 | 30 | #include <TDirectory.h> |
6ad05e72 | 31 | #include <TFile.h> |
32 | #include <TH1D.h> | |
809a4336 | 33 | #include <TH1F.h> |
34 | #include <TH1I.h> | |
35 | #include <TH2F.h> | |
36 | #include <TIterator.h> | |
37 | #include <TList.h> | |
38 | #include <TLegend.h> | |
39 | #include <TMath.h> | |
40 | #include <TObjArray.h> | |
41 | #include <TParticle.h> | |
42 | #include <TProfile.h> | |
75d81601 | 43 | #include <TString.h> |
809a4336 | 44 | #include <TTree.h> |
45 | ||
46 | #include "AliCFContainer.h" | |
47 | #include "AliCFManager.h" | |
722347d8 | 48 | #include "AliCFEffGrid.h" |
809a4336 | 49 | #include "AliESDEvent.h" |
50 | #include "AliESDInputHandler.h" | |
51 | #include "AliESDtrack.h" | |
52 | #include "AliESDtrackCuts.h" | |
53 | #include "AliLog.h" | |
54 | #include "AliAnalysisManager.h" | |
55 | #include "AliMCEvent.h" | |
56 | #include "AliMCEventHandler.h" | |
57 | #include "AliMCParticle.h" | |
58 | #include "AliPID.h" | |
259c3296 | 59 | #include "AliStack.h" |
809a4336 | 60 | |
61 | #include "AliHFEpid.h" | |
62 | #include "AliHFEcuts.h" | |
259c3296 | 63 | #include "AliHFEmcQA.h" |
64 | #include "AliHFEsecVtx.h" | |
809a4336 | 65 | #include "AliAnalysisTaskHFE.h" |
66 | ||
67 | //____________________________________________________________ | |
68 | AliAnalysisTaskHFE::AliAnalysisTaskHFE(): | |
259c3296 | 69 | AliAnalysisTask("PID efficiency Analysis", "") |
75d81601 | 70 | , fQAlevel(0) |
71 | , fPIDdetectors("") | |
0792aa82 | 72 | , fPIDstrategy(0) |
73 | , fESD(0x0) | |
74 | , fMC(0x0) | |
75 | , fCFM(0x0) | |
76 | , fCorrelation(0x0) | |
77 | , fPIDperformance(0x0) | |
78 | , fPID(0x0) | |
79 | , fCuts(0x0) | |
80 | , fSecVtx(0x0) | |
81 | , fMCQA(0x0) | |
82 | , fNEvents(0x0) | |
83 | , fNElectronTracksEvent(0x0) | |
84 | , fQA(0x0) | |
85 | , fOutput(0x0) | |
86 | , fHistMCQA(0x0) | |
87 | , fHistSECVTX(0x0) | |
88 | { | |
89 | // | |
50685501 | 90 | // Dummy constructor |
0792aa82 | 91 | // |
92 | DefineInput(0, TChain::Class()); | |
93 | DefineOutput(0, TH1I::Class()); | |
94 | DefineOutput(1, TList::Class()); | |
95 | DefineOutput(2, TList::Class()); | |
96 | ||
97 | // Initialize cuts | |
98 | fPID = new AliHFEpid; | |
99 | ||
100 | SetHasMCData(); | |
101 | } | |
102 | ||
103 | //____________________________________________________________ | |
104 | AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name): | |
105 | AliAnalysisTask(name, "") | |
106 | , fQAlevel(0) | |
107 | , fPIDdetectors("") | |
108 | , fPIDstrategy(0) | |
259c3296 | 109 | , fESD(0x0) |
110 | , fMC(0x0) | |
111 | , fCFM(0x0) | |
112 | , fCorrelation(0x0) | |
78ea5ef4 | 113 | , fPIDperformance(0x0) |
259c3296 | 114 | , fPID(0x0) |
809a4336 | 115 | , fCuts(0x0) |
259c3296 | 116 | , fSecVtx(0x0) |
dbe3abbe | 117 | , fMCQA(0x0) |
259c3296 | 118 | , fNEvents(0x0) |
75d81601 | 119 | , fNElectronTracksEvent(0x0) |
259c3296 | 120 | , fQA(0x0) |
121 | , fOutput(0x0) | |
122 | , fHistMCQA(0x0) | |
123 | , fHistSECVTX(0x0) | |
809a4336 | 124 | { |
dbe3abbe | 125 | // |
126 | // Default constructor | |
127 | // | |
128 | DefineInput(0, TChain::Class()); | |
129 | DefineOutput(0, TH1I::Class()); | |
130 | DefineOutput(1, TList::Class()); | |
259c3296 | 131 | DefineOutput(2, TList::Class()); |
809a4336 | 132 | |
133 | // Initialize cuts | |
809a4336 | 134 | fPID = new AliHFEpid; |
722347d8 | 135 | |
136 | SetHasMCData(); | |
809a4336 | 137 | } |
138 | ||
dbe3abbe | 139 | //____________________________________________________________ |
140 | AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref): | |
75d81601 | 141 | AliAnalysisTask(ref) |
142 | , fQAlevel(ref.fQAlevel) | |
143 | , fPIDdetectors(ref.fPIDdetectors) | |
0792aa82 | 144 | , fPIDstrategy(ref.fPIDstrategy) |
75d81601 | 145 | , fESD(ref.fESD) |
146 | , fMC(ref.fMC) | |
147 | , fCFM(ref.fCFM) | |
148 | , fCorrelation(ref.fCorrelation) | |
78ea5ef4 | 149 | , fPIDperformance(ref.fPIDperformance) |
75d81601 | 150 | , fPID(ref.fPID) |
151 | , fCuts(ref.fCuts) | |
152 | , fSecVtx(ref.fSecVtx) | |
153 | , fMCQA(ref.fMCQA) | |
154 | , fNEvents(ref.fNEvents) | |
155 | , fNElectronTracksEvent(ref.fNElectronTracksEvent) | |
156 | , fQA(ref.fQA) | |
157 | , fOutput(ref.fOutput) | |
158 | , fHistMCQA(ref.fHistMCQA) | |
159 | , fHistSECVTX(ref.fHistSECVTX) | |
dbe3abbe | 160 | { |
161 | // | |
162 | // Copy Constructor | |
163 | // | |
dbe3abbe | 164 | } |
165 | ||
166 | //____________________________________________________________ | |
167 | AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){ | |
168 | // | |
169 | // Assignment operator | |
170 | // | |
75d81601 | 171 | if(this == &ref) return *this; |
172 | AliAnalysisTask::operator=(ref); | |
173 | fQAlevel = ref.fQAlevel; | |
174 | fPIDdetectors = ref.fPIDdetectors; | |
0792aa82 | 175 | fPIDstrategy = ref.fPIDstrategy; |
75d81601 | 176 | fESD = ref.fESD; |
177 | fMC = ref.fMC; | |
178 | fCFM = ref.fCFM; | |
179 | fCorrelation = ref.fCorrelation; | |
78ea5ef4 | 180 | fPIDperformance = ref.fPIDperformance; |
75d81601 | 181 | fPID = ref.fPID; |
182 | fCuts = ref.fCuts; | |
183 | fSecVtx = ref.fSecVtx; | |
184 | fMCQA = ref.fMCQA; | |
185 | fNEvents = ref.fNEvents; | |
186 | fNElectronTracksEvent = ref.fNElectronTracksEvent; | |
187 | fQA = ref.fQA; | |
188 | fOutput = ref.fOutput; | |
189 | fHistMCQA = ref.fHistMCQA; | |
190 | fHistSECVTX = ref.fHistSECVTX; | |
dbe3abbe | 191 | return *this; |
192 | } | |
193 | ||
809a4336 | 194 | //____________________________________________________________ |
195 | AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){ | |
dbe3abbe | 196 | // |
197 | // Destructor | |
198 | // | |
199 | if(fESD) delete fESD; | |
200 | if(fMC) delete fMC; | |
201 | if(fPID) delete fPID; | |
202 | if(fQA){ | |
259c3296 | 203 | fQA->Clear(); |
204 | delete fQA; | |
205 | } | |
206 | if(fOutput){ | |
207 | fOutput->Clear(); | |
208 | delete fOutput; | |
209 | } | |
210 | if(fHistMCQA){ | |
211 | fHistMCQA->Clear(); | |
212 | delete fHistMCQA; | |
213 | } | |
214 | if(fHistSECVTX){ | |
215 | fHistSECVTX->Clear(); | |
216 | delete fHistSECVTX; | |
217 | } | |
259c3296 | 218 | if(fSecVtx) delete fSecVtx; |
dbe3abbe | 219 | if(fMCQA) delete fMCQA; |
220 | if(fNEvents) delete fNEvents; | |
78ea5ef4 | 221 | if(fCorrelation){ |
222 | fCorrelation->Clear(); | |
223 | delete fCorrelation; | |
224 | } | |
225 | if(fPIDperformance) delete fPIDperformance; | |
809a4336 | 226 | } |
227 | ||
228 | //____________________________________________________________ | |
229 | void AliAnalysisTaskHFE::ConnectInputData(Option_t *){ | |
50685501 | 230 | // |
231 | // Connecting the input | |
232 | // Called once per worker | |
233 | // | |
75d81601 | 234 | /* TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0)); |
dbe3abbe | 235 | if(!esdchain){ |
236 | AliError("ESD chain empty"); | |
237 | return; | |
238 | } else { | |
239 | esdchain->SetBranchStatus("Tracks", 1); | |
240 | } | |
75d81601 | 241 | */ |
dbe3abbe | 242 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
243 | if(!esdH){ | |
244 | AliError("No ESD input handler"); | |
245 | return; | |
246 | } else { | |
247 | fESD = esdH->GetEvent(); | |
248 | } | |
249 | AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
250 | if(!mcH){ | |
251 | AliError("No MC truth handler"); | |
252 | return; | |
253 | } else { | |
254 | fMC = mcH->MCEvent(); | |
255 | } | |
809a4336 | 256 | } |
257 | ||
258 | //____________________________________________________________ | |
259 | void AliAnalysisTaskHFE::CreateOutputObjects(){ | |
50685501 | 260 | // |
261 | // Creating output container and output objects | |
262 | // Here we also Initialize the correction framework container and | |
263 | // the objects for | |
264 | // - PID | |
265 | // - MC QA | |
266 | // - SecVtx | |
267 | // QA histograms are created if requested | |
268 | // Called once per worker | |
269 | // | |
dbe3abbe | 270 | 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 | 271 | fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100); |
dbe3abbe | 272 | // First Step: TRD alone |
273 | if(!fQA) fQA = new TList; | |
274 | fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0); | |
275 | fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1); | |
276 | fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2); | |
277 | fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3); | |
278 | fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4); | |
279 | fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5); | |
280 | fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6); | |
809a4336 | 281 | |
259c3296 | 282 | if(!fOutput) fOutput = new TList; |
809a4336 | 283 | // Initialize correction Framework and Cuts |
284 | fCFM = new AliCFManager; | |
285 | MakeParticleContainer(); | |
286 | // Temporary fix: Initialize particle cuts with 0x0 | |
287 | for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++) | |
288 | fCFM->SetParticleCutsList(istep, 0x0); | |
0792aa82 | 289 | if(!fCuts){ |
290 | AliWarning("Cuts not available. Default cuts will be used"); | |
291 | fCuts = new AliHFEcuts; | |
292 | fCuts->CreateStandardCuts(); | |
809a4336 | 293 | } |
0792aa82 | 294 | fCuts->Initialize(fCFM); |
295 | if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms()); | |
296 | ||
259c3296 | 297 | // add output objects to the List |
298 | fOutput->AddAt(fCFM->GetParticleContainer(), 0); | |
299 | fOutput->AddAt(fCorrelation, 1); | |
78ea5ef4 | 300 | fOutput->AddAt(fPIDperformance, 2); |
75d81601 | 301 | fOutput->AddAt(fNElectronTracksEvent, 3); |
809a4336 | 302 | |
303 | // Initialize PID | |
75d81601 | 304 | if(IsQAOn(kPIDqa)){ |
305 | AliInfo("PID QA switched on"); | |
722347d8 | 306 | //fPID->SetDebugLevel(2); |
809a4336 | 307 | fPID->SetQAOn(); |
259c3296 | 308 | fQA->Add(fPID->GetQAhistograms()); |
809a4336 | 309 | } |
722347d8 | 310 | fPID->SetHasMCData(HasMCData()); |
0792aa82 | 311 | if(!fPIDdetectors.Length() && ! fPIDstrategy) AddPIDdetector("TPC"); |
312 | if(fPIDstrategy) | |
313 | fPID->InitializePID(Form("Strategy%d", fPIDstrategy)); | |
314 | else | |
315 | fPID->InitializePID(fPIDdetectors.Data()); // Only restrictions to TPC allowed | |
259c3296 | 316 | |
dbe3abbe | 317 | // mcQA---------------------------------- |
75d81601 | 318 | if (IsQAOn(kMCqa)) { |
dbe3abbe | 319 | AliInfo("MC QA on"); |
320 | if(!fMCQA) fMCQA = new AliHFEmcQA; | |
259c3296 | 321 | if(!fHistMCQA) fHistMCQA = new TList(); |
75d81601 | 322 | fHistMCQA->SetName("MCqa"); |
dbe3abbe | 323 | fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm |
324 | fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty | |
325 | fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm | |
326 | fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty | |
0792aa82 | 327 | fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm |
328 | fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty | |
329 | fMCQA->CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_"); // create histograms for charm | |
330 | fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_"); // create histograms for beauty | |
331 | fMCQA->CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_"); // create histograms for charm | |
332 | fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_"); // create histograms for beauty | |
75d81601 | 333 | TIter next(gDirectory->GetList()); |
334 | TObject *obj; | |
335 | int counter = 0; | |
dbe3abbe | 336 | TString objname; |
75d81601 | 337 | while ((obj = next.Next())) { |
338 | objname = obj->GetName(); | |
dbe3abbe | 339 | TObjArray *toks = objname.Tokenize("_"); |
340 | if (toks->GetEntriesFast()){ | |
341 | TObjString *fpart = (TObjString *)(toks->UncheckedAt(0)); | |
75d81601 | 342 | if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++); |
dbe3abbe | 343 | } |
344 | } | |
259c3296 | 345 | fQA->Add(fHistMCQA); |
dbe3abbe | 346 | } |
259c3296 | 347 | |
dbe3abbe | 348 | // secvtx---------------------------------- |
259c3296 | 349 | if (IsSecVtxOn()) { |
dbe3abbe | 350 | AliInfo("Secondary Vertex Analysis on"); |
351 | fSecVtx = new AliHFEsecVtx; | |
352 | ||
259c3296 | 353 | if(!fHistSECVTX) fHistSECVTX = new TList(); |
75d81601 | 354 | fHistSECVTX->SetName("SecVtx"); |
259c3296 | 355 | fSecVtx->CreateHistograms("secvtx_"); |
75d81601 | 356 | TIter next(gDirectory->GetList()); |
357 | TObject *obj; | |
358 | int counter = 0; | |
dbe3abbe | 359 | TString objname; |
75d81601 | 360 | while ((obj = next.Next())) { |
361 | objname = obj->GetName(); | |
dbe3abbe | 362 | TObjArray *toks = objname.Tokenize("_"); |
363 | if (toks->GetEntriesFast()){ | |
364 | TObjString *fpart = (TObjString *)(toks->UncheckedAt(0)); | |
75d81601 | 365 | if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++); |
259c3296 | 366 | } |
367 | } | |
75d81601 | 368 | fOutput->Add(fHistSECVTX); |
dbe3abbe | 369 | } |
809a4336 | 370 | } |
371 | ||
372 | //____________________________________________________________ | |
373 | void AliAnalysisTaskHFE::Exec(Option_t *){ | |
dbe3abbe | 374 | // |
375 | // Run the analysis | |
376 | // | |
377 | if(!fESD){ | |
378 | AliError("No ESD Event"); | |
379 | return; | |
380 | } | |
381 | if(!fMC){ | |
382 | AliError("No MC Event"); | |
383 | return; | |
384 | } | |
722347d8 | 385 | if(!fCuts){ |
386 | AliError("HFE cuts not available"); | |
387 | return; | |
388 | } | |
809a4336 | 389 | |
dbe3abbe | 390 | //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC); |
391 | Double_t container[6]; | |
78ea5ef4 | 392 | // container for the output THnSparse |
393 | Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B'] | |
809a4336 | 394 | |
dbe3abbe | 395 | // Loop over the Monte Carlo tracks to see whether we have overlooked any track |
396 | AliMCParticle *mctrack = 0x0; | |
397 | Int_t nElectrons = 0; | |
809a4336 | 398 | |
dbe3abbe | 399 | if (IsSecVtxOn()) { |
400 | fSecVtx->SetEvent(fESD); | |
401 | fSecVtx->SetStack(fMC->Stack()); | |
402 | } | |
259c3296 | 403 | |
dbe3abbe | 404 | // run mc QA ------------------------------------------------ |
75d81601 | 405 | if (IsQAOn(kMCqa)) { |
dbe3abbe | 406 | AliDebug(2, "Running MC QA"); |
259c3296 | 407 | |
dbe3abbe | 408 | fMCQA->SetStack(fMC->Stack()); |
409 | fMCQA->Init(); | |
259c3296 | 410 | |
259c3296 | 411 | Int_t nMCTracks = fMC->Stack()->GetNtrack(); |
412 | ||
259c3296 | 413 | // loop over all tracks for decayed electrons |
414 | for (Int_t igen = 0; igen < nMCTracks; igen++){ | |
722347d8 | 415 | TParticle* mcpart = fMC->Stack()->Particle(igen); |
416 | fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm); | |
417 | fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty); | |
418 | fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm); | |
419 | fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty); | |
420 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0); // no accept cut | |
421 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut | |
422 | if (TMath::Abs(mcpart->Eta()) < 0.9) { | |
423 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9 | |
424 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9 | |
425 | } | |
426 | if (TMath::Abs(GetRapidity(mcpart)) < 0.5) { | |
427 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5 | |
428 | fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5 | |
429 | } | |
259c3296 | 430 | } |
722347d8 | 431 | fMCQA->EndOfEventAna(AliHFEmcQA::kCharm); |
432 | fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty); | |
259c3296 | 433 | |
434 | } // end of MC QA loop | |
435 | // ----------------------------------------------------------------- | |
436 | ||
dbe3abbe | 437 | // |
438 | // Loop MC | |
439 | // | |
0792aa82 | 440 | fCuts->SetEventInfo(fMC); |
dbe3abbe | 441 | for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){ |
78ea5ef4 | 442 | mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(imc)); |
dbe3abbe | 443 | |
444 | container[0] = mctrack->Pt(); | |
445 | container[1] = mctrack->Eta(); | |
809a4336 | 446 | container[2] = mctrack->Phi(); |
447 | ||
dbe3abbe | 448 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue; |
449 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated); | |
722347d8 | 450 | if(IsSignalElectron(mctrack)) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal); |
dbe3abbe | 451 | (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi()); |
452 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue; | |
453 | // find the label in the vector | |
454 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance); | |
455 | nElectrons++; | |
456 | } | |
457 | (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons); | |
458 | ||
459 | // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD); | |
460 | ||
461 | ||
462 | // | |
463 | // Loop ESD | |
464 | // | |
75d81601 | 465 | Int_t nElectronCandidates = 0; |
dbe3abbe | 466 | AliESDtrack *track = 0x0, *htrack = 0x0; |
259c3296 | 467 | Int_t pid = 0; |
dbe3abbe | 468 | // For double counted tracks |
78ea5ef4 | 469 | LabelContainer cont(fESD->GetNumberOfTracks()); |
dbe3abbe | 470 | Bool_t alreadyseen = kFALSE; |
471 | ||
78ea5ef4 | 472 | Bool_t signal = kTRUE; |
473 | ||
0792aa82 | 474 | fCuts->SetEventInfo(fESD); |
dbe3abbe | 475 | for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){ |
476 | ||
477 | track = fESD->GetTrack(itrack); | |
478 | ||
479 | container[0] = track->Pt(); | |
480 | container[1] = track->Eta(); | |
809a4336 | 481 | container[2] = track->Phi(); |
dbe3abbe | 482 | |
78ea5ef4 | 483 | dataE[0] = track->Pt(); |
484 | dataE[1] = track->Eta(); | |
485 | dataE[2] = track->Phi(); | |
486 | dataE[3] = -1; | |
487 | dataE[4] = -1; | |
488 | ||
489 | signal = kTRUE; | |
dbe3abbe | 490 | |
722347d8 | 491 | // RecKine: ITSTPC cuts |
492 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) continue; | |
dbe3abbe | 493 | |
722347d8 | 494 | // Check if it is electrons near the vertex |
78ea5ef4 | 495 | if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue; |
0792aa82 | 496 | TParticle* mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel())); |
78ea5ef4 | 497 | |
dbe3abbe | 498 | container[3] = mctrack->Pt(); |
499 | container[4] = mctrack->Eta(); | |
500 | container[5] = mctrack->Phi(); | |
501 | ||
502 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; | |
503 | ||
504 | if(signal) { | |
78ea5ef4 | 505 | alreadyseen = cont.Find(TMath::Abs(track->GetLabel())); |
506 | cont.Append(TMath::Abs(track->GetLabel())); | |
dbe3abbe | 507 | |
722347d8 | 508 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecKineITSTPC); |
509 | fCFM->GetParticleContainer()->Fill(&container[0], (AliHFEcuts::kStepRecKineITSTPC + 2*(AliHFEcuts::kNcutESDSteps + 1))); | |
dbe3abbe | 510 | if(alreadyseen) { |
722347d8 | 511 | fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITSTPC + (AliHFEcuts::kNcutESDSteps + 1))); |
dbe3abbe | 512 | } |
513 | } | |
514 | ||
dbe3abbe | 515 | // Check TRD criterions (outside the correction framework) |
516 | if(track->GetTRDncls()){ | |
517 | (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls()); | |
518 | (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts | |
519 | (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality()); | |
520 | (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls()); | |
521 | } | |
522 | ||
523 | ||
524 | // RecPrim | |
525 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue; | |
dbe3abbe | 526 | if(signal) { |
722347d8 | 527 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1)); |
528 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim); | |
dbe3abbe | 529 | if(alreadyseen) { |
722347d8 | 530 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1); |
dbe3abbe | 531 | } |
532 | } | |
533 | ||
534 | ||
535 | // HFEcuts: ITS layers cuts | |
536 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue; | |
dbe3abbe | 537 | if(signal) { |
722347d8 | 538 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1)); |
539 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS); | |
dbe3abbe | 540 | if(alreadyseen) { |
722347d8 | 541 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1); |
dbe3abbe | 542 | } |
543 | } | |
544 | ||
78ea5ef4 | 545 | // HFEcuts: Nb of tracklets TRD0 |
dbe3abbe | 546 | if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue; |
dbe3abbe | 547 | if(signal) { |
722347d8 | 548 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)); |
549 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD); | |
dbe3abbe | 550 | if(alreadyseen) { |
722347d8 | 551 | fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + (AliHFEcuts::kNcutESDSteps + 1))); |
dbe3abbe | 552 | } |
78ea5ef4 | 553 | // dimensions 3&4&5 : pt,eta,phi (MC) |
722347d8 | 554 | ((THnSparseF *)fCorrelation->At(0))->Fill(container); |
dbe3abbe | 555 | } |
556 | ||
0792aa82 | 557 | if (IsQAOn(kMCqa)) { |
558 | // mc qa for after the reconstruction cuts | |
559 | AliDebug(2, "Running MC QA"); | |
560 | fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 3); // charm | |
561 | fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 3); // beauty | |
562 | } | |
dbe3abbe | 563 | |
809a4336 | 564 | // track accepted, do PID |
722347d8 | 565 | AliHFEpidObject hfetrack; |
566 | hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis; | |
567 | hfetrack.fRecTrack = track; | |
568 | if(HasMCData()) hfetrack.fMCtrack = mctrack; | |
569 | if(!fPID->IsSelected(&hfetrack)) continue; | |
75d81601 | 570 | nElectronCandidates++; |
dbe3abbe | 571 | |
0792aa82 | 572 | if (IsQAOn(kMCqa)) { |
573 | // mc qa for after the reconstruction and pid cuts | |
574 | AliDebug(2, "Running MC QA"); | |
575 | fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 4); // charm | |
576 | fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 4); // beauty | |
577 | } | |
dbe3abbe | 578 | |
579 | // Fill Containers | |
dbe3abbe | 580 | if(signal) { |
722347d8 | 581 | fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD +1 + 2*(AliHFEcuts::kNcutESDSteps + 1)); |
582 | fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 1); | |
dbe3abbe | 583 | if(alreadyseen) { |
722347d8 | 584 | fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + (AliHFEcuts::kNcutESDSteps + 1))); |
dbe3abbe | 585 | } |
259c3296 | 586 | // dimensions 3&4&5 : pt,eta,phi (MC) |
722347d8 | 587 | ((THnSparseF *)fCorrelation->At(1))->Fill(container); |
dbe3abbe | 588 | } |
589 | ||
dbe3abbe | 590 | // Track selected: distinguish between true and fake |
78ea5ef4 | 591 | AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode())); |
dbe3abbe | 592 | if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){ |
78ea5ef4 | 593 | Int_t type = IsSignalElectron(track); |
594 | AliDebug(1, Form("Type: %d\n", type)); | |
595 | if(type){ | |
596 | dataE[4] = type; // beauty[1] or charm[2] | |
597 | dataE[3] = 2; // signal electron | |
598 | } | |
599 | else{ | |
600 | dataE[3] = 1; // not a signal electron | |
601 | dataE[4] = 0; | |
602 | } | |
259c3296 | 603 | // pair analysis [mj] |
604 | if (IsSecVtxOn()) { | |
dbe3abbe | 605 | AliDebug(2, "Running Secondary Vertex Analysis"); |
606 | fSecVtx->InitAnaPair(); | |
259c3296 | 607 | for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){ |
608 | htrack = fESD->GetTrack(jtrack); | |
609 | if ( itrack == jtrack ) continue; | |
610 | //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue; | |
611 | if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue; | |
dbe3abbe | 612 | fSecVtx->AnaPair(track, htrack, jtrack); |
259c3296 | 613 | } |
dbe3abbe | 614 | // based on the partner of e info, you run secandary vertexing function |
615 | fSecVtx->RunSECVTX(track); | |
616 | } // end of pair analysis | |
78ea5ef4 | 617 | } |
618 | else { | |
259c3296 | 619 | // Fill THnSparse with the information for Fake Electrons |
78ea5ef4 | 620 | dataE[3] = 0; |
621 | dataE[4] = 0; | |
622 | } | |
623 | // fill the performance THnSparse, if the mc origin could be defined | |
624 | if(dataE[3] > -1){ | |
625 | AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4])); | |
626 | fPIDperformance->Fill(dataE); | |
259c3296 | 627 | } |
dbe3abbe | 628 | } |
78ea5ef4 | 629 | |
630 | ||
dbe3abbe | 631 | fNEvents->Fill(1); |
75d81601 | 632 | fNElectronTracksEvent->Fill(nElectronCandidates); |
dbe3abbe | 633 | |
634 | // Done!!! | |
635 | PostData(0, fNEvents); | |
636 | PostData(1, fOutput); | |
637 | PostData(2, fQA); | |
809a4336 | 638 | } |
639 | ||
640 | //____________________________________________________________ | |
641 | void AliAnalysisTaskHFE::Terminate(Option_t *){ | |
dbe3abbe | 642 | // |
643 | // Terminate not implemented at the moment | |
644 | // | |
6ad05e72 | 645 | if(IsRunningPostProcess()){ |
646 | fOutput = dynamic_cast<TList *>(GetOutputData(1)); | |
647 | if(!fOutput){ | |
648 | AliError("Results not available"); | |
649 | return; | |
650 | } | |
651 | PostProcess(); | |
652 | } | |
653 | } | |
654 | ||
655 | //____________________________________________________________ | |
656 | void AliAnalysisTaskHFE::Load(TString filename){ | |
657 | // | |
658 | // Load Results into the task | |
659 | // | |
660 | fQA = NULL; fOutput = NULL; fNEvents = NULL; | |
661 | TFile *input = TFile::Open(filename.Data()); | |
662 | if(!input || input->IsZombie()){ | |
663 | AliError("Cannot read file"); | |
664 | return; | |
665 | } | |
666 | TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents")); | |
667 | if(htmp) | |
668 | fNEvents = dynamic_cast<TH1I *>(htmp->Clone()); | |
669 | else | |
670 | AliError("Event Counter histogram not found"); | |
671 | TList *ltmp = dynamic_cast<TList *>(input->Get("Results")); | |
672 | if(ltmp) | |
673 | fOutput = dynamic_cast<TList *>(ltmp->Clone()); | |
674 | else | |
675 | AliError("Output Histograms not found"); | |
676 | ltmp = dynamic_cast<TList *>(input->Get("QA")); | |
677 | if(ltmp) | |
678 | fQA = dynamic_cast<TList *>(ltmp->Clone()); | |
679 | else | |
680 | AliError("QA histograms not found"); | |
681 | input->Close(); | |
682 | delete input; | |
683 | Int_t nObjects = 0; | |
684 | if(fNEvents) nObjects++; | |
685 | if(fOutput) nObjects++; | |
686 | if(fQA) nObjects++; | |
687 | AliInfo(Form("Loaded %d Objects into task", nObjects)); | |
688 | } | |
689 | ||
690 | //____________________________________________________________ | |
691 | void AliAnalysisTaskHFE::PostProcess(){ | |
692 | // | |
693 | // Plotting Ratio histograms | |
694 | // + All electrons / all candidates (Purity for Electrons) | |
695 | // + All signal electrons / all electrons (Purity for signals) | |
696 | // For this the following pt-histograms have to be projected from the THnSparse | |
697 | // + All Electron candidates | |
698 | // + All Real electrons | |
699 | // + All Signal Electrons | |
700 | // + All misidentified electrons | |
722347d8 | 701 | // Additionally we draw Efficiency histograms: |
702 | // + Reconstructible Electrons | |
703 | // + Signal Electrons | |
704 | // + Selected Electrons | |
705 | // + Selected Electrons in Acceptance | |
6ad05e72 | 706 | // |
722347d8 | 707 | |
6ad05e72 | 708 | fPIDperformance = dynamic_cast<THnSparseF *>(fOutput->FindObject("PIDperformance")); |
709 | if(!fPIDperformance){ | |
710 | AliError("PID Performance histogram not found in the output container"); | |
711 | return; | |
712 | } | |
6ad05e72 | 713 | // Make projection |
714 | // always project to pt dimension | |
715 | // get the histograms under our control | |
716 | TH1 *allCandidates = NULL, *allElectrons = NULL, *allSignals = NULL, *allFakes = NULL; | |
717 | allCandidates = fPIDperformance->Projection(0); | |
718 | allCandidates->SetName("hAllCandidates"); | |
719 | allCandidates->SetTitle("All Candidates"); | |
722347d8 | 720 | allCandidates->Sumw2(); |
6ad05e72 | 721 | Int_t firstDim3 = fPIDperformance->GetAxis(3)->GetFirst(), lastDim3 = fPIDperformance->GetAxis(3)->GetLast(); |
722 | fPIDperformance->GetAxis(3)->SetRange(firstDim3 + 1, lastDim3); | |
723 | allElectrons = fPIDperformance->Projection(0); | |
722347d8 | 724 | allElectrons->Sumw2(); |
6ad05e72 | 725 | allElectrons->SetName("hAllElectrons"); |
726 | allElectrons->SetTitle("All Electrons"); | |
727 | Int_t firstDim4 = fPIDperformance->GetAxis(4)->GetFirst(), lastDim4 = fPIDperformance->GetAxis(4)->GetLast(); | |
728 | fPIDperformance->GetAxis(4)->SetRange(firstDim4 + 1, lastDim4); | |
729 | allSignals = fPIDperformance->Projection(0); | |
722347d8 | 730 | allSignals->Sumw2(); |
6ad05e72 | 731 | allSignals->SetName("hAllSignals"); |
732 | allSignals->SetTitle("All Signal Electrons"); | |
733 | fPIDperformance->GetAxis(4)->SetRange(firstDim4, lastDim4); // Reset 4th axis | |
722347d8 | 734 | fPIDperformance->GetAxis(3)->SetRange(firstDim3, firstDim3); // Select fakes |
6ad05e72 | 735 | allFakes = fPIDperformance->Projection(0); |
722347d8 | 736 | allFakes->Sumw2(); |
6ad05e72 | 737 | allFakes->SetName("hAllFakes"); |
738 | allFakes->SetTitle("All Fakes"); | |
739 | fPIDperformance->GetAxis(3)->SetRange(firstDim3, lastDim3); // Reset also 3rd axis | |
740 | ||
741 | // Make Ratios | |
742 | TH1D *electronPurity = dynamic_cast<TH1D *>(allElectrons->Clone()); | |
743 | electronPurity->Divide(allCandidates); | |
744 | electronPurity->SetName("electronPurity"); | |
745 | electronPurity->SetTitle("Electron Purity"); | |
722347d8 | 746 | electronPurity->GetXaxis()->SetTitle("p_{T} / GeV/c"); |
6ad05e72 | 747 | electronPurity->GetYaxis()->SetTitle("Purity / %"); |
748 | TH1D *signalPurity = dynamic_cast<TH1D *>(allSignals->Clone()); | |
749 | signalPurity->Divide(allElectrons); | |
750 | signalPurity->SetName("signalPurity"); | |
751 | signalPurity->SetTitle("Purity of Electrons coming from Heavy flavours"); | |
722347d8 | 752 | signalPurity->GetXaxis()->SetTitle("p_{T} / GeV/c"); |
6ad05e72 | 753 | signalPurity->GetYaxis()->SetTitle("Purity / %"); |
754 | TH1D *fakeContamination = dynamic_cast<TH1D *>(allFakes->Clone()); | |
755 | fakeContamination->Divide(allCandidates); | |
756 | fakeContamination->SetName("fakeContamination"); | |
757 | fakeContamination->SetTitle("Contamination of misidentified hadrons"); | |
722347d8 | 758 | fakeContamination->GetXaxis()->SetTitle("p_{T} / GeV/c"); |
6ad05e72 | 759 | fakeContamination->GetYaxis()->SetTitle("Purity / %"); |
722347d8 | 760 | |
761 | // Graphics settings | |
762 | const Int_t nHistos = 7; | |
763 | TH1 *hptr[7] = {allCandidates, allElectrons, allSignals, allFakes, electronPurity, signalPurity, fakeContamination}, *histo; | |
764 | for(Int_t ihist = 0; ihist < nHistos; ihist++){ | |
765 | (histo = hptr[ihist])->SetStats(kFALSE); | |
766 | histo->SetLineColor(kBlack); | |
767 | histo->SetMarkerColor(kBlue); | |
768 | histo->SetMarkerStyle(22); | |
769 | } | |
770 | ||
771 | // Make percent output | |
772 | electronPurity->Scale(100); | |
773 | signalPurity->Scale(100); | |
774 | fakeContamination->Scale(100); | |
6ad05e72 | 775 | |
776 | // Draw output | |
777 | TCanvas *cSpectra = new TCanvas("cSpectra","pt Spectra", 800, 600); | |
778 | cSpectra->Divide(2,2); | |
779 | TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600); | |
780 | cRatios->Divide(2,2); | |
781 | cSpectra->cd(1); | |
722347d8 | 782 | gPad->SetLogy(); |
783 | allCandidates->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
784 | allCandidates->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c"); | |
785 | allCandidates->Draw("e"); | |
6ad05e72 | 786 | cSpectra->cd(2); |
722347d8 | 787 | gPad->SetLogy(); |
788 | allElectrons->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
789 | allElectrons->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c"); | |
790 | allElectrons->Draw("e"); | |
6ad05e72 | 791 | cSpectra->cd(3); |
722347d8 | 792 | gPad->SetLogy(); |
793 | allSignals->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
794 | allSignals->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c"); | |
795 | allSignals->Draw("e"); | |
6ad05e72 | 796 | cSpectra->cd(4); |
722347d8 | 797 | gPad->SetLogy(); |
798 | allFakes->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
799 | allFakes->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c"); | |
800 | allFakes->Draw("e"); | |
6ad05e72 | 801 | cRatios->cd(1); |
722347d8 | 802 | electronPurity->Draw("e"); |
6ad05e72 | 803 | cRatios->cd(2); |
722347d8 | 804 | signalPurity->Draw("e"); |
6ad05e72 | 805 | cRatios->cd(3); |
722347d8 | 806 | fakeContamination->Draw("e"); |
807 | ||
808 | // Now draw the Efficiency | |
809 | // We show: | |
810 | // + InAcceptance / Generated | |
811 | // + Signal / Generated | |
812 | // + Selected / Generated | |
813 | // + Selected / InAcceptance (Reconstructible) | |
814 | TCanvas *cEff = new TCanvas("cEff", "Efficiency", 800, 600); | |
815 | cEff->Divide(2,2); | |
816 | AliCFContainer *effc = dynamic_cast<AliCFContainer *>(fOutput->At(0)); | |
817 | if(!effc){ | |
818 | AliError("Efficiency Container not found in Output Container"); | |
819 | return; | |
820 | } | |
821 | AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *effc); | |
822 | effCalc->CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated); | |
823 | TH1 *effReconstructibleP = effCalc->Project(0); | |
824 | effReconstructibleP->SetName("effReconstructibleP"); | |
825 | effReconstructibleP->SetTitle("Efficiency of reconstructible tracks"); | |
826 | effReconstructibleP->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
827 | effReconstructibleP->GetYaxis()->SetTitle("Efficiency"); | |
828 | effReconstructibleP->GetYaxis()->SetRangeUser(0.,1.); | |
829 | effReconstructibleP->SetMarkerStyle(22); | |
830 | effReconstructibleP->SetMarkerColor(kBlue); | |
831 | effReconstructibleP->SetLineColor(kBlack); | |
832 | effReconstructibleP->SetStats(kFALSE); | |
833 | cEff->cd(1); | |
834 | effReconstructibleP->Draw("e"); | |
835 | effCalc->CalculateEfficiency(AliHFEcuts::kStepMCsignal, AliHFEcuts::kStepMCGenerated); | |
836 | TH1 *effSignal = effCalc->Project(0); | |
837 | effSignal->SetName("effSignal"); | |
838 | effSignal->SetTitle("Efficiency of Signal Electrons"); | |
839 | effSignal->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
840 | effSignal->GetYaxis()->SetTitle("Efficiency"); | |
841 | effSignal->GetYaxis()->SetRangeUser(0., 1.); | |
842 | effSignal->SetMarkerStyle(22); | |
843 | effSignal->SetMarkerColor(kBlue); | |
844 | effSignal->SetLineColor(kBlack); | |
845 | effSignal->SetStats(kFALSE); | |
846 | cEff->cd(2); | |
847 | effSignal->Draw("e"); | |
848 | effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCGenerated); | |
849 | TH1 *effPIDP = effCalc->Project(0); | |
850 | effPIDP->SetName("effPIDP"); | |
851 | effPIDP->SetTitle("Efficiency of selected tracks"); | |
852 | effPIDP->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
853 | effPIDP->GetYaxis()->SetTitle("Efficiency"); | |
854 | effPIDP->GetYaxis()->SetRangeUser(0.,1.); | |
855 | effPIDP->SetMarkerStyle(22); | |
856 | effPIDP->SetMarkerColor(kBlue); | |
857 | effPIDP->SetLineColor(kBlack); | |
858 | effPIDP->SetStats(kFALSE); | |
859 | cEff->cd(3); | |
860 | effPIDP->Draw("e"); | |
861 | effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCInAcceptance); | |
862 | TH1 *effPIDAcc = effCalc->Project(0); | |
863 | effPIDAcc->SetName("effPIDAcc"); | |
864 | effPIDAcc->SetTitle("Efficiency of selected tracks in acceptance"); | |
865 | effPIDAcc->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
866 | effPIDAcc->GetYaxis()->SetTitle("Efficiency"); | |
867 | effPIDAcc->GetYaxis()->SetRangeUser(0.,1.); | |
868 | effPIDAcc->SetMarkerStyle(22); | |
869 | effPIDAcc->SetMarkerColor(kBlue); | |
870 | effPIDAcc->SetLineColor(kBlack); | |
871 | effPIDAcc->SetStats(kFALSE); | |
872 | cEff->cd(4); | |
873 | effPIDAcc->Draw("e"); | |
874 | delete effCalc; | |
6ad05e72 | 875 | |
876 | // cleanup | |
877 | //delete allCandidates; delete allElectrons; delete allSignals; delete allFakes; | |
878 | //delete electronPurity; delete signalPurity; delete fakeContamination; | |
809a4336 | 879 | } |
880 | ||
881 | //____________________________________________________________ | |
882 | void AliAnalysisTaskHFE::MakeParticleContainer(){ | |
883 | // | |
884 | // Create the particle container for the correction framework manager and | |
885 | // link it | |
886 | // | |
259c3296 | 887 | const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi |
dbe3abbe | 888 | const Double_t kPtmin = 0.1, kPtmax = 10.; |
259c3296 | 889 | const Double_t kEtamin = -0.9, kEtamax = 0.9; |
890 | const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi(); | |
809a4336 | 891 | |
892 | //arrays for the number of bins in each dimension | |
259c3296 | 893 | Int_t iBin[kNvar]; |
dbe3abbe | 894 | iBin[0] = 40; //bins in pt |
809a4336 | 895 | iBin[1] = 8; //bins in eta |
896 | iBin[2] = 18; // bins in phi | |
897 | ||
898 | //arrays for lower bounds : | |
259c3296 | 899 | Double_t* binEdges[kNvar]; |
900 | for(Int_t ivar = 0; ivar < kNvar; ivar++) | |
901 | binEdges[ivar] = new Double_t[iBin[ivar] + 1]; | |
809a4336 | 902 | |
903 | //values for bin lower bounds | |
dbe3abbe | 904 | 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 | 905 | for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i; |
906 | for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i; | |
809a4336 | 907 | |
908 | //one "container" for MC | |
722347d8 | 909 | AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 2*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin); |
dbe3abbe | 910 | |
809a4336 | 911 | //setting the bin limits |
dbe3abbe | 912 | for(Int_t ivar = 0; ivar < kNvar; ivar++) |
913 | container -> SetBinLimits(ivar, binEdges[ivar]); | |
809a4336 | 914 | fCFM->SetParticleContainer(container); |
259c3296 | 915 | |
916 | //create correlation matrix for unfolding | |
917 | Int_t thnDim[2*kNvar]; | |
918 | for (int k=0; k<kNvar; k++) { | |
919 | //first half : reconstructed | |
920 | //second half : MC | |
921 | thnDim[k] = iBin[k]; | |
922 | thnDim[k+kNvar] = iBin[k]; | |
923 | } | |
924 | ||
78ea5ef4 | 925 | if(!fCorrelation) fCorrelation = new TList(); |
926 | fCorrelation->SetName("correlation"); | |
927 | ||
722347d8 | 928 | THnSparseF *correlation0 = new THnSparseF("correlationstepbeforePID","THnSparse with correlations",2*kNvar,thnDim); |
929 | THnSparseF *correlation1 = new THnSparseF("correlationstepafterPID","THnSparse with correlations",2*kNvar,thnDim); | |
259c3296 | 930 | for (int k=0; k<kNvar; k++) { |
78ea5ef4 | 931 | correlation0->SetBinEdges(k,binEdges[k]); |
932 | correlation0->SetBinEdges(k+kNvar,binEdges[k]); | |
933 | correlation1->SetBinEdges(k,binEdges[k]); | |
934 | correlation1->SetBinEdges(k+kNvar,binEdges[k]); | |
259c3296 | 935 | } |
78ea5ef4 | 936 | correlation0->Sumw2(); |
937 | correlation1->Sumw2(); | |
722347d8 | 938 | |
78ea5ef4 | 939 | fCorrelation->AddAt(correlation0,0); |
940 | fCorrelation->AddAt(correlation1,1); | |
78ea5ef4 | 941 | |
259c3296 | 942 | // Add a histogram for Fake electrons |
78ea5ef4 | 943 | const Int_t nDim=5; |
944 | Int_t nBin[nDim] = {40, 8, 18, 3, 3}; | |
945 | Double_t* binEdges2[nDim]; | |
946 | for(Int_t ivar = 0; ivar < nDim; ivar++) | |
947 | binEdges2[ivar] = new Double_t[nBin[ivar] + 1]; | |
948 | ||
949 | //values for bin lower bounds | |
950 | for(Int_t i=0; i<=nBin[0]; i++) binEdges2[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/nBin[0]*(Double_t)i); | |
951 | for(Int_t i=0; i<=nBin[1]; i++) binEdges2[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/nBin[1]*(Double_t)i; | |
952 | for(Int_t i=0; i<=nBin[2]; i++) binEdges2[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/nBin[2]*(Double_t)i; | |
953 | for(Int_t i=0; i<=nBin[3]; i++) binEdges2[3][i] = i; | |
954 | for(Int_t i=0; i<=nBin[4]; i++) binEdges2[4][i] = i; | |
955 | ||
956 | fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad] type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin); | |
957 | for(Int_t idim = 0; idim < nDim; idim++) | |
958 | fPIDperformance->SetBinEdges(idim, binEdges2[idim]); | |
809a4336 | 959 | } |
dbe3abbe | 960 | |
50685501 | 961 | //____________________________________________________________ |
75d81601 | 962 | void AliAnalysisTaskHFE::AddPIDdetector(Char_t *detector){ |
50685501 | 963 | // |
964 | // Adding PID detector to the task | |
965 | // | |
75d81601 | 966 | if(!fPIDdetectors.Length()) |
967 | fPIDdetectors = detector; | |
968 | else | |
969 | fPIDdetectors += Form(":%s", detector); | |
970 | } | |
971 | ||
972 | //____________________________________________________________ | |
50685501 | 973 | void AliAnalysisTaskHFE::PrintStatus() const { |
78ea5ef4 | 974 | // |
975 | // Print Analysis status | |
976 | // | |
75d81601 | 977 | printf("\n\tAnalysis Settings\n\t========================================\n\n"); |
978 | printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO"); | |
979 | printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO"); | |
980 | printf("\n"); | |
981 | printf("\tParticle Identification Detectors:\n"); | |
982 | TObjArray *detectors = fPIDdetectors.Tokenize(":"); | |
983 | for(Int_t idet = 0; idet < detectors->GetEntries(); idet++) | |
984 | printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data()); | |
985 | printf("\n"); | |
986 | printf("\tQA: \n"); | |
987 | printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" : "NO"); | |
0792aa82 | 988 | printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsInDebugMode()) ? "YES" : "NO"); |
75d81601 | 989 | printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO"); |
990 | printf("\n"); | |
991 | } | |
78ea5ef4 | 992 | |
993 | //____________________________________________________________ | |
994 | AliAnalysisTaskHFE::LabelContainer::LabelContainer(Int_t capacity): | |
995 | fContainer(NULL), | |
996 | fBegin(NULL), | |
997 | fEnd(NULL), | |
998 | fLast(NULL), | |
999 | fCurrent(NULL) | |
1000 | { | |
1001 | // | |
1002 | // Default constructor | |
1003 | // | |
1004 | fContainer = new Int_t[capacity]; | |
1005 | fBegin = &fContainer[0]; | |
1006 | fEnd = &fContainer[capacity - 1]; | |
1007 | fLast = fCurrent = fBegin; | |
1008 | } | |
1009 | ||
1010 | //____________________________________________________________ | |
1011 | Bool_t AliAnalysisTaskHFE::LabelContainer::Append(Int_t label){ | |
1012 | // | |
1013 | // Add Label to the container | |
1014 | // | |
1015 | if(fLast > fEnd) return kFALSE; | |
1016 | *fLast++ = label; | |
1017 | return kTRUE; | |
1018 | } | |
1019 | ||
1020 | //____________________________________________________________ | |
50685501 | 1021 | Bool_t AliAnalysisTaskHFE::LabelContainer::Find(Int_t label) const { |
78ea5ef4 | 1022 | // |
1023 | // Find track in the list of labels | |
1024 | // | |
1025 | for(Int_t *entry = fBegin; entry <= fLast; entry++) | |
1026 | if(*entry == label) return kTRUE; | |
1027 | return kFALSE; | |
1028 | } | |
1029 | ||
1030 | //____________________________________________________________ | |
1031 | Int_t AliAnalysisTaskHFE::LabelContainer::Next(){ | |
1032 | // | |
1033 | // Mimic iterator | |
1034 | // | |
1035 | if(fCurrent > fLast) return -1; | |
1036 | return *fCurrent++; | |
1037 | } | |
1038 | ||
1039 | //____________________________________________________________ | |
722347d8 | 1040 | Int_t AliAnalysisTaskHFE::IsSignalElectron(AliVParticle *fTrack) const{ |
78ea5ef4 | 1041 | // |
1042 | // Checks whether the identified electron track is coming from heavy flavour | |
1043 | // returns 0 in case of no signal, 1 in case of charm and 2 in case of Bottom | |
1044 | // | |
1045 | enum{ | |
1046 | kNoSignal = 0, | |
722347d8 | 1047 | kCharm = 1, |
1048 | kBeauty = 2 | |
78ea5ef4 | 1049 | }; |
722347d8 | 1050 | AliMCParticle *mctrack = NULL; |
1051 | TString objname = fTrack->IsA()->GetName(); | |
1052 | if(!objname.CompareTo("AliESDtrack")){ | |
1053 | AliDebug(2, "Checking signal for ESD track"); | |
1054 | AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack); | |
1055 | mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel()))); | |
1056 | } | |
1057 | else if(!objname.CompareTo("AliAODtrack")){ | |
1058 | AliError("AOD Analysis not implemented yet"); | |
1059 | return kNoSignal; | |
1060 | } | |
1061 | else if(!objname.CompareTo("AliMCParticle")){ | |
1062 | AliDebug(2, "Checking signal for MC track"); | |
1063 | mctrack = dynamic_cast<AliMCParticle *>(fTrack); | |
1064 | } | |
1065 | else{ | |
1066 | AliError("Input object not supported"); | |
1067 | return kNoSignal; | |
1068 | } | |
78ea5ef4 | 1069 | if(!mctrack) return kNoSignal; |
1070 | TParticle *ecand = mctrack->Particle(); | |
1071 | if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron | |
1072 | Int_t motherLabel = TMath::Abs(ecand->GetFirstMother()); | |
1073 | AliDebug(3, Form("mother label: %d\n", motherLabel)); | |
1074 | if(!motherLabel) return kNoSignal; // mother track unknown | |
1075 | AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel)); | |
1076 | if(!motherTrack) return kNoSignal; | |
1077 | TParticle *mparticle = motherTrack->Particle(); | |
1078 | Int_t pid = TMath::Abs(mparticle->GetPdgCode()); | |
1079 | AliDebug(3, Form("PDG code: %d\n", pid)); | |
1080 | ||
1081 | // identify signal according to Pdg Code | |
1082 | if((pid % 1000) / 100 == 4) return kCharm; // charmed meson, 3rd position in pdg code == 4 | |
1083 | if(pid / 1000 == 4) return kCharm; // charmed baryon, 4th position in pdg code == 4 | |
1084 | if((pid % 1000) / 100 == 5) return kBeauty; // beauty meson, 3rd position in pdg code == 5 | |
1085 | if(pid / 1000 == 5) return kBeauty; // beauty baryon, 4th position in pdg code == 5 | |
1086 | return kNoSignal; | |
1087 | } | |
722347d8 | 1088 | |
1089 | //__________________________________________ | |
50685501 | 1090 | Float_t AliAnalysisTaskHFE::GetRapidity(TParticle *part) const |
722347d8 | 1091 | { |
50685501 | 1092 | // |
1093 | // return rapidity | |
1094 | // | |
1095 | Float_t rapidity; | |
1096 | if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999; | |
1097 | else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); | |
1098 | return rapidity; | |
722347d8 | 1099 | } |
1100 |