]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ANALYSIS/AliAnalysisTaskPIDqa.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Coverity fixes
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskPIDqa.cxx
... / ...
CommitLineData
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/* $Id: AliAnalysisTaskPIDqa.cxx 43811 2010-09-23 14:13:31Z wiechula $ */
17#include <TList.h>
18#include <TVectorD.h>
19#include <TObjArray.h>
20#include <TH2.h>
21#include <TFile.h>
22#include <TPRegexp.h>
23#include <TChain.h>
24#include <TF1.h>
25#include <TSpline.h>
26
27#include <AliAnalysisManager.h>
28#include <AliInputEventHandler.h>
29#include <AliVEventHandler.h>
30#include <AliVEvent.h>
31#include <AliVParticle.h>
32#include <AliVTrack.h>
33#include <AliLog.h>
34#include <AliPID.h>
35#include <AliPIDResponse.h>
36#include <AliITSPIDResponse.h>
37#include <AliTPCPIDResponse.h>
38#include <AliTRDPIDResponse.h>
39#include <AliTOFPIDResponse.h>
40
41#include <AliESDEvent.h>
42
43#include "AliAnalysisTaskPIDqa.h"
44
45
46ClassImp(AliAnalysisTaskPIDqa)
47
48//______________________________________________________________________________
49AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
50AliAnalysisTaskSE(),
51fPIDResponse(0x0),
52fListQA(0x0),
53fListQAits(0x0),
54fListQAitsSA(0x0),
55fListQAitsPureSA(0x0),
56fListQAtpc(0x0),
57fListQAtrd(0x0),
58fListQAtof(0x0),
59fListQAemcal(0x0),
60fListQAhmpid(0x0),
61fListQAtofhmpid(0x0),
62fListQAtpctof(0x0)
63{
64 //
65 // Dummy constructor
66 //
67}
68
69//______________________________________________________________________________
70AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
71AliAnalysisTaskSE(name),
72fPIDResponse(0x0),
73fListQA(0x0),
74fListQAits(0x0),
75fListQAitsSA(0x0),
76fListQAitsPureSA(0x0),
77fListQAtpc(0x0),
78fListQAtrd(0x0),
79fListQAtof(0x0),
80fListQAemcal(0x0),
81fListQAhmpid(0x0),
82fListQAtofhmpid(0x0),
83fListQAtpctof(0x0)
84{
85 //
86 // Default constructor
87 //
88 DefineInput(0,TChain::Class());
89 DefineOutput(1,TList::Class());
90}
91
92//______________________________________________________________________________
93AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
94{
95 //
96 // Destructor
97 //
98 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA;
99}
100
101//______________________________________________________________________________
102void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
103{
104 //
105 // Create the output QA objects
106 //
107
108 AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
109
110 //input hander
111 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
112 AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
113 if (!inputHandler) AliFatal("Input handler needed");
114
115 //pid response object
116 fPIDResponse=inputHandler->GetPIDResponse();
117 if (!fPIDResponse) AliError("PIDResponse object was not created");
118
119 //
120 fListQA=new TList;
121 fListQA->SetOwner();
122
123 fListQAits=new TList;
124 fListQAits->SetOwner();
125 fListQAits->SetName("ITS");
126
127 fListQAitsSA=new TList;
128 fListQAitsSA->SetOwner();
129 fListQAitsSA->SetName("ITS_SA");
130
131 fListQAitsPureSA=new TList;
132 fListQAitsPureSA->SetOwner();
133 fListQAitsPureSA->SetName("ITS_PureSA");
134
135 fListQAtpc=new TList;
136 fListQAtpc->SetOwner();
137 fListQAtpc->SetName("TPC");
138
139 fListQAtrd=new TList;
140 fListQAtrd->SetOwner();
141 fListQAtrd->SetName("TRD");
142
143 fListQAtof=new TList;
144 fListQAtof->SetOwner();
145 fListQAtof->SetName("TOF");
146
147 fListQAemcal=new TList;
148 fListQAemcal->SetOwner();
149 fListQAemcal->SetName("EMCAL");
150
151 fListQAhmpid=new TList;
152 fListQAhmpid->SetOwner();
153 fListQAhmpid->SetName("HMPID");
154
155 fListQAtpctof=new TList;
156 fListQAtpctof->SetOwner();
157 fListQAtpctof->SetName("TPC_TOF");
158
159 fListQAtofhmpid=new TList;
160 fListQAtofhmpid->SetOwner();
161 fListQAtofhmpid->SetName("TOF_HMPID");
162
163 fListQA->Add(fListQAits);
164 fListQA->Add(fListQAitsSA);
165 fListQA->Add(fListQAitsPureSA);
166 fListQA->Add(fListQAtpc);
167 fListQA->Add(fListQAtrd);
168 fListQA->Add(fListQAtof);
169 fListQA->Add(fListQAemcal);
170 fListQA->Add(fListQAhmpid);
171 fListQA->Add(fListQAtpctof);
172 fListQA->Add(fListQAtofhmpid);
173
174 SetupITSqa();
175 SetupTPCqa();
176 SetupTRDqa();
177 SetupTOFqa();
178 SetupEMCALqa();
179 SetupHMPIDqa();
180 SetupTPCTOFqa();
181 SetupTOFHMPIDqa();
182
183 PostData(1,fListQA);
184}
185
186
187//______________________________________________________________________________
188void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
189{
190 //
191 // Setup the PID response functions and fill the QA histograms
192 //
193
194 AliVEvent *event=InputEvent();
195 if (!event||!fPIDResponse) return;
196
197
198 FillITSqa();
199 FillTPCqa();
200 FillTRDqa();
201 FillTOFqa();
202 FillEMCALqa();
203 FillHMPIDqa();
204
205 //combined detector QA
206 FillTPCTOFqa();
207 FillTOFHMPIDqa();
208
209 PostData(1,fListQA);
210}
211
212//______________________________________________________________________________
213void AliAnalysisTaskPIDqa::FillITSqa()
214{
215 //
216 // Fill PID qa histograms for the ITS
217 //
218
219 AliVEvent *event=InputEvent();
220
221 Int_t ntracks=event->GetNumberOfTracks();
222 for(Int_t itrack = 0; itrack < ntracks; itrack++){
223 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
224 ULong_t status=track->GetStatus();
225 // not that nice. status bits not in virtual interface
226 // ITS refit + ITS pid selection
227 if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
228 ! ( (status & AliVTrack::kITSpid )==AliVTrack::kITSpid ) )) continue;
229 Double_t mom=track->P();
230
231 TList *theList = 0x0;
232 if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
233 //ITS+TPC tracks
234 theList=fListQAits;
235 }else{
236 if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){
237 //ITS Standalone tracks
238 theList=fListQAitsSA;
239 }else{
240 //ITS Pure Standalone tracks
241 theList=fListQAitsPureSA;
242 }
243 }
244
245
246 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
247 TH2 *h=(TH2*)theList->At(ispecie);
248 if (!h) continue;
249 Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
250 h->Fill(mom,nSigma);
251 }
252 TH2 *h=(TH2*)theList->At(AliPID::kSPECIES);
253 if (h) {
254 Double_t sig=track->GetITSsignal();
255 h->Fill(mom,sig);
256 }
257 }
258}
259
260//______________________________________________________________________________
261void AliAnalysisTaskPIDqa::FillTPCqa()
262{
263 //
264 // Fill PID qa histograms for the TPC
265 //
266
267 AliVEvent *event=InputEvent();
268
269 Int_t ntracks=event->GetNumberOfTracks();
270 for(Int_t itrack = 0; itrack < ntracks; itrack++){
271 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
272
273 //
274 //basic track cuts
275 //
276 ULong_t status=track->GetStatus();
277 // not that nice. status bits not in virtual interface
278 // TPC refit + ITS refit + TPC pid
279 if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
280 !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
281
282 // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
283 //|| !( (status & AliVTrack::kTPCpid ) == AliVTrack::kTPCpid )
284 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
285 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
286 if (track->GetTPCNclsF()>0) {
287 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
288 }
289
290 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
291
292 Double_t mom=track->GetTPCmomentum();
293
294 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
295 TH2 *h=(TH2*)fListQAtpc->At(ispecie);
296 if (!h) continue;
297 Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
298 h->Fill(mom,nSigma);
299 }
300
301 TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
302 if (h) {
303 Double_t sig=track->GetTPCsignal();
304 h->Fill(mom,sig);
305 }
306 }
307}
308
309//______________________________________________________________________________
310void AliAnalysisTaskPIDqa::FillTRDqa()
311{
312 //
313 // Fill PID qa histograms for the TRD
314 //
315 AliVEvent *event=InputEvent();
316 Int_t ntracks = event->GetNumberOfTracks();
317 for(Int_t itrack = 0; itrack < ntracks; itrack++){
318 AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
319
320 //
321 //basic track cuts
322 //
323 ULong_t status=track->GetStatus();
324 // not that nice. status bits not in virtual interface
325 // TPC refit + ITS refit + TPC pid + TRD out
326 if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
327 !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
328// !( (status & AliVTrack::kTPCpid ) == AliVTrack::kTPCpid ) || //removes light nuclei. So it is out for the moment
329 !( (status & AliVTrack::kTRDout ) == AliVTrack::kTRDout )) continue;
330
331 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
332 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
333 if (track->GetTPCNclsF()>0) {
334 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
335 }
336
337 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
338
339 Double_t likelihoods[AliPID::kSPECIES];
340 if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
341 Int_t ntracklets = 0;
342 Double_t momentum = -1.;
343 for(Int_t itl = 0; itl < 6; itl++)
344 if(track->GetTRDmomentum(itl) > 0.){
345 ntracklets++;
346 if(momentum < 0) momentum = track->GetTRDmomentum(itl);
347 }
348 for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
349 TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
350 if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
351 }
352 }
353}
354
355//______________________________________________________________________________
356void AliAnalysisTaskPIDqa::FillTOFqa()
357{
358 //
359 // Fill TOF information
360 //
361 AliVEvent *event=InputEvent();
362
363 Int_t ntracks=event->GetNumberOfTracks();
364 Int_t tracksAtTof = 0;
365 for(Int_t itrack = 0; itrack < ntracks; itrack++){
366 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
367
368 //
369 //basic track cuts
370 //
371 ULong_t status=track->GetStatus();
372 // TPC refit + ITS refit +
373 // TOF out + TOFpid +
374 // kTIME
375 // (we don't use kTOFmismatch because it depends on TPC....)
376 if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
377 !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
378 !((status & AliVTrack::kTOFout ) == AliVTrack::kTOFout ) ||
379 !((status & AliVTrack::kTOFpid ) == AliVTrack::kTOFpid ) ||
380 !((status & AliVTrack::kTIME ) == AliVTrack::kTIME ) ) continue;
381
382 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
383 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
384 if (track->GetTPCNclsF()>0) {
385 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
386 }
387
388 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
389
390 tracksAtTof++;
391
392 Double_t mom=track->P();
393
394 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
395 TH2 *h=(TH2*)fListQAtof->At(ispecie);
396 if (!h) continue;
397 Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
398 h->Fill(mom,nSigma);
399 }
400
401 TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
402 if (h) {
403 Double_t sig=track->GetTOFsignal()/1000.;
404 h->Fill(mom,sig);
405 }
406
407 Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
408 ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
409
410 if (mom >= 0.75 && mom <= 1.25 ) {
411 Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
412 if (mask == 0) {
413 ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
414 } else if (mask == 1) {
415 ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
416 } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
417 ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
418 } else {
419 ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
420 }
421 }
422
423 Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
424 ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
425
426 AliESDEvent *esd = dynamic_cast<AliESDEvent *>(event);
427 if (esd) {
428 Double_t startTime = esd->GetT0TOF(0);
429 if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTime);
430 else {
431 startTime = esd->GetT0TOF(1);
432 if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTime);
433 startTime = esd->GetT0TOF(2);
434 if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTime);
435 }
436 }
437 }
438 if (tracksAtTof > 0) {
439 ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
440 Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
441 if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
442 }
443}
444
445
446//______________________________________________________________________________
447void AliAnalysisTaskPIDqa::FillEMCALqa()
448{
449 //
450 // Fill PID qa histograms for the EMCAL
451 //
452
453 AliVEvent *event=InputEvent();
454
455 Int_t ntracks=event->GetNumberOfTracks();
456 for(Int_t itrack = 0; itrack < ntracks; itrack++){
457 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
458
459 //
460 //basic track cuts
461 //
462 ULong_t status=track->GetStatus();
463 // not that nice. status bits not in virtual interface
464 // TPC refit + ITS refit +
465 // TOF out + TOFpid +
466 // kTIME
467 if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
468
469 Double_t pt=track->Pt();
470
471 //EMCAL nSigma (only for electrons at the moment)
472 TH2 *h=(TH2*)fListQAemcal->At(0);
473 if (!h) continue;
474 Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
475 h->Fill(pt,nSigma);
476
477 //EMCAL signal (E/p vs. pT)
478 h=(TH2*)fListQAemcal->At(1);
479 if (h) {
480
481 Int_t nMatchClus = track->GetEMCALcluster();
482 Double_t mom = track->P();
483 Double_t eop = -1.;
484
485 if(nMatchClus > -1){
486
487 AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
488
489 if(matchedClus){
490
491 // matched cluster is EMCAL
492 if(matchedClus->IsEMCAL()){
493
494 Double_t fClsE = matchedClus->E();
495 eop = fClsE/mom;
496
497 h->Fill(pt,eop);
498
499 }
500 }
501 }
502 }
503 }
504}
505
506
507//______________________________________________________________________________
508void AliAnalysisTaskPIDqa::FillHMPIDqa()
509{
510 //
511 // Fill PID qa histograms for the HMPID
512 //
513
514 AliVEvent *event=InputEvent();
515
516 Int_t ntracks=event->GetNumberOfTracks();
517 for(Int_t itrack = 0; itrack < ntracks; itrack++){
518 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
519
520 //
521 //basic track cuts
522 //
523 ULong_t status=track->GetStatus();
524 // not that nice. status bits not in virtual interface
525 // TPC refit + ITS refit +
526 // TOF out + TOFpid +
527 // kTIME
528 if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
529 !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
530
531 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
532 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
533 if (track->GetTPCNclsF()>0) {
534 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
535 }
536
537 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
538
539 Double_t mom = track->P();
540 Double_t ckovAngle = track->GetHMPIDsignal();
541
542 TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(0);;
543
544 hThetavsMom->Fill(mom,ckovAngle);
545
546 }
547}
548//______________________________________________________________________________
549void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
550{
551 //
552 // Fill PID qa histograms for the HMPID
553 //
554
555 AliVEvent *event=InputEvent();
556
557 Int_t ntracks=event->GetNumberOfTracks();
558 for(Int_t itrack = 0; itrack < ntracks; itrack++){
559 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
560
561 //
562 //basic track cuts
563 //
564 ULong_t status=track->GetStatus();
565 // not that nice. status bits not in virtual interface
566 // TPC refit + ITS refit +
567 // TOF out + TOFpid +
568 // kTIME
569 if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
570 !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
571 !((status & AliVTrack::kTOFout ) == AliVTrack::kTOFout ) ||
572 !((status & AliVTrack::kTOFpid ) == AliVTrack::kTOFpid ) ||
573 !((status & AliVTrack::kTIME ) == AliVTrack::kTIME ) ) continue;
574
575 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
576 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
577 if (track->GetTPCNclsF()>0) {
578 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
579 }
580
581 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
582
583 Double_t mom = track->P();
584 Double_t ckovAngle = track->GetHMPIDsignal();
585
586 Double_t nSigmaTOF[3];
587 TH1F *h[3];
588
589 for (Int_t ispecie=2; ispecie<AliPID::kSPECIES; ++ispecie){
590 //TOF nSigma
591 nSigmaTOF[ispecie]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
592 h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
593
594 if(TMath::Abs(nSigmaTOF[0])<2) h[0]->Fill(mom,ckovAngle);
595
596 if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3) h[1]->Fill(mom,ckovAngle);
597
598 if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3) h[2]->Fill(mom,ckovAngle);
599
600 }
601
602}
603
604//______________________________________________________________________________
605void AliAnalysisTaskPIDqa::FillTPCTOFqa()
606{
607 //
608 // Fill PID qa histograms for the TOF
609 // Here also the TPC histograms after TOF selection are filled
610 //
611
612 AliVEvent *event=InputEvent();
613
614 Int_t ntracks=event->GetNumberOfTracks();
615 for(Int_t itrack = 0; itrack < ntracks; itrack++){
616 AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
617
618 //
619 //basic track cuts
620 //
621 ULong_t status=track->GetStatus();
622 // not that nice. status bits not in virtual interface
623 // TPC refit + ITS refit +
624 // TOF out + TOFpid +
625 // kTIME
626 if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
627 !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
628// !( (status & AliVTrack::kTPCpid ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
629 !((status & AliVTrack::kTOFout ) == AliVTrack::kTOFout ) ||
630 !((status & AliVTrack::kTOFpid ) == AliVTrack::kTOFpid ) ||
631 !((status & AliVTrack::kTIME ) == AliVTrack::kTIME ) ) continue;
632
633 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
634 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
635 if (track->GetTPCNclsF()>0) {
636 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
637 }
638
639 if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
640
641
642 Double_t mom=track->P();
643 Double_t momTPC=track->GetTPCmomentum();
644
645 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
646 //TOF nSigma
647 Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
648 Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
649
650 //TPC after TOF cut
651 TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
652 if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
653
654 //TOF after TPC cut
655 h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);
656 if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
657
658 //EMCAL after TOF and TPC cut
659 h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);
660 if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
661
662 Int_t nMatchClus = track->GetEMCALcluster();
663 Double_t pt = track->Pt();
664 Double_t eop = -1.;
665
666 if(nMatchClus > -1){
667
668 AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
669
670 if(matchedClus){
671
672 // matched cluster is EMCAL
673 if(matchedClus->IsEMCAL()){
674
675 Double_t fClsE = matchedClus->E();
676 eop = fClsE/mom;
677
678 h->Fill(pt,eop);
679
680
681 }
682 }
683 }
684 }
685 }
686 }
687}
688
689//______________________________________________________________________________
690void AliAnalysisTaskPIDqa::SetupITSqa()
691{
692 //
693 // Create the ITS qa objects
694 //
695
696 TVectorD *vX=MakeLogBinning(200,.1,30);
697
698 //ITS+TPC tracks
699 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
700 TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
701 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
702 vX->GetNrows()-1,vX->GetMatrixArray(),
703 200,-10,10);
704 fListQAits->Add(hNsigmaP);
705 }
706 TH2F *hSig = new TH2F("hSigP_ITS",
707 "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
708 vX->GetNrows()-1,vX->GetMatrixArray(),
709 300,0,300);
710 fListQAits->Add(hSig);
711
712 //ITS Standalone tracks
713 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
714 TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
715 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
716 vX->GetNrows()-1,vX->GetMatrixArray(),
717 200,-10,10);
718 fListQAitsSA->Add(hNsigmaPSA);
719 }
720 TH2F *hSigSA = new TH2F("hSigP_ITSSA",
721 "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
722 vX->GetNrows()-1,vX->GetMatrixArray(),
723 300,0,300);
724 fListQAitsSA->Add(hSigSA);
725
726 //ITS Pure Standalone tracks
727 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
728 TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
729 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
730 vX->GetNrows()-1,vX->GetMatrixArray(),
731 200,-10,10);
732 fListQAitsPureSA->Add(hNsigmaPPureSA);
733 }
734 TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
735 "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
736 vX->GetNrows()-1,vX->GetMatrixArray(),
737 300,0,300);
738 fListQAitsPureSA->Add(hSigPureSA);
739
740 delete vX;
741}
742
743//______________________________________________________________________________
744void AliAnalysisTaskPIDqa::SetupTPCqa()
745{
746 //
747 // Create the TPC qa objects
748 //
749
750 TVectorD *vX=MakeLogBinning(200,.1,30);
751
752 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
753 TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
754 Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
755 vX->GetNrows()-1,vX->GetMatrixArray(),
756 200,-10,10);
757 fListQAtpc->Add(hNsigmaP);
758 }
759
760
761 TH2F *hSig = new TH2F("hSigP_TPC",
762 "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
763 vX->GetNrows()-1,vX->GetMatrixArray(),
764 300,0,300);
765 fListQAtpc->Add(hSig);
766
767 delete vX;
768}
769
770//______________________________________________________________________________
771void AliAnalysisTaskPIDqa::SetupTRDqa()
772{
773 //
774 // Create the TRD qa objects
775 //
776 TVectorD *vX=MakeLogBinning(200,.1,30);
777 for(Int_t itl = 0; itl < 6; ++itl){
778 for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
779 TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
780 Form("TRD Likelihood to be %s %s for tracks having %d %s; p (GeV/c); TRD %s Likelihood", ispecie == 0 ? "an" : "a", AliPID::ParticleName(ispecie), itl+1, itl == 0 ? "tracklet" : "tracklets", AliPID::ParticleName(ispecie)),
781 vX->GetNrows()-1, vX->GetMatrixArray(),
782 100, 0., 1.);
783 fListQAtrd->Add(hLikeP);
784 }
785 }
786 delete vX;
787}
788
789//______________________________________________________________________________
790void AliAnalysisTaskPIDqa::SetupTOFqa()
791{
792 //
793 // Create the TOF qa objects
794 //
795
796 TVectorD *vX=MakeLogBinning(200,.1,30);
797
798 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
799 TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
800 Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
801 vX->GetNrows()-1,vX->GetMatrixArray(),
802 200,-10,10);
803 fListQAtof->Add(hNsigmaP);
804 }
805
806 // for Kaons PID we differentiate on Time Zero
807 TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
808 fListQAtof->Add(hnSigT0Fill);
809 TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
810 fListQAtof->Add(hnSigT0T0);
811 TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
812 fListQAtof->Add(hnSigT0TOF);
813 TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
814 fListQAtof->Add(hnSigT0Best);
815
816 TH2F *hSig = new TH2F("hSigP_TOF",
817 "TOF signal vs. p;p [GeV]; TOF signal [ns]",
818 vX->GetNrows()-1,vX->GetMatrixArray(),
819 300,0,30);
820
821 delete vX;
822
823 fListQAtof->Add(hSig);
824
825 TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
826 fListQAtof->Add(hStartTimeMaskTOF);
827 TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
828 fListQAtof->Add(hStartTimeResTOF);
829
830 TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",20,0,20);
831 fListQAtof->Add(hnTracksAtTOF);
832 TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",20,0,20);
833 fListQAtof->Add(hT0MakerEff);
834
835 // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
836 TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
837 fListQAtof->Add(hStartTimeAT0);
838 TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
839 fListQAtof->Add(hStartTimeCT0);
840 TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
841 fListQAtof->Add(hStartTimeACT0);
842}
843
844//______________________________________________________________________________
845void AliAnalysisTaskPIDqa::SetupEMCALqa()
846{
847 //
848 // Create the EMCAL qa objects
849 //
850
851 TVectorD *vX=MakeLogBinning(200,.1,30);
852
853 TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
854 Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
855 vX->GetNrows()-1,vX->GetMatrixArray(),
856 200,-10,10);
857 fListQAemcal->Add(hNsigmaPt);
858
859 TH2F *hSigPt = new TH2F("hSigPt_EMCAL",
860 "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
861 vX->GetNrows()-1,vX->GetMatrixArray(),
862 200,0,2);
863 fListQAemcal->Add(hSigPt);
864
865 delete vX;
866}
867
868//______________________________________________________________________________
869void AliAnalysisTaskPIDqa::SetupHMPIDqa()
870{
871 //
872 // Create the HMPID qa objects
873 //
874
875 TH2F *hCkovAnglevsMom = new TH2F("hCkovAnglevsMom", "Cherenkov angle vs momnetum",500,0,5.,500,0,1);
876 fListQAhmpid->Add(hCkovAnglevsMom);
877
878}
879
880//______________________________________________________________________________
881void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
882{
883 //
884 // Create the HMPID qa objects
885 //
886
887 TH2F *hCkovAnglevsMomPion = new TH2F("hCkovAnglevsMom_pion", "Cherenkov angle vs momnetum for pions",500,0,5.,500,0,1);
888 fListQAtofhmpid->Add(hCkovAnglevsMomPion);
889
890 TH2F *hCkovAnglevsMomKaon = new TH2F("hCkovAnglevsMom_kaon", "Cherenkov angle vs momnetum for kaons",500,0,5.,500,0,1);
891 fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
892
893 TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momnetum for protons",500,0,5.,500,0,1);
894 fListQAtofhmpid->Add(hCkovAnglevsMomProton);
895
896
897}
898
899//______________________________________________________________________________
900void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
901{
902 //
903 // Create the qa objects for TPC + TOF combination
904 //
905
906 TVectorD *vX=MakeLogBinning(200,.1,30);
907
908 //TPC signals after TOF cut
909 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
910 TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
911 Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
912 vX->GetNrows()-1,vX->GetMatrixArray(),
913 200,-10,10);
914 fListQAtpctof->Add(hNsigmaP);
915 }
916
917 //TOF signals after TPC cut
918 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
919 TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
920 Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
921 vX->GetNrows()-1,vX->GetMatrixArray(),
922 200,-10,10);
923 fListQAtpctof->Add(hNsigmaP);
924 }
925
926 //EMCAL signal after TOF and TPC cut
927 for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
928 TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
929 Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
930 vX->GetNrows()-1,vX->GetMatrixArray(),
931 200,0,2);
932 fListQAtpctof->Add(heopPt);
933 }
934
935 delete vX;
936}
937
938//______________________________________________________________________________
939TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
940{
941 //
942 // Make logarithmic binning
943 // the user has to delete the array afterwards!!!
944 //
945
946 //check limits
947 if (xmin<1e-20 || xmax<1e-20){
948 AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
949 return MakeLinBinning(nbinsX, xmin, xmax);
950 }
951 if (xmax<xmin){
952 Double_t tmp=xmin;
953 xmin=xmax;
954 xmax=tmp;
955 }
956 TVectorD *binLim=new TVectorD(nbinsX+1);
957 Double_t first=xmin;
958 Double_t last=xmax;
959 Double_t expMax=TMath::Log(last/first);
960 for (Int_t i=0; i<nbinsX+1; ++i){
961 (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
962 }
963 return binLim;
964}
965
966//______________________________________________________________________________
967TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
968{
969 //
970 // Make linear binning
971 // the user has to delete the array afterwards!!!
972 //
973 if (xmax<xmin){
974 Double_t tmp=xmin;
975 xmin=xmax;
976 xmax=tmp;
977 }
978 TVectorD *binLim=new TVectorD(nbinsX+1);
979 Double_t first=xmin;
980 Double_t last=xmax;
981 Double_t binWidth=(last-first)/nbinsX;
982 for (Int_t i=0; i<nbinsX+1; ++i){
983 (*binLim)[i]=first+binWidth*(Double_t)i;
984 }
985 return binLim;
986}
987
988//_____________________________________________________________________________
989TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
990{
991 //
992 // Make arbitrary binning, bins separated by a ','
993 //
994 TString limits(bins);
995 if (limits.IsNull()){
996 AliError("Bin Limit string is empty, cannot add the variable");
997 return 0x0;
998 }
999
1000 TObjArray *arr=limits.Tokenize(",");
1001 Int_t nLimits=arr->GetEntries();
1002 if (nLimits<2){
1003 AliError("Need at leas 2 bin limits, cannot add the variable");
1004 delete arr;
1005 return 0x0;
1006 }
1007
1008 TVectorD *binLimits=new TVectorD(nLimits);
1009 for (Int_t iLim=0; iLim<nLimits; ++iLim){
1010 (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
1011 }
1012
1013 delete arr;
1014 return binLimits;
1015}
1016