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