Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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 #include <AliAODEvent.h>
43 #include <AliESDv0.h>
44 #include <AliAODv0.h>
45 #include <AliESDv0KineCuts.h>
46
47 #include "AliAnalysisTaskPIDqa.h"
48
49
50 ClassImp(AliAnalysisTaskPIDqa)
51
52 //______________________________________________________________________________
53 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
54 AliAnalysisTaskSE(),
55 fPIDResponse(0x0),
56 fV0cuts(0x0),
57 fV0electrons(0x0),
58 fV0pions(0x0),
59 fV0kaons(0x0),
60 fV0protons(0x0),
61 fListQA(0x0),
62 fListQAits(0x0),
63 fListQAitsSA(0x0),
64 fListQAitsPureSA(0x0),
65 fListQAtpc(0x0),
66 fListQAtrd(0x0),
67 fListQAtof(0x0),
68 fListQAt0(0x0),
69 fListQAemcal(0x0),
70 fListQAhmpid(0x0),
71 fListQAtofhmpid(0x0),
72 fListQAtpctof(0x0),
73 fListQAV0(0x0),
74 fListQAinfo(0x0)
75 {
76   //
77   // Dummy constructor
78   //
79 }
80
81 //______________________________________________________________________________
82 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
83 AliAnalysisTaskSE(name),
84 fPIDResponse(0x0),
85 fV0cuts(0x0),
86 fV0electrons(0x0),
87 fV0pions(0x0),
88 fV0kaons(0x0),
89 fV0protons(0x0),
90 fListQA(0x0),
91 fListQAits(0x0),
92 fListQAitsSA(0x0),
93 fListQAitsPureSA(0x0),
94 fListQAtpc(0x0),
95 fListQAtrd(0x0),
96 fListQAtof(0x0),
97 fListQAt0(0x0),
98 fListQAemcal(0x0),
99 fListQAhmpid(0x0),
100 fListQAtofhmpid(0x0),
101 fListQAtpctof(0x0),
102 fListQAV0(0x0),
103 fListQAinfo(0x0)
104 {
105   //
106   // Default constructor
107   //
108   DefineInput(0,TChain::Class());
109   DefineOutput(1,TList::Class());
110 }
111
112 //______________________________________________________________________________
113 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
114 {
115   //
116   // Destructor
117   //
118
119   delete fV0cuts;
120   delete fV0electrons;
121   delete fV0pions;
122   delete fV0kaons;
123   delete fV0protons;
124
125   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA;
126 }
127
128 //______________________________________________________________________________
129 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
130 {
131   //
132   // Create the output QA objects
133   //
134
135   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
136
137   //input hander
138   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
139   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
140   if (!inputHandler) AliFatal("Input handler needed");
141
142   //pid response object
143   fPIDResponse=inputHandler->GetPIDResponse();
144   if (!fPIDResponse) AliError("PIDResponse object was not created");
145   
146   // V0 Kine cuts 
147   fV0cuts = new AliESDv0KineCuts;
148  
149   // V0 PID Obj arrays
150   fV0electrons = new TObjArray;
151   fV0pions     = new TObjArray;
152   fV0kaons     = new TObjArray;
153   fV0protons   = new TObjArray;
154
155   //
156   fListQA=new TList;
157   fListQA->SetOwner();
158   
159   fListQAits=new TList;
160   fListQAits->SetOwner();
161   fListQAits->SetName("ITS");
162
163   fListQAitsSA=new TList;
164   fListQAitsSA->SetOwner();
165   fListQAitsSA->SetName("ITS_SA");
166
167   fListQAitsPureSA=new TList;
168   fListQAitsPureSA->SetOwner();
169   fListQAitsPureSA->SetName("ITS_PureSA");
170
171   fListQAtpc=new TList;
172   fListQAtpc->SetOwner();
173   fListQAtpc->SetName("TPC");
174   
175   fListQAtrd=new TList;
176   fListQAtrd->SetOwner();
177   fListQAtrd->SetName("TRD");
178   
179   fListQAtof=new TList;
180   fListQAtof->SetOwner();
181   fListQAtof->SetName("TOF");
182
183   fListQAt0=new TList;
184   fListQAt0->SetOwner();
185   fListQAt0->SetName("T0");
186   
187   fListQAemcal=new TList;
188   fListQAemcal->SetOwner();
189   fListQAemcal->SetName("EMCAL");
190   
191   fListQAhmpid=new TList;
192   fListQAhmpid->SetOwner();
193   fListQAhmpid->SetName("HMPID");
194   
195   fListQAtpctof=new TList;
196   fListQAtpctof->SetOwner();
197   fListQAtpctof->SetName("TPC_TOF");
198
199   fListQAtofhmpid=new TList;
200   fListQAtofhmpid->SetOwner();
201   fListQAtofhmpid->SetName("TOF_HMPID");
202   
203   fListQAV0=new TList;
204   fListQAV0->SetOwner();
205   fListQAV0->SetName("V0decay");
206
207   fListQAinfo=new TList;
208   fListQAinfo->SetOwner();
209   fListQAinfo->SetName("QAinfo");
210   
211   fListQA->Add(fListQAits);
212   fListQA->Add(fListQAitsSA);
213   fListQA->Add(fListQAitsPureSA);
214   fListQA->Add(fListQAtpc);
215   fListQA->Add(fListQAtrd);
216   fListQA->Add(fListQAtof);
217   fListQA->Add(fListQAt0);
218   fListQA->Add(fListQAemcal);
219   fListQA->Add(fListQAhmpid);
220   fListQA->Add(fListQAtpctof);
221   fListQA->Add(fListQAtofhmpid);
222   fListQA->Add(fListQAV0);
223   fListQA->Add(fListQAinfo);
224
225   SetupITSqa();
226   SetupTPCqa();
227   SetupTRDqa();
228   SetupTOFqa();
229   SetupT0qa();
230   SetupEMCALqa();
231   SetupHMPIDqa();
232   SetupTPCTOFqa();
233   SetupTOFHMPIDqa();
234   SetupV0qa();
235   SetupQAinfo();
236   
237   PostData(1,fListQA);
238 }
239
240
241 //______________________________________________________________________________
242 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
243 {
244   //
245   // Setup the PID response functions and fill the QA histograms
246   //
247
248   AliVEvent *event=InputEvent();
249   if (!event||!fPIDResponse) return;
250
251   // Start with the V0 task (only possible for ESDs?)
252   FillV0PIDlist();
253   
254   FillITSqa();
255   FillTPCqa();
256   FillTRDqa();
257   FillTOFqa();
258   FillEMCALqa();
259   FillHMPIDqa();
260   FillT0qa();
261   
262   //combined detector QA
263   FillTPCTOFqa();
264   FillTOFHMPIDqa();
265   
266   // Clear the V0 PID arrays
267   ClearV0PIDlist();
268
269   //QA info
270   FillQAinfo();
271   
272   PostData(1,fListQA);
273 }
274
275 //______________________________________________________________________________
276 void  AliAnalysisTaskPIDqa::FillV0PIDlist(){
277
278   //
279   // Fill the PID object arrays holding the pointers to identified particle tracks
280   //
281
282   // Dynamic cast to ESD events (DO NOTHING for AOD events)
283   AliESDEvent *event = dynamic_cast<AliESDEvent *>(InputEvent());
284   if ( !event )  return;
285   
286   if(TString(event->GetBeamType())=="Pb-Pb" || TString(event->GetBeamType())=="A-A"){
287     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb); 
288   }
289   else{
290     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP); 
291   }
292
293   // V0 selection
294   // set event
295   fV0cuts->SetEvent(event);
296
297   // loop over V0 particles
298   for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
299
300     AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0);
301  
302     if(!v0) continue;
303     if(v0->GetOnFlyStatus()) continue; 
304   
305     // Get the particle selection 
306     Bool_t foundV0 = kFALSE;
307     Int_t pdgV0, pdgP, pdgN;
308
309     foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
310     if(!foundV0) continue;
311     
312     Int_t iTrackP = v0->GetPindex();  // positive track
313     Int_t iTrackN = v0->GetNindex();  // negative track
314
315     // v0 Armenteros plot (QA)
316     Float_t armVar[2] = {0.0,0.0};
317     fV0cuts->Armenteros(v0, armVar);
318
319     TH2 *h=(TH2*)fListQAV0->At(0);
320     if (!h) continue;
321     h->Fill(armVar[0],armVar[1]);
322
323     // fill the Object arrays
324     // positive particles
325     if( pdgP == -11){
326       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
327     }
328     else if( pdgP == 211){
329       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
330     }
331     else if( pdgP == 321){
332       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackP));
333     }
334     else if( pdgP == 2212){
335       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
336     }
337
338     // negative particles
339     if( pdgN == 11){
340       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
341     }
342     else if( pdgN == -211){
343       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
344     }
345     else if( pdgN == -321){
346       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackN));
347     }
348     else if( pdgN == -2212){
349       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
350     }
351   
352
353   }
354 }
355 //______________________________________________________________________________
356 void  AliAnalysisTaskPIDqa::ClearV0PIDlist(){
357
358   //
359   // Clear the PID object arrays
360   //
361
362   fV0electrons->Clear();
363   fV0pions->Clear();
364   fV0kaons->Clear();
365   fV0protons->Clear();
366
367 }
368 //______________________________________________________________________________
369 void AliAnalysisTaskPIDqa::FillITSqa()
370 {
371   //
372   // Fill PID qa histograms for the ITS
373   //
374
375   AliVEvent *event=InputEvent();
376   
377   Int_t ntracks=event->GetNumberOfTracks();
378   for(Int_t itrack = 0; itrack < ntracks; itrack++){
379     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
380     ULong_t status=track->GetStatus();
381     // not that nice. status bits not in virtual interface
382     // ITS refit + ITS pid selection
383     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
384            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
385     Double_t mom=track->P();
386     
387     TList *theList = 0x0;
388     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
389       //ITS+TPC tracks
390       theList=fListQAits;
391     }else{
392       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){ 
393         //ITS Standalone tracks
394         theList=fListQAitsSA;
395       }else{
396         //ITS Pure Standalone tracks
397         theList=fListQAitsPureSA;
398       }
399     }
400     
401     
402     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
403       TH2 *h=(TH2*)theList->At(ispecie);
404       if (!h) continue;
405       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
406       h->Fill(mom,nSigma);
407     }
408     TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC);
409     if (h) {
410       Double_t sig=track->GetITSsignal();
411       h->Fill(mom,sig);
412     }
413   }
414 }
415
416 //______________________________________________________________________________
417 void AliAnalysisTaskPIDqa::FillTPCqa()
418 {
419   //
420   // Fill PID qa histograms for the TPC
421   //
422   
423   AliVEvent *event=InputEvent();
424   
425   Int_t ntracks=event->GetNumberOfTracks();
426   for(Int_t itrack = 0; itrack < ntracks; itrack++){
427     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
428     
429     //
430     //basic track cuts
431     //
432     ULong_t status=track->GetStatus();
433     // not that nice. status bits not in virtual interface
434     // TPC refit + ITS refit + TPC pid
435     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
436         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
437
438     // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
439     //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
440     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
441     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
442     if (track->GetTPCNclsF()>0) {
443       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
444     }
445     
446     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
447     
448     Double_t mom=track->GetTPCmomentum();
449     // the default scenario
450     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
451       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
452       if (!h) continue;
453       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
454       h->Fill(mom,nSigma);
455     }
456     // the "hybrid" scenario
457     if (track->GetTPCdEdxInfo()){
458       for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
459         TH2 *h=(TH2*)fListQAtpc->At(ispecie+AliPID::kSPECIESC);
460         if (!h) continue;
461         Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxHybrid);
462         h->Fill(mom,nSigma);
463       }
464
465       // the "OROC" scenario
466       for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
467         TH2 *h=(TH2*)fListQAtpc->At(ispecie+2*AliPID::kSPECIESC);
468         if (!h) continue;
469         Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxOROC);
470         //TSpline3* spline = fPIDResponse->GetTPCResponse().GetCurrentResponseFunction();
471         //std::cout<<ispecie<<" "<<nSigma<<" phi:"<<track->Phi()<<". "<<std::endl;
472         //if (spline) {cout<<spline->GetName()<<endl;}
473         //else {cout<<"NULL spline"<<endl;}
474         h->Fill(mom,nSigma);
475         }
476     }
477     
478     TH2 *h=(TH2*)fListQAtpc->At(3*AliPID::kSPECIESC);
479
480     if (h) {
481       Double_t sig=track->GetTPCsignal();
482       h->Fill(mom,sig);
483     }
484   }
485 }
486
487 //______________________________________________________________________________
488 void AliAnalysisTaskPIDqa::FillTRDqa()
489 {
490   //
491   // Fill PID qa histograms for the TRD
492   //
493   AliVEvent *event=InputEvent();
494   Int_t ntracks = event->GetNumberOfTracks();
495   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
496     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
497
498     //
499     //basic track cuts
500     //
501     ULong_t status=track->GetStatus();
502     // not that nice. status bits not in virtual interface
503     // TPC refit + ITS refit + TPC pid + TRD out
504     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
505         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
506 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) || //removes light nuclei. So it is out for the moment
507         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
508     
509     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
510     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
511     if (track->GetTPCNclsF()>0) {
512       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
513     }
514     
515     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
516
517     Double_t likelihoods[AliPID::kSPECIES];
518     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
519     Int_t ntracklets = 0;
520     Double_t momentum = -1.;
521     for(Int_t itl = 0; itl < 6; itl++)
522       if(track->GetTRDmomentum(itl) > 0.){
523         ntracklets++;
524         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
525     } 
526     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
527       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
528       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
529     }
530   }
531 }
532
533 //______________________________________________________________________________
534 void AliAnalysisTaskPIDqa::FillTOFqa()
535 {
536   //
537   // Fill TOF information
538   //
539   AliVEvent *event=InputEvent();
540
541   Int_t ntracks=event->GetNumberOfTracks();
542   Int_t tracksAtTof = 0;
543   for(Int_t itrack = 0; itrack < ntracks; itrack++){
544     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
545
546     //
547     //basic track cuts
548     //
549     ULong_t status=track->GetStatus();
550     // TPC refit + ITS refit +
551     // TOF out + kTIME
552     // kTIME
553     // (we don't use kTOFmismatch because it depends on TPC and kTOFpid because it prevents light nuclei
554     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
555         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
556         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
557         //        !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
558         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
559
560     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
561     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
562     if (track->GetTPCNclsF()>0) {
563       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
564     }
565
566     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
567
568     tracksAtTof++;
569
570     Double_t mom=track->P();
571
572     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
573       TH2 *h=(TH2*)fListQAtof->At(ispecie);
574       if (!h) continue;
575       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
576       h->Fill(mom,nSigma);
577     }
578
579     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
580     if (h) {
581       Double_t sig=track->GetTOFsignal()/1000.;
582       h->Fill(mom,sig);
583     }
584
585     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
586     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
587
588     if (mom >= 0.75 && mom <= 1.25 ) {
589       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
590       if (mask == 0) {
591         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
592       } else if (mask == 1) {
593         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
594       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
595         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
596       } else {
597         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
598       }
599       if (mask & 0x1) { //at least TOF-T0 present
600         Double_t delta=0;
601         (void)fPIDResponse->GetSignalDelta((AliPIDResponse::EDetector)AliPIDResponse::kTOF,track,(AliPID::EParticleType)AliPID::kPion,delta);
602         ((TH1F*)fListQAtof->FindObject("hDelta_TOF_Pion"))->Fill(delta);
603       }
604     }
605
606     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
607     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
608
609     Double_t startTimeT0 = event->GetT0TOF(0);
610     if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
611     else {
612       startTimeT0 = event->GetT0TOF(1);
613       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
614       startTimeT0 = event->GetT0TOF(2);
615       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
616     }
617   }
618   if (tracksAtTof > 0) {
619     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
620     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
621     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
622   }
623 }
624
625 //______________________________________________________________________________
626 void AliAnalysisTaskPIDqa::FillT0qa()
627 {
628   //
629   // Fill TOF information
630   //
631   AliVEvent *event=InputEvent();
632
633   Int_t ntracks=event->GetNumberOfTracks();
634
635   Int_t tracksAtT0 = 0;
636
637   for(Int_t itrack = 0; itrack < ntracks; itrack++){
638     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
639
640     //
641     //basic track cuts
642     //
643     ULong_t status=track->GetStatus();
644     // TPC refit + ITS refit +
645     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
646         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
647     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
648     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
649     if (track->GetTPCNclsF()>0) {
650       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
651     }
652     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
653
654     tracksAtT0++;
655   }
656
657   Bool_t t0A = kFALSE;
658   Bool_t t0C = kFALSE;
659   Bool_t t0And = kFALSE;
660   Double_t startTimeT0 = event->GetT0TOF(0);     // AND
661   if (startTimeT0 < 90000) {
662     t0And = kTRUE;
663     ((TH1F*)fListQAt0->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
664     }
665   startTimeT0 = event->GetT0TOF(1);             // T0A 
666   if (startTimeT0 < 90000) {
667     t0A = kTRUE;
668     ((TH1F*)fListQAt0->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
669     
670   }
671   startTimeT0 = event->GetT0TOF(2);             // T0C 
672   if (startTimeT0 < 90000) {
673     t0C = kTRUE;
674     ((TH1F*)fListQAt0->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
675   }
676   
677   ((TH1F* )fListQAt0->FindObject("hnTracksAt_T0"))->Fill(tracksAtT0);
678   if (t0A) ((TH1F*)fListQAt0->FindObject("hT0AEff"))->Fill(tracksAtT0);
679   if (t0C) ((TH1F*)fListQAt0->FindObject("hT0CEff"))->Fill(tracksAtT0);
680   if (t0And) ((TH1F*)fListQAt0->FindObject("hT0AndEff"))->Fill(tracksAtT0);
681   if (t0A || t0C) ((TH1F*)fListQAt0->FindObject("hT0OrEff"))->Fill(tracksAtT0);
682 }
683
684
685 //______________________________________________________________________________
686 void AliAnalysisTaskPIDqa::FillEMCALqa()
687 {
688   //
689   // Fill PID qa histograms for the EMCAL
690   //
691
692   AliVEvent *event=InputEvent();
693   
694   Int_t ntracks=event->GetNumberOfTracks();
695   for(Int_t itrack = 0; itrack < ntracks; itrack++){
696     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
697     
698     //
699     //basic track cuts
700     //
701     ULong_t status=track->GetStatus();
702     // not that nice. status bits not in virtual interface
703     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
704
705     Double_t pt=track->Pt();
706    
707     //EMCAL nSigma (only for electrons at the moment)
708     TH2 *h=(TH2*)fListQAemcal->At(0);
709     if (!h) continue;
710     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
711     h->Fill(pt,nSigma);
712     
713   }
714
715    //EMCAL signal (E/p vs. pT) for electrons from V0
716   for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
717     AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
718
719     //
720     //basic track cuts
721     //
722     ULong_t status=track->GetStatus();
723     // not that nice. status bits not in virtual interface
724     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
725
726     Double_t pt=track->Pt();
727
728     TH2 *h=(TH2*)fListQAemcal->At(1);
729     if (h) {
730
731       Int_t nMatchClus = track->GetEMCALcluster();
732       Double_t mom     = track->P();
733       Double_t eop     = -1.;
734
735       if(nMatchClus > -1){
736     
737         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
738
739         if(matchedClus){
740
741           // matched cluster is EMCAL
742           if(matchedClus->IsEMCAL()){
743
744             Double_t fClsE       = matchedClus->E();
745             eop                  = fClsE/mom;
746
747             h->Fill(pt,eop);
748
749           }
750         }
751       }
752     }
753   }
754
755    //EMCAL signal (E/p vs. pT) for pions from V0
756   for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
757     AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
758
759     //
760     //basic track cuts
761     //
762     ULong_t status=track->GetStatus();
763     // not that nice. status bits not in virtual interface
764     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
765
766     Double_t pt=track->Pt();
767
768     TH2 *h=(TH2*)fListQAemcal->At(2);
769     if (h) {
770
771       Int_t nMatchClus = track->GetEMCALcluster();
772       Double_t mom     = track->P();
773       Double_t eop     = -1.;
774
775       if(nMatchClus > -1){
776     
777         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
778
779         if(matchedClus){
780
781           // matched cluster is EMCAL
782           if(matchedClus->IsEMCAL()){
783
784             Double_t fClsE       = matchedClus->E();
785             eop                  = fClsE/mom;
786
787             h->Fill(pt,eop);
788
789           }
790         }
791       }
792     }
793   }
794
795    //EMCAL signal (E/p vs. pT) for protons from V0
796   for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
797     AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
798
799     //
800     //basic track cuts
801     //
802     ULong_t status=track->GetStatus();
803     // not that nice. status bits not in virtual interface
804     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
805
806     Double_t pt=track->Pt();
807
808     TH2 *hP=(TH2*)fListQAemcal->At(3);
809     TH2 *hAP=(TH2*)fListQAemcal->At(4);
810     if (hP && hAP) {
811
812       Int_t nMatchClus = track->GetEMCALcluster();
813       Double_t mom     = track->P();
814       Int_t charge     = track->Charge();             
815       Double_t eop     = -1.;
816
817       if(nMatchClus > -1){
818     
819         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
820
821         if(matchedClus){
822
823           // matched cluster is EMCAL
824           if(matchedClus->IsEMCAL()){
825
826             Double_t fClsE       = matchedClus->E();
827             eop                  = fClsE/mom;
828
829             if(charge > 0)      hP->Fill(pt,eop);
830             else if(charge < 0) hAP->Fill(pt,eop);
831
832           }
833         }
834       }
835     }
836   }
837
838 }
839
840
841 //______________________________________________________________________________
842 void AliAnalysisTaskPIDqa::FillHMPIDqa()
843 {
844   //
845   // Fill PID qa histograms for the HMPID
846   //
847   
848   AliVEvent *event=InputEvent();
849   
850   Int_t ntracks=event->GetNumberOfTracks();
851   for(Int_t itrack = 0; itrack < ntracks; itrack++){
852     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
853     
854     //
855     //basic track cuts
856     //
857     const ULong_t status=track->GetStatus();
858     // not that nice. status bits not in virtual interface
859     // TPC refit + ITS refit +
860     // TOF out + TOFpid +
861     // kTIME
862     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
863         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
864
865     const Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
866     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
867     if (track->GetTPCNclsF()>0) {
868       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
869     }
870
871     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
872     
873     const Double_t mom = track->P();
874     const Double_t ckovAngle = track->GetHMPIDsignal();
875
876     Int_t nhists=0;
877     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
878       if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
879       TH2 *h=(TH2*)fListQAhmpid->At(nhists);
880       if (!h) {++nhists; continue;}
881       const Double_t nSigma=fPIDResponse->NumberOfSigmasHMPID(track, (AliPID::EParticleType)ispecie);
882       h->Fill(mom,nSigma);
883       ++nhists;
884     }
885     
886     TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(AliPID::kSPECIESC);
887     
888     if (hThetavsMom) hThetavsMom->Fill(mom,ckovAngle);
889   
890   }
891 }
892 //______________________________________________________________________________
893 void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
894 {
895   //
896   // Fill PID qa histograms for the HMPID
897   //
898   
899   AliVEvent *event=InputEvent();
900   
901   Int_t ntracks=event->GetNumberOfTracks();
902   for(Int_t itrack = 0; itrack < ntracks; itrack++){
903     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
904     
905     //
906     //basic track cuts
907     //
908     ULong_t status=track->GetStatus();
909     // not that nice. status bits not in virtual interface
910     // TPC refit + ITS refit +
911     // TOF out + TOFpid +
912     // kTIME
913     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
914         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
915         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
916         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
917         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
918
919     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
920     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
921     if (track->GetTPCNclsF()>0) {
922       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
923     }
924
925     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
926     
927     Double_t mom = track->P();
928     Double_t ckovAngle = track->GetHMPIDsignal();
929     
930     Double_t nSigmaTOF[3]; 
931     TH1F *h[3];
932     
933     for (Int_t ispecie=2; ispecie<5; ++ispecie){
934       //TOF nSigma
935       nSigmaTOF[ispecie-2]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
936       h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
937       
938     if(TMath::Abs(nSigmaTOF[0])<2)                                                              h[0]->Fill(mom,ckovAngle);
939     
940     if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3)                                h[1]->Fill(mom,ckovAngle);
941
942     if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3)  h[2]->Fill(mom,ckovAngle);
943       
944   }
945   
946 }
947
948 //______________________________________________________________________________
949 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
950 {
951   //
952   // Fill PID qa histograms for the TOF
953   //   Here also the TPC histograms after TOF selection are filled
954   //
955
956   AliVEvent *event=InputEvent();
957
958   Int_t ntracks=event->GetNumberOfTracks();
959   for(Int_t itrack = 0; itrack < ntracks; itrack++){
960     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
961
962     //
963     //basic track cuts
964     //
965     ULong_t status=track->GetStatus();
966     // not that nice. status bits not in virtual interface
967     // TPC refit + ITS refit +
968     // TOF out + TOFpid +
969     // kTIME
970     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
971         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
972 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
973         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
974         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
975         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
976
977     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
978     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
979     if (track->GetTPCNclsF()>0) {
980       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
981     }
982
983     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
984
985
986     Double_t mom=track->P();
987     Double_t momTPC=track->GetTPCmomentum();
988
989     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
990       //TOF nSigma
991       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
992       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
993
994       //TPC after TOF cut
995       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
996       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
997
998       //TOF after TPC cut
999       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
1000       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
1001
1002       //EMCAL after TOF and TPC cut
1003       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
1004       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
1005
1006         Int_t nMatchClus = track->GetEMCALcluster();
1007         Double_t pt      = track->Pt();
1008         Double_t eop     = -1.;
1009         
1010         if(nMatchClus > -1){
1011           
1012           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1013           
1014           if(matchedClus){
1015             
1016             // matched cluster is EMCAL
1017             if(matchedClus->IsEMCAL()){
1018               
1019               Double_t fClsE       = matchedClus->E();
1020               eop                  = fClsE/mom;
1021
1022               h->Fill(pt,eop);
1023  
1024               
1025             }
1026           }
1027         }
1028       }
1029     }
1030   }
1031 }
1032
1033 //_____________________________________________________________________________
1034 void AliAnalysisTaskPIDqa::FillQAinfo()
1035 {
1036   //
1037   // Fill the QA information
1038   //
1039
1040
1041   //TPC QA info
1042   TObjArray *arrTPC=static_cast<TObjArray*>(fListQAinfo->At(0));
1043   if (fPIDResponse && arrTPC){
1044     AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
1045     // fill spline names
1046     if (!arrTPC->UncheckedAt(0)){
1047       
1048       TObjArray *arrTPCsplineNames=new TObjArray(AliPID::kSPECIESC);
1049       arrTPCsplineNames->SetOwner();
1050       arrTPCsplineNames->SetName("TPC_spline_names");
1051       arrTPC->AddAt(arrTPCsplineNames,0);
1052       
1053       for (Int_t iresp=0; iresp<AliPID::kSPECIESC; ++iresp){
1054         const TObject *o=tpcResp.GetResponseFunction((AliPID::EParticleType)iresp);
1055         if (!o) continue;
1056         arrTPCsplineNames->Add(new TObjString(Form("%02d: %s",iresp, o->GetName())));
1057       }
1058     }
1059
1060     // tpc response config
1061     if (!arrTPC->UncheckedAt(1)){
1062       
1063       TObjArray *arrTPCconfigInfo=new TObjArray;
1064       arrTPCconfigInfo->SetOwner();
1065       arrTPCconfigInfo->SetName("TPC_config_info");
1066       arrTPC->AddAt(arrTPCconfigInfo,1);
1067
1068       TObjString *ostr=0x0;
1069       ostr=new TObjString;
1070       ostr->String().Form("Eta Corr map: %s", tpcResp.GetEtaCorrMap()?tpcResp.GetEtaCorrMap()->GetName():"none");
1071       arrTPCconfigInfo->Add(ostr);
1072
1073       ostr=new TObjString;
1074       ostr->String().Form("Sigma Par map: %s", tpcResp.GetSigmaPar1Map()?tpcResp.GetSigmaPar1Map()->GetName():"none");
1075       arrTPCconfigInfo->Add(ostr);
1076
1077       ostr=new TObjString;
1078       ostr->String().Form("MIP: %.2f", tpcResp.GetMIP());
1079       arrTPCconfigInfo->Add(ostr);
1080       
1081       ostr=new TObjString;
1082       ostr->String().Form("Res: Def %.3g (%.3g) : AllHigh %.3g (%.3g) : OROC high %.3g (%.3g)",
1083                           tpcResp.GetRes0(AliTPCPIDResponse::kDefault), tpcResp.GetResN2(AliTPCPIDResponse::kDefault),
1084                           tpcResp.GetRes0(AliTPCPIDResponse::kALLhigh), tpcResp.GetResN2(AliTPCPIDResponse::kALLhigh),
1085                           tpcResp.GetRes0(AliTPCPIDResponse::kOROChigh), tpcResp.GetResN2(AliTPCPIDResponse::kOROChigh)
1086                          );
1087       arrTPCconfigInfo->Add(ostr);
1088     }
1089   }
1090 }
1091
1092 //______________________________________________________________________________
1093 void AliAnalysisTaskPIDqa::SetupITSqa()
1094 {
1095   //
1096   // Create the ITS qa objects
1097   //
1098   
1099   TVectorD *vX=MakeLogBinning(200,.1,30);
1100   
1101   //ITS+TPC tracks
1102   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1103     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
1104                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1105                               vX->GetNrows()-1,vX->GetMatrixArray(),
1106                               200,-10,10);
1107     fListQAits->Add(hNsigmaP);
1108   }
1109   TH2F *hSig = new TH2F("hSigP_ITS",
1110                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1111                         vX->GetNrows()-1,vX->GetMatrixArray(),
1112                         300,0,300);
1113   fListQAits->Add(hSig);
1114
1115   //ITS Standalone tracks
1116   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1117     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
1118                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1119                                 vX->GetNrows()-1,vX->GetMatrixArray(),
1120                                 200,-10,10);
1121     fListQAitsSA->Add(hNsigmaPSA);
1122   }
1123   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
1124                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1125                           vX->GetNrows()-1,vX->GetMatrixArray(),
1126                           300,0,300);
1127   fListQAitsSA->Add(hSigSA);
1128   
1129   //ITS Pure Standalone tracks
1130   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1131     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
1132                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1133                                     vX->GetNrows()-1,vX->GetMatrixArray(),
1134                                     200,-10,10);
1135     fListQAitsPureSA->Add(hNsigmaPPureSA);
1136   }
1137   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
1138                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1139                               vX->GetNrows()-1,vX->GetMatrixArray(),
1140                               300,0,300);
1141   fListQAitsPureSA->Add(hSigPureSA);
1142   
1143   delete vX;  
1144 }
1145
1146 //______________________________________________________________________________
1147 void AliAnalysisTaskPIDqa::SetupTPCqa()
1148 {
1149   //
1150   // Create the TPC qa objects
1151   //
1152   
1153   TVectorD *vX=MakeLogBinning(200,.1,30);
1154   
1155   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1156     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
1157                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1158                               vX->GetNrows()-1,vX->GetMatrixArray(),
1159                               200,-10,10);
1160     fListQAtpc->Add(hNsigmaP);
1161   }
1162
1163   // the "hybrid" scenario
1164   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1165     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_Hybrid",AliPID::ParticleName(ispecie)),
1166                               Form("TPC n#sigma %s vs. p (Hybrid gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1167                               vX->GetNrows()-1,vX->GetMatrixArray(),
1168                               200,-10,10);
1169     fListQAtpc->Add(hNsigmaP);
1170   }
1171    
1172   // the "OROC high" scenario
1173   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1174     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_OROChigh",AliPID::ParticleName(ispecie)),
1175                               Form("TPC n#sigma %s vs. p (OROChigh gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1176                               vX->GetNrows()-1,vX->GetMatrixArray(),
1177                               200,-10,10);
1178     fListQAtpc->Add(hNsigmaP);
1179   }
1180   
1181   
1182   
1183   TH2F *hSig = new TH2F("hSigP_TPC",
1184                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
1185                         vX->GetNrows()-1,vX->GetMatrixArray(),
1186                         500,0,1000);
1187   fListQAtpc->Add(hSig); //3*AliPID::kSPECIESC
1188
1189   delete vX;  
1190 }
1191
1192 //______________________________________________________________________________
1193 void AliAnalysisTaskPIDqa::SetupTRDqa()
1194 {
1195   //
1196   // Create the TRD qa objects
1197   //
1198   TVectorD *vX=MakeLogBinning(200,.1,30);
1199   for(Int_t itl = 0; itl < 6; ++itl){
1200     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1201       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
1202                               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)),
1203                               vX->GetNrows()-1, vX->GetMatrixArray(),
1204                               100, 0., 1.);
1205       fListQAtrd->Add(hLikeP);
1206     }
1207   }
1208   delete vX;
1209 }
1210
1211 //______________________________________________________________________________
1212 void AliAnalysisTaskPIDqa::SetupTOFqa()
1213 {
1214   //
1215   // Create the TOF qa objects
1216   //
1217   
1218   TVectorD *vX=MakeLogBinning(200,.1,30);
1219
1220   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1221     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
1222                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1223                               vX->GetNrows()-1,vX->GetMatrixArray(),
1224                               200,-10,10);
1225     fListQAtof->Add(hNsigmaP);
1226   }
1227
1228   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
1229   fListQAtof->Add(hnSigT0Fill);
1230   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
1231   fListQAtof->Add(hnSigT0T0);
1232   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
1233   fListQAtof->Add(hnSigT0TOF);
1234   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
1235   fListQAtof->Add(hnSigT0Best);
1236   TH1F *hnDeltaPi = new TH1F("hDelta_TOF_Pion","DeltaT (Pion) [0.75-1.25 GeV/c]",50,-500,500);
1237   fListQAtof->Add(hnDeltaPi);
1238   
1239   TH2F *hSig = new TH2F("hSigP_TOF",
1240                         "TOF signal vs. p;p [GeV]; TOF signal [ns]",
1241                         vX->GetNrows()-1,vX->GetMatrixArray(),
1242                         300,0,30);
1243
1244   delete vX;
1245   
1246   fListQAtof->Add(hSig);
1247
1248   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
1249   fListQAtof->Add(hStartTimeMaskTOF);
1250   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
1251   fListQAtof->Add(hStartTimeResTOF);
1252
1253   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",100,0,100);
1254   fListQAtof->Add(hnTracksAtTOF);
1255   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",100,0,100);
1256   fListQAtof->Add(hT0MakerEff);
1257
1258   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
1259   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
1260   fListQAtof->Add(hStartTimeAT0);
1261   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
1262   fListQAtof->Add(hStartTimeCT0);
1263   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
1264   fListQAtof->Add(hStartTimeACT0);
1265 }
1266
1267
1268 //______________________________________________________________________________
1269 void AliAnalysisTaskPIDqa::SetupT0qa()
1270 {
1271   //
1272   // Create the T0 qa objects
1273   //
1274   
1275   // these are similar to plots inside TOFqa, but these are for all events
1276   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
1277   fListQAt0->Add(hStartTimeAT0);
1278   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
1279   fListQAt0->Add(hStartTimeCT0);
1280   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
1281   fListQAt0->Add(hStartTimeACT0);
1282
1283   TH1F *hnTracksAtT0 = new TH1F("hnTracksAt_T0","Tracks for events selected for T0",100,0,100);
1284   fListQAt0->Add(hnTracksAtT0);
1285   TH1F *hT0AEff = new TH1F("hT0AEff","Events with T0A vs nTracks",100,0,100);
1286   fListQAt0->Add(hT0AEff);
1287   TH1F *hT0CEff = new TH1F("hT0CEff","Events with T0C vs nTracks",100,0,100);
1288   fListQAt0->Add(hT0CEff);
1289   TH1F *hT0AndEff = new TH1F("hT0AndEff","Events with T0AC (AND) vs nTracks",100,0,100);
1290   fListQAt0->Add(hT0AndEff);
1291   TH1F *hT0OrEff = new TH1F("hT0OrEff","Events with T0AC (OR) vs nTracks",100,0,100);
1292   fListQAt0->Add(hT0OrEff);
1293
1294
1295 }
1296
1297 //______________________________________________________________________________
1298 void AliAnalysisTaskPIDqa::SetupEMCALqa()
1299 {
1300   //
1301   // Create the EMCAL qa objects
1302   //
1303
1304   TVectorD *vX=MakeLogBinning(200,.1,30);
1305   
1306   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
1307                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
1308                              vX->GetNrows()-1,vX->GetMatrixArray(),
1309                              200,-10,10);
1310   fListQAemcal->Add(hNsigmaPt);  
1311   
1312
1313   TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele",
1314                         "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1315                         vX->GetNrows()-1,vX->GetMatrixArray(),
1316                         200,0,2);
1317   fListQAemcal->Add(hSigPtEle);
1318
1319   TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions",
1320                         "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1321                         vX->GetNrows()-1,vX->GetMatrixArray(),
1322                         200,0,2);
1323   fListQAemcal->Add(hSigPtPions);
1324
1325   TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons",
1326                         "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1327                         vX->GetNrows()-1,vX->GetMatrixArray(),
1328                         200,0,2);
1329   fListQAemcal->Add(hSigPtProtons);
1330
1331   TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons",
1332                         "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1333                         vX->GetNrows()-1,vX->GetMatrixArray(),
1334                         200,0,2);
1335   fListQAemcal->Add(hSigPtAntiProtons);
1336
1337   delete vX;  
1338 }
1339
1340 //______________________________________________________________________________
1341 void AliAnalysisTaskPIDqa::SetupHMPIDqa()
1342 {
1343   //
1344   // Create the HMPID qa objects
1345   //
1346
1347   TVectorD *vX=MakeLogBinning(200,.1,30);
1348
1349   // nSigmas
1350   Int_t nhists=0;
1351   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
1352     if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
1353     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_HMPID_%s",AliPID::ParticleName(ispecie)),
1354                               Form("HMPID n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1355                               vX->GetNrows()-1,vX->GetMatrixArray(),
1356                               200,-10,10);
1357     fListQAhmpid->AddAt(hNsigmaP, nhists);
1358     ++nhists;
1359   }
1360   
1361   // cherenkov angle
1362   TH2F *hCkovAnglevsMom   = new TH2F("hCkovAnglevsMom",  "Cherenkov angle vs momentum",
1363                                      vX->GetNrows()-1,vX->GetMatrixArray(),
1364                                      500,0,1);
1365   fListQAhmpid->AddAt(hCkovAnglevsMom,nhists);
1366   
1367   delete vX;
1368 }
1369
1370 //______________________________________________________________________________
1371 void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
1372 {
1373   //
1374   // Create the HMPID qa objects
1375   //
1376   
1377   TH2F *hCkovAnglevsMomPion   = new TH2F("hCkovAnglevsMom_pion",  "Cherenkov angle vs momentum for pions",500,0,5.,500,0,1);
1378   fListQAtofhmpid->Add(hCkovAnglevsMomPion);
1379   
1380   TH2F *hCkovAnglevsMomKaon   = new TH2F("hCkovAnglevsMom_kaon",  "Cherenkov angle vs momentum for kaons",500,0,5.,500,0,1);
1381   fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
1382   
1383   TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momentum for protons",500,0,5.,500,0,1);
1384   fListQAtofhmpid->Add(hCkovAnglevsMomProton);
1385   
1386   
1387 }  
1388
1389 //______________________________________________________________________________
1390 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
1391 {
1392   //
1393   // Create the qa objects for TPC + TOF combination
1394   //
1395   
1396   TVectorD *vX=MakeLogBinning(200,.1,30);
1397
1398   //TPC signals after TOF cut
1399   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1400     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
1401                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1402                               vX->GetNrows()-1,vX->GetMatrixArray(),
1403                               200,-10,10);
1404     fListQAtpctof->Add(hNsigmaP);
1405   }
1406
1407   //TOF signals after TPC cut
1408   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1409     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
1410                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1411                               vX->GetNrows()-1,vX->GetMatrixArray(),
1412                               200,-10,10);
1413     fListQAtpctof->Add(hNsigmaP);
1414   }
1415
1416   //EMCAL signal after TOF and TPC cut
1417   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
1418     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
1419                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
1420                             vX->GetNrows()-1,vX->GetMatrixArray(),
1421                             200,0,2);
1422     fListQAtpctof->Add(heopPt);
1423   }
1424
1425   delete vX;
1426 }
1427 //______________________________________________________________________________
1428 void AliAnalysisTaskPIDqa::SetupV0qa()
1429 {
1430   //
1431   // Create the qa objects for V0 Kine cuts
1432   //
1433   
1434   TH2F *hArmenteros  = new TH2F("hArmenteros",  "Armenteros plot",200,-1.,1.,200,0.,0.4);
1435   fListQAV0->Add(hArmenteros);
1436  
1437 }
1438
1439 //_____________________________________________________________________________
1440 void AliAnalysisTaskPIDqa::SetupQAinfo(){
1441   //
1442   // Setup the info of QA objects
1443   //
1444
1445   TObjArray *arr=new TObjArray;
1446   arr->SetName("TPC_info");
1447   fListQAinfo->Add(arr);
1448 }
1449
1450 //______________________________________________________________________________
1451 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
1452 {
1453   //
1454   // Make logarithmic binning
1455   // the user has to delete the array afterwards!!!
1456   //
1457   
1458   //check limits
1459   if (xmin<1e-20 || xmax<1e-20){
1460     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
1461     return MakeLinBinning(nbinsX, xmin, xmax);
1462   }
1463   if (xmax<xmin){
1464     Double_t tmp=xmin;
1465     xmin=xmax;
1466     xmax=tmp;
1467   }
1468   TVectorD *binLim=new TVectorD(nbinsX+1);
1469   Double_t first=xmin;
1470   Double_t last=xmax;
1471   Double_t expMax=TMath::Log(last/first);
1472   for (Int_t i=0; i<nbinsX+1; ++i){
1473     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
1474   }
1475   return binLim;
1476 }
1477
1478 //______________________________________________________________________________
1479 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
1480 {
1481   //
1482   // Make linear binning
1483   // the user has to delete the array afterwards!!!
1484   //
1485   if (xmax<xmin){
1486     Double_t tmp=xmin;
1487     xmin=xmax;
1488     xmax=tmp;
1489   }
1490   TVectorD *binLim=new TVectorD(nbinsX+1);
1491   Double_t first=xmin;
1492   Double_t last=xmax;
1493   Double_t binWidth=(last-first)/nbinsX;
1494   for (Int_t i=0; i<nbinsX+1; ++i){
1495     (*binLim)[i]=first+binWidth*(Double_t)i;
1496   }
1497   return binLim;
1498 }
1499
1500 //_____________________________________________________________________________
1501 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
1502 {
1503   //
1504   // Make arbitrary binning, bins separated by a ','
1505   //
1506   TString limits(bins);
1507   if (limits.IsNull()){
1508     AliError("Bin Limit string is empty, cannot add the variable");
1509     return 0x0;
1510   }
1511   
1512   TObjArray *arr=limits.Tokenize(",");
1513   Int_t nLimits=arr->GetEntries();
1514   if (nLimits<2){
1515     AliError("Need at leas 2 bin limits, cannot add the variable");
1516     delete arr;
1517     return 0x0;
1518   }
1519   
1520   TVectorD *binLimits=new TVectorD(nLimits);
1521   for (Int_t iLim=0; iLim<nLimits; ++iLim){
1522     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
1523   }
1524   
1525   delete arr;
1526   return binLimits;
1527 }
1528