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