7b8c04a885592f9affd8843aef159d7e2bc84125
[u/mrichter/AliRoot.git] / ANALYSIS / ANALYSISalice / 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 #include <AliTPCdEdxInfo.h>
41
42 #include <AliESDEvent.h>
43 #include <AliAODEvent.h>
44 #include <AliESDv0.h>
45 #include <AliAODv0.h>
46 #include <AliESDv0KineCuts.h>
47 #include <AliESDtrackCuts.h>
48
49 #include <AliMCEvent.h>
50
51 #include "AliAnalysisTaskPIDqa.h"
52
53
54 ClassImp(AliAnalysisTaskPIDqa)
55
56 //______________________________________________________________________________
57 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
58 AliAnalysisTaskSE(),
59 fPIDResponse(0x0),
60 fV0cuts(0x0),
61 fV0electrons(0x0),
62 fV0pions(0x0),
63 fV0kaons(0x0),
64 fV0protons(0x0),
65 fListQA(0x0),
66 fListQAits(0x0),
67 fListQAitsSA(0x0),
68 fListQAitsPureSA(0x0),
69 fListQAtpc(0x0),
70 fListQAtpcBasic(0x0),
71 fListQAtpcMCtruth(0x0),
72 fListQAtpcHybrid(0x0),
73 fListQAtpcOROChigh(0x0),
74 fListQAtpcV0(0x0),
75 fListQAtrd(0x0),
76 fListQAtrdNsig(0x0),
77 fListQAtrdNsigTPCTOF(0x0),
78 fListQAtof(0x0),
79 fListQAt0(0x0),
80 fListQAemcal(0x0),
81 fListQAhmpid(0x0),
82 fListQAtofhmpid(0x0),
83 fListQAtpctof(0x0),
84 fListQAV0(0x0),
85 fListQAinfo(0x0)
86 {
87   //
88   // Dummy constructor
89   //
90 }
91
92 //______________________________________________________________________________
93 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
94 AliAnalysisTaskSE(name),
95 fPIDResponse(0x0),
96 fV0cuts(0x0),
97 fV0electrons(0x0),
98 fV0pions(0x0),
99 fV0kaons(0x0),
100 fV0protons(0x0),
101 fListQA(0x0),
102 fListQAits(0x0),
103 fListQAitsSA(0x0),
104 fListQAitsPureSA(0x0),
105 fListQAtpc(0x0),
106 fListQAtpcBasic(0x0),
107 fListQAtpcMCtruth(0x0),
108 fListQAtpcHybrid(0x0),
109 fListQAtpcOROChigh(0x0),
110 fListQAtpcV0(0x0),
111 fListQAtrd(0x0),
112 fListQAtrdNsig(0x0),
113 fListQAtrdNsigTPCTOF(0x0),
114 fListQAtof(0x0),
115 fListQAt0(0x0),
116 fListQAemcal(0x0),
117 fListQAhmpid(0x0),
118 fListQAtofhmpid(0x0),
119 fListQAtpctof(0x0),
120 fListQAV0(0x0),
121 fListQAinfo(0x0)
122 {
123   //
124   // Default constructor
125   //
126   DefineInput(0,TChain::Class());
127   DefineOutput(1,TList::Class());
128 }
129
130 //______________________________________________________________________________
131 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
132 {
133   //
134   // Destructor
135   //
136
137   delete fV0cuts;
138   delete fV0electrons;
139   delete fV0pions;
140   delete fV0kaons;
141   delete fV0protons;
142
143   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA;
144 }
145
146 //______________________________________________________________________________
147 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
148 {
149   //
150   // Create the output QA objects
151   //
152
153   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
154
155   //input hander
156   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
157   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
158   if (!inputHandler) AliFatal("Input handler needed");
159
160   //pid response object
161   fPIDResponse=inputHandler->GetPIDResponse();
162   if (!fPIDResponse) AliError("PIDResponse object was not created");
163   
164   // V0 Kine cuts 
165   fV0cuts = new AliESDv0KineCuts;
166  
167   // V0 PID Obj arrays
168   fV0electrons = new TObjArray;
169   fV0pions     = new TObjArray;
170   fV0kaons     = new TObjArray;
171   fV0protons   = new TObjArray;
172
173   //
174   fListQA=new TList;
175   fListQA->SetOwner();
176   
177   fListQAits=new TList;
178   fListQAits->SetOwner();
179   fListQAits->SetName("ITS");
180
181   fListQAitsSA=new TList;
182   fListQAitsSA->SetOwner();
183   fListQAitsSA->SetName("ITS_SA");
184
185   fListQAitsPureSA=new TList;
186   fListQAitsPureSA->SetOwner();
187   fListQAitsPureSA->SetName("ITS_PureSA");
188
189   fListQAtpc=new TList;
190   fListQAtpc->SetOwner();
191   fListQAtpc->SetName("TPC");
192
193   fListQAtrd=new TList;
194   fListQAtrd->SetOwner();
195   fListQAtrd->SetName("TRD");
196
197   fListQAtrdNsig=new TList;
198   fListQAtrdNsig->SetOwner();
199   fListQAtrdNsig->SetName("TRDnSigma");
200   
201   fListQAtrdNsigTPCTOF=new TList;
202   fListQAtrdNsigTPCTOF->SetOwner();
203   fListQAtrdNsigTPCTOF->SetName("TRDnSigma_TPCTOF");
204   
205   fListQAtof=new TList;
206   fListQAtof->SetOwner();
207   fListQAtof->SetName("TOF");
208
209   fListQAt0=new TList;
210   fListQAt0->SetOwner();
211   fListQAt0->SetName("T0");
212   
213   fListQAemcal=new TList;
214   fListQAemcal->SetOwner();
215   fListQAemcal->SetName("EMCAL");
216   
217   fListQAhmpid=new TList;
218   fListQAhmpid->SetOwner();
219   fListQAhmpid->SetName("HMPID");
220   
221   fListQAtpctof=new TList;
222   fListQAtpctof->SetOwner();
223   fListQAtpctof->SetName("TPC_TOF");
224
225   fListQAtofhmpid=new TList;
226   fListQAtofhmpid->SetOwner();
227   fListQAtofhmpid->SetName("TOF_HMPID");
228   
229   fListQAV0=new TList;
230   fListQAV0->SetOwner();
231   fListQAV0->SetName("V0decay");
232
233   fListQAinfo=new TList;
234   fListQAinfo->SetOwner();
235   fListQAinfo->SetName("QAinfo");
236   
237   fListQA->Add(fListQAits);
238   fListQA->Add(fListQAitsSA);
239   fListQA->Add(fListQAitsPureSA);
240   fListQA->Add(fListQAtpc);
241   fListQA->Add(fListQAtrd);
242   fListQA->Add(fListQAtof);
243   fListQA->Add(fListQAt0);
244   fListQA->Add(fListQAemcal);
245   fListQA->Add(fListQAhmpid);
246   fListQA->Add(fListQAtpctof);
247   fListQA->Add(fListQAtofhmpid);
248   fListQA->Add(fListQAV0);
249   fListQA->Add(fListQAinfo);
250
251   SetupITSqa();
252 //  SetupTPCqa(kFALSE, kTRUE, kFALSE);
253   SetupTRDqa();
254   SetupTOFqa();
255   SetupT0qa();
256   SetupEMCALqa();
257   SetupHMPIDqa();
258   SetupTPCTOFqa();
259   SetupTOFHMPIDqa();
260   SetupV0qa();
261   SetupQAinfo();
262   
263   PostData(1,fListQA);
264 }
265
266
267 //______________________________________________________________________________
268 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
269 {
270   //
271   // Setup the PID response functions and fill the QA histograms
272   //
273
274   AliVEvent *event=InputEvent();
275   if (!event||!fPIDResponse) return;
276
277   // Start with the V0 task (only possible for ESDs?)
278   FillV0PIDlist();
279   
280   FillITSqa();
281   FillTPCqa();
282   FillTRDqa();
283   FillTOFqa();
284   FillEMCALqa();
285   FillHMPIDqa();
286   FillT0qa();
287   
288   //combined detector QA
289   FillTPCTOFqa();
290   FillTOFHMPIDqa();
291   
292   // Clear the V0 PID arrays
293   ClearV0PIDlist();
294
295   //QA info
296   FillQAinfo();
297   
298   PostData(1,fListQA);
299 }
300
301 //______________________________________________________________________________
302 void  AliAnalysisTaskPIDqa::FillV0PIDlist(){
303
304   //
305   // Fill the PID object arrays holding the pointers to identified particle tracks
306   //
307
308   // Dynamic cast to ESD events (DO NOTHING for AOD events)
309   AliESDEvent *event = dynamic_cast<AliESDEvent *>(InputEvent());
310   if ( !event )  return;
311   
312   if(TString(event->GetBeamType())=="Pb-Pb" || TString(event->GetBeamType())=="A-A"){
313     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb); 
314   }
315   else{
316     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP); 
317   }
318
319   // V0 selection
320   // set event
321   fV0cuts->SetEvent(event);
322
323   // loop over V0 particles
324   for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
325
326     AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0);
327  
328     if(!v0) continue;
329     if(v0->GetOnFlyStatus()) continue; 
330   
331     // Get the particle selection 
332     Bool_t foundV0 = kFALSE;
333     Int_t pdgV0, pdgP, pdgN;
334
335     foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
336     if(!foundV0) continue;
337     
338     Int_t iTrackP = v0->GetPindex();  // positive track
339     Int_t iTrackN = v0->GetNindex();  // negative track
340
341     // v0 Armenteros plot (QA)
342     Float_t armVar[2] = {0.0,0.0};
343     fV0cuts->Armenteros(v0, armVar);
344
345     TH2 *h=(TH2*)fListQAV0->At(0);
346     if (!h) continue;
347     h->Fill(armVar[0],armVar[1]);
348
349     // fill the Object arrays
350     // positive particles
351     if( pdgP == -11){
352       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
353     }
354     else if( pdgP == 211){
355       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
356     }
357     else if( pdgP == 321){
358       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackP));
359     }
360     else if( pdgP == 2212){
361       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
362     }
363
364     // negative particles
365     if( pdgN == 11){
366       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
367     }
368     else if( pdgN == -211){
369       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
370     }
371     else if( pdgN == -321){
372       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackN));
373     }
374     else if( pdgN == -2212){
375       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
376     }
377   
378
379   }
380 }
381 //______________________________________________________________________________
382 void  AliAnalysisTaskPIDqa::ClearV0PIDlist(){
383
384   //
385   // Clear the PID object arrays
386   //
387
388   fV0electrons->Clear();
389   fV0pions->Clear();
390   fV0kaons->Clear();
391   fV0protons->Clear();
392
393 }
394 //______________________________________________________________________________
395 void AliAnalysisTaskPIDqa::FillITSqa()
396 {
397   //
398   // Fill PID qa histograms for the ITS
399   //
400
401   AliVEvent *event=InputEvent();
402   
403   Int_t ntracks=event->GetNumberOfTracks();
404   for(Int_t itrack = 0; itrack < ntracks; itrack++){
405     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
406     ULong_t status=track->GetStatus();
407     // not that nice. status bits not in virtual interface
408     // ITS refit + ITS pid selection
409     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
410            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
411     Double_t mom=track->P();
412     
413     TList *theList = 0x0;
414     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
415       //ITS+TPC tracks
416       theList=fListQAits;
417     }else{
418       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){ 
419         //ITS Standalone tracks
420         theList=fListQAitsSA;
421       }else{
422         //ITS Pure Standalone tracks
423         theList=fListQAitsPureSA;
424       }
425     }
426     
427     
428     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
429       TH2 *h=(TH2*)theList->At(ispecie);
430       if (!h) continue;
431       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
432       h->Fill(mom,nSigma);
433     }
434     TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC);
435     if (h) {
436       Double_t sig=track->GetITSsignal();
437       h->Fill(mom,sig);
438     }
439   }
440 }
441
442
443 //______________________________________________________________________________
444 void AliAnalysisTaskPIDqa::FillTPCHistogramsSignal(TList *sublist, Int_t scenario, AliVTrack *track, Int_t mult)
445 {
446   //
447   // Fill PID qa histograms for the TPC: Fill the histograms for the TPC signal for different settings
448   //
449
450   AliMCEvent *eventMC=MCEvent();  // MC event for MC truth PID
451
452   Double_t mom=0.;      // track momentum
453   Double_t eta=0.;      // track eta
454   Double_t sig=0.;      // TPC dE/dx signal
455   Double_t sigStd=0.;         // TPC dE/dx signal (standard = all ROCs)
456   Double_t sigIROC=0.;        // TPC dE/dx signal (IROC) 
457   Double_t sigOROCmedium=0.;  // TPC dE/dx signal (OROCmedium) 
458   Double_t sigOROClong=0.;    // TPC dE/dx signal (OROClong) 
459   Double_t eleLineDist=0.;    // difference between TPC signal and electron expectation
460   Int_t trackLabel=0;   // label of the AliVTrack to identify the corresponding MCtrack
461   Int_t pdgCode=0;      // pdgcode of MC track for MC truth scenario
462   Int_t pdgCodeAbs=0;   // absolute value of pdgcode to get both particles and antiparticles
463   Int_t iSigMax=1;      // number of TPC signals (std = 1, set automatically higher if available)
464   Int_t nSpecies=0;     // number of particle species under study
465   Int_t count=0;        // counter for the number of plot sets for all species (i.e. nsigma vs. p, eta and mult)
466
467   mom=track->GetTPCmomentum();
468   eta=track->Eta();
469   sigStd=track->GetTPCsignal();
470
471   eleLineDist=sigStd-fPIDResponse->GetTPCResponse().GetExpectedSignal(track,AliPID::kElectron);
472
473   // Get number of particle species (less for V0 candidates = scenarios 40-44)
474   if (scenario > 39) nSpecies=(Int_t)AliPID::kSPECIES;
475   else nSpecies=(Int_t)AliPID::kSPECIESC;
476
477   // Set number of plot sets for all species
478   // (i.e. only nsigma vs. p => count=1; also vs. eta and mult => count=3)
479   if ( scenario == 1 || scenario > 39) count=3;
480   else count=1;
481
482   // Get MC track ( --> can be deleted if TPC signal is NOT filled for scenario=1 (MC truth)
483   if (eventMC) {
484     trackLabel=TMath::Abs(track->GetLabel());
485     AliVTrack *mcTrack=(AliVTrack*)eventMC->GetTrack(trackLabel);
486     pdgCode=mcTrack->PdgCode();
487     pdgCodeAbs=TMath::Abs(pdgCode);
488   }
489
490   // Get TPC dE/dx info and different TPC signals (IROC, OROCmedium, OROClong)
491   AliTPCdEdxInfo* fTPCdEdxInfo = 0x0;
492   fTPCdEdxInfo = track->GetTPCdEdxInfo();
493
494   if (fTPCdEdxInfo) {
495     sigIROC=fTPCdEdxInfo->GetTPCsignalShortPad();
496     sigOROCmedium=fTPCdEdxInfo->GetTPCsignalMediumPad();
497     sigOROClong=fTPCdEdxInfo->GetTPCsignalLongPad();
498     iSigMax=4;
499
500     //printf("mom = %.3f  sigStd = %.3f  sigIROC = %.3f  sigOROCmedium = %.3f  sigOROClong = %.3f \n",mom,sigStd,sigIROC,sigOROCmedium,sigOROClong);
501   }
502
503
504   // TPC signal for all particles vs. momentum (standard, IROC, OROCmedium, OROClong)
505   TH2 *h1std=(TH2*)sublist->At(count*nSpecies+4); 
506   if (h1std) {
507     h1std->Fill(mom,sigStd);
508   }
509
510   TH2 *h1iroc=(TH2*)sublist->At(count*nSpecies+5);
511   if ( h1iroc && sigIROC ) {
512     h1iroc->Fill(mom,sigIROC);
513   }
514
515   TH2 *h1orocm=(TH2*)sublist->At(count*nSpecies+6);
516   if  (h1orocm && sigOROCmedium ) {
517     h1orocm->Fill(mom,sigOROCmedium);
518   }
519
520   TH2 *h1orocl=(TH2*)sublist->At(count*nSpecies+7);
521   if ( h1orocl && sigOROClong ) {
522     h1orocl->Fill(mom,sigOROClong);
523   }
524
525
526   // - Beginn: MIP pions: TPC signal vs. eta, TPC signal vs. mult -
527   if (mom>0.45 && mom<0.5 && sigStd>40 && sigStd<60) {
528
529     Bool_t isPionMC=kTRUE;
530
531     if (scenario == 1) {
532       if ( pdgCodeAbs != 211 && pdgCodeAbs != 111 ) isPionMC=kFALSE;
533     }
534
535     // MIP pions: TPC signal vs. eta (standard, IROC, OROCmedium, OROClong)
536     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
537       if (iSig==0) sig=sigStd;
538       else if (iSig==1) sig=sigIROC;
539       else if (iSig==2) sig=sigOROCmedium;
540       else if (iSig==3) sig=sigOROClong;
541
542       TH2 *h2=(TH2*)sublist->At(count*nSpecies+8+iSig);
543       if ( h2 && isPionMC ) {
544         h2->Fill(eta,sig);
545       }
546     }
547
548     // MIP pions: TPC signal vs. mult (standard, IROC, OROCmedium, OROClong)
549     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
550       if (iSig==0) sig=sigStd;
551       else if (iSig==1) sig=sigIROC;
552       else if (iSig==2) sig=sigOROCmedium;
553       else if (iSig==3) sig=sigOROClong;
554
555       TH2 *h3=(TH2*)sublist->At(count*nSpecies+12+iSig);
556       if ( h3 && isPionMC && mult > 0 ) {
557         h3->Fill(mult,sig);
558       }
559     }
560   } // - End: MIP pions -
561
562   // - Beginn: Electrons: TPC signal vs. eta, TPC signal vs. mult -
563   if (mom>0.32 && mom<0.38 && eleLineDist>-10. && eleLineDist<15.) {
564
565     Bool_t isElectronMC=kTRUE;
566
567     if (scenario == 1) {
568       if ( pdgCodeAbs != 11 ) isElectronMC=kFALSE;
569     }
570
571     // Electrons: TPC signal vs. eta (standard, IROC, OROCmedium, OROClong)
572     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
573       if (iSig==0) sig=sigStd;
574       else if (iSig==1) sig=sigIROC;
575       else if (iSig==2) sig=sigOROCmedium;
576       else if (iSig==3) sig=sigOROClong;
577
578       TH2 *h4=(TH2*)sublist->At(count*nSpecies+16+iSig);
579       if ( h4 && isElectronMC ) {
580         h4->Fill(eta,sig);
581       }
582     }
583
584     // Electrons: TPC signal vs. mult (standard, IROC, OROCmedium, OROClong)
585     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
586       if (iSig==0) sig=sigStd;
587       else if (iSig==1) sig=sigIROC;
588       else if (iSig==2) sig=sigOROCmedium;
589       else if (iSig==3) sig=sigOROClong;
590
591       TH2 *h5=(TH2*)sublist->At(count*nSpecies+20+iSig);
592       if ( h5 && isElectronMC && mult > 0 ) {
593         h5->Fill(mult,sig);
594       }
595     }
596   } // - End: Electrons -
597
598 }
599
600 //______________________________________________________________________________
601 void AliAnalysisTaskPIDqa::FillTPCHistogramsNsigma(TList *sublist, Int_t scenario, AliVTrack *track, Int_t mult)
602 {
603   //
604   // Fill PID qa histograms for the TPC: Fill the histograms for TPC Nsigma for different settings
605   //
606
607   AliMCEvent *eventMC=MCEvent();  // MC event for MC truth PID
608
609   Double_t mom=0.;      // track momentum
610   Double_t eta=0.;      // track eta
611   Double_t nSigma=0.;   // number of sigmas wrt. expected signal
612   Double_t sig=0.;      // TPC dE/dx signal
613   Double_t eleLineDist=0.;  // difference between TPC signal and electron expectation
614   Int_t trackLabel=0;   // label of the AliVTrack to identify the corresponding MCtrack
615   Int_t pdgCode=0;      // pdgcode of MC track for MC truth scenario
616   Int_t pdgCodeAbs=0;   // absolute value of pdgcode to get both particles and antiparticles
617   Int_t nSpecies=0;     // number of particle species under study
618   Int_t count=0;        // counter for the number of plot sets for all species (i.e. vs. p, eta and mult)
619
620   mom=track->GetTPCmomentum();
621   eta=track->Eta();
622   sig=track->GetTPCsignal();
623
624   eleLineDist=sig-fPIDResponse->GetTPCResponse().GetExpectedSignal(track,AliPID::kElectron);
625
626   // Get number of particle species (less for V0 candidates = scenarios 40-44)
627   if (scenario > 39) nSpecies=(Int_t)AliPID::kSPECIES;
628   else nSpecies=(Int_t)AliPID::kSPECIESC;
629
630   // Set number of plot sets for all species
631   // (i.e. only vs. p => count=1; also vs. eta and mult => count=3)
632   if ( scenario == 1 || scenario > 39 ) count=3;
633   else count=1;
634
635   // Get MC track
636   if (eventMC) {
637     trackLabel=TMath::Abs(track->GetLabel());
638     AliVTrack *mcTrack=(AliVTrack*)eventMC->GetTrack(trackLabel);
639     pdgCode=mcTrack->PdgCode();
640     pdgCodeAbs=TMath::Abs(pdgCode);
641   }
642
643
644   // - Beginn: Nsigma vs. p, vs. eta and vs. multiplicity for different particle species -
645   for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
646
647     TH2 *h=(TH2*)sublist->At(ispecie);
648     if (!h) continue;
649
650     if (scenario == 1) {
651       if ( ispecie == 0 && pdgCodeAbs != 11 ) continue;  // Electron
652       if ( ispecie == 1 && pdgCodeAbs != 13 ) continue;  // Muon
653       if ( ispecie == 2 && pdgCodeAbs != 211 && pdgCodeAbs!=111 ) continue;  // Pion
654       if ( ispecie == 3 && pdgCodeAbs != 321 && pdgCodeAbs!=311 ) continue;  // Kaon
655       if ( ispecie == 4 && pdgCodeAbs != 2212 ) continue;  // Proton
656       if ( ispecie == 5 && pdgCodeAbs != 1000010020 ) continue;  // Deuteron
657       if ( ispecie == 6 && pdgCodeAbs != 1000010030 ) continue;  // Triton
658       if ( ispecie == 7 && pdgCodeAbs != 1000020030 ) continue;  // Helium-3
659       if ( ispecie == 8 && pdgCodeAbs != 1000020040 ) continue;  // Alpha
660     }
661     else if (scenario > 39) {
662       if ( ispecie == 0 && scenario != 40 ) continue;  // Electron
663       if ( ispecie == 1 ) continue;  // Muon
664       if ( ispecie == 2 && scenario != 42 ) continue;  // Pion
665       if ( ispecie == 3 && scenario != 43 ) continue;  // Kaon
666       if ( ispecie == 4 && scenario != 44 ) continue;  // Proton
667     }
668
669     if (scenario == 2) {
670       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxHybrid);
671     }
672     else if (scenario == 3) {
673       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxOROC);
674     }
675     else {
676       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
677     }
678
679     h->Fill(mom,nSigma);
680
681     if (count == 3) {
682       TH2 *hEta=(TH2*)sublist->At(ispecie+nSpecies);
683       TH2 *hMult=(TH2*)sublist->At(ispecie+2*nSpecies);
684  
685       if ( hEta ) hEta->Fill(eta,nSigma);
686       if ( hMult && mult > 0 ) hMult->Fill(mult,nSigma);
687     }
688   } // - End: different particle species -
689
690
691   // -- Beginn: Fill histograms for MIP pions and electrons (only for some scenarios) --
692   if ( scenario == 0 || scenario == 2 || scenario == 3 ) {
693
694     // - Beginn: MIP pions: Nsigma vs. eta, Nsigma vs. mult -
695     if (mom>0.45 && mom<0.5 && sig>40 && sig<60) {
696
697       Bool_t isPionMC=kTRUE;
698
699       TH2 *h1=(TH2*)sublist->At(count*nSpecies);
700       if (h1) {
701         if (scenario == 1) {
702           if ( pdgCodeAbs != 211 && pdgCodeAbs != 111 ) isPionMC=kFALSE;
703           if (isPionMC) {
704             nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
705           }
706         }
707         else if (scenario == 2) {
708           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion, AliTPCPIDResponse::kdEdxHybrid);
709         }
710         else if (scenario == 3) {
711           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion, AliTPCPIDResponse::kdEdxOROC);
712         }
713         else nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
714
715         if (isPionMC) h1->Fill(eta,nSigma);
716       }
717
718       TH2 *h2m=(TH2*)sublist->At(count*nSpecies+1);
719       if ( h2m && isPionMC && mult > 0 ) {
720         h2m->Fill(mult,nSigma);
721       }
722    
723     } // - End: MIP pions -
724
725     // - Beginn: Electrons: Nsigma vs. eta, Nsigma vs. mult -
726     if (mom>0.32 && mom<0.38 && eleLineDist>-10. && eleLineDist<15.) {
727
728       Bool_t isElectronMC=kTRUE;
729
730       TH2 *h3=(TH2*)sublist->At(count*nSpecies+2);
731       if (h3) {
732         if (scenario == 1) {
733           if ( pdgCodeAbs != 11 ) isElectronMC=kFALSE;
734           if (isElectronMC) {
735             nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
736           }
737         }
738         if (scenario == 2) {
739           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxHybrid);
740         }
741         else if (scenario == 3) {
742           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxOROC);
743         }
744         else nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
745
746         if (isElectronMC) h3->Fill(eta,nSigma);
747       }
748
749       TH2 *h4m=(TH2*)sublist->At(count*nSpecies+3);
750       if ( h4m && isElectronMC && mult > 0 ) {
751         h4m->Fill(mult,nSigma);
752       }
753
754     } // - End: Electrons -
755   } // -- End: Fill histograms for MIP pions and electrons --
756
757 }
758
759 //______________________________________________________________________________
760 void AliAnalysisTaskPIDqa::FillTPCqa()
761 {
762   //
763   // Fill PID qa histograms for the TPC
764   //
765
766   // switches for the different scenarios
767   Bool_t scBasic=1;     // default/basic
768   Bool_t scMCtruth=1;   // for MC truth tracks
769   Bool_t scHybrid=1;    // for hybrid PID (only LHC11h)
770   Bool_t scOROChigh=1;  // only OROC signal (only LHC11h)
771   Bool_t scV0=1;        // for V0 candidates (only for ESDs available)
772   Int_t scCounter=0;    // counter of scenarios, used for the histograms at the end of FillTPCqa
773
774   // input handler
775   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
776   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
777   if (!inputHandler) AliFatal("Input handler needed");
778
779   AliVEvent *event=InputEvent();
780
781   // ESD or AOD event needed to get reference multiplicity (not in AliVEvent)
782   AliAODEvent *fAODevent = 0x0;   // AOD event
783   AliESDEvent *fESDevent = 0x0;   // ESD event
784   AliESDtrackCuts *esdTrackCuts = 0x0;  // ESD track Cuts (ref mult is in AliESDtrackCuts)
785
786   Double_t eta=0.;    // track eta
787   Int_t mult=0;       // event multiplicity (TPConlyRefMult)
788   //Int_t nacc=0;       // counter for accepted multiplicity
789
790   // Check for MC
791   scMCtruth=(MCEvent()!=0x0);
792
793   // Check if period is data LHC11h by checking if
794   // the splines for ALLhigh have been set by AliPIDResponse
795   AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
796   if (tpcResp.GetResponseFunction(AliPID::kPion, AliTPCPIDResponse::kALLhigh)==0x0) {
797     scHybrid   = kFALSE;
798     scOROChigh = kFALSE;
799   }
800
801   // Check if "ESD" or "AOD" and get the corresponding event and the beam type (or centrality)
802   TString analysisType = inputHandler->GetDataType(); // can be "ESD" or "AOD"
803   if (analysisType == "ESD") {
804     fESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
805     esdTrackCuts = new AliESDtrackCuts("esdTrackCuts");
806     //printf("\n--- New event - event type = ESD \n");
807   }
808   else if (analysisType == "AOD") {
809     fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
810     //printf("\n--- New event - event type = AOD \n");
811
812     // disable V0 scenario, because V0s are not available for AODs in this current implementation
813     scV0=0;
814   }
815
816   // Check if Basic list is already created
817   // If not: Go to SetupTPCqa and creat lists and histograms
818   if(!fListQAtpcBasic) {
819     //printf("\n--- No list QA TPC Basic found -> go to SetupTPCqa! ---\n");
820     SetupTPCqa(scMCtruth, scHybrid, scV0);
821   }
822
823   // Get the number of scenarios by counting those, which are switched on
824   if (scBasic) scCounter++;
825   if (scMCtruth) scCounter++;
826   if (scHybrid) scCounter++;
827   if (scOROChigh) scCounter++;
828   if (scV0) scCounter++;
829
830   // Get reference multiplicity for ESDs
831   if ( analysisType == "ESD" && esdTrackCuts ) {
832     mult=esdTrackCuts->GetReferenceMultiplicity(fESDevent,kTRUE);
833   }
834
835   // Get reference multiplicity for AODs
836   if ( analysisType == "AOD" && fAODevent ) {
837     AliAODHeader * header=dynamic_cast<AliAODHeader*>(fAODevent->GetHeader());
838     if(!header) AliFatal("Not a standard AOD");
839     mult=header->GetTPConlyRefMultiplicity();
840   }
841
842   /*if (mult < 0) {
843     printf("Reference multiplicity not available \n");
844     //return;
845   }*/
846
847   //printf("The multiplicity is = %i ",mult);
848
849
850   // -- Begin: track loop --
851   Int_t ntracks=event->GetNumberOfTracks();
852   for(Int_t itrack = 0; itrack < ntracks; itrack++){
853     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
854
855     //
856     //basic track cuts
857     //
858     ULong_t status=track->GetStatus();
859     // not that nice. status bits not in virtual interface
860     // TPC refit + ITS refit + TPC pid
861     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
862         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
863
864     // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
865     //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
866     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
867     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
868     if (track->GetTPCNclsF()>0) {
869       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
870     }
871
872     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
873
874     eta=track->Eta();
875     if ( TMath::Abs(eta)>0.9 ) continue;
876
877     //nacc++; // counter for accepted multiplicity
878
879     // the default ("basic") scenario
880     if (scBasic == 1) {
881       FillTPCHistogramsNsigma(fListQAtpcBasic,0,track,mult);
882       FillTPCHistogramsSignal(fListQAtpcBasic,0,track,mult);
883     }
884
885     // only MC truth identified particles
886     if (scMCtruth == 1) {
887       FillTPCHistogramsNsigma(fListQAtpcMCtruth,1,track,mult);
888     }
889
890     // the "hybrid" scenario (only for LHC11h)
891     if (scHybrid == 1) {
892       FillTPCHistogramsNsigma(fListQAtpcHybrid,2,track,mult);
893     }
894
895     // the "OROC high" scenario (only for LHC11h)
896     if (scOROChigh == 1) {
897       FillTPCHistogramsNsigma(fListQAtpcOROChigh,3,track,mult);
898     }
899
900   } // -- End: track loop --
901
902
903   // -- Begin: track loops for V0 candidates --
904   if (scV0 == 1) {
905
906     // - Begin: track loop for electrons from V0 -
907     for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
908       AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
909
910       //
911       //basic track cuts
912       //
913       ULong_t status=track->GetStatus();
914       // not that nice. status bits not in virtual interface
915       // TPC refit + ITS refit + TPC pid
916       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
917           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
918
919       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
920       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
921       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
922       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
923       if (track->GetTPCNclsF()>0) {
924         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
925       }
926
927       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
928
929       eta=track->Eta();
930       if ( TMath::Abs(eta)>0.9 ) continue;
931
932       // fill histograms for V0 candidates
933       FillTPCHistogramsNsigma(fListQAtpcV0,40,track,mult);
934
935     } // - End: track loop for electrons from V0 -
936
937
938     // - Begin: track loop for pions from V0 -
939     for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
940       AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
941
942       //
943       //basic track cuts
944       //
945       ULong_t status=track->GetStatus();
946       // not that nice. status bits not in virtual interface
947       // TPC refit + ITS refit + TPC pid
948       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
949           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
950
951       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
952       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
953       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
954       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
955       if (track->GetTPCNclsF()>0) {
956         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
957       }
958
959       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
960
961       eta=track->Eta();
962       if ( TMath::Abs(eta)>0.9 ) continue;
963
964       // fill histograms for V0 candidates
965       FillTPCHistogramsNsigma(fListQAtpcV0,42,track,mult);
966
967     } // - End: track loop for pions from V0 -
968
969
970     // - Begin: track loop for kaons from V0 -
971     for(Int_t itrack = 0; itrack < fV0kaons->GetEntries(); itrack++){
972       AliVTrack *track=(AliVTrack*)fV0kaons->At(itrack);
973
974       //
975       //basic track cuts
976       //
977       ULong_t status=track->GetStatus();
978       // not that nice. status bits not in virtual interface
979       // TPC refit + ITS refit + TPC pid
980       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
981           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
982
983       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
984       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
985       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
986       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
987       if (track->GetTPCNclsF()>0) {
988         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
989       }
990
991       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
992
993       eta=track->Eta();
994       if ( TMath::Abs(eta)>0.9 ) continue;
995
996       // fill histograms for V0 candidates
997       FillTPCHistogramsNsigma(fListQAtpcV0,43,track,mult);
998
999     } // - End: track loop for kaons from V0 -
1000
1001
1002     // - Begin: track loop for protons from V0 -
1003     for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1004       AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1005
1006       //
1007       //basic track cuts
1008       //
1009       ULong_t status=track->GetStatus();
1010       // not that nice. status bits not in virtual interface
1011       // TPC refit + ITS refit + TPC pid
1012       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1013           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1014
1015       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
1016       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
1017       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1018       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1019       if (track->GetTPCNclsF()>0) {
1020         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1021       }
1022
1023       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1024
1025       eta=track->Eta();
1026       if ( TMath::Abs(eta)>0.9 ) continue;
1027
1028       // fill histograms for V0 candidates
1029       FillTPCHistogramsNsigma(fListQAtpcV0,44,track,mult);
1030
1031     } // - End: track loop for protons from V0 -
1032
1033   } // -- End: track loops for V0 candidates --
1034
1035
1036   // Multiplicity distribution
1037   TH1 *hm=(TH1*)fListQAtpc->At(scCounter);
1038   if (hm) {
1039     hm->Fill(mult);
1040   }
1041
1042   //printf("\nAccepted multiplicity = %i \n --- END of event --- \n",nacc);
1043
1044 }
1045
1046 //______________________________________________________________________________
1047 void AliAnalysisTaskPIDqa::FillTRDqa()
1048 {
1049   //
1050   // Fill PID qa histograms for the TRD
1051   //
1052   AliVEvent *event=InputEvent();
1053   Int_t ntracks = event->GetNumberOfTracks();
1054   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
1055     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
1056
1057     //
1058     //basic track cuts
1059     //
1060     ULong_t status=track->GetStatus();
1061     // not that nice. status bits not in virtual interface
1062     // TPC refit + ITS refit + TPC pid + TRD out
1063     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1064         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1065 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) || //removes light nuclei. So it is out for the moment
1066         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
1067     
1068     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1069     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1070     if (track->GetTPCNclsF()>0) {
1071       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1072     }
1073     
1074     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1075
1076     Double_t likelihoods[AliPID::kSPECIES];
1077     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
1078     Int_t ntracklets = 0;
1079     Double_t momentum = -1.;
1080     for(Int_t itl = 0; itl < 6; itl++) {
1081       if(track->GetTRDmomentum(itl) > 0.) {
1082         ntracklets++;
1083         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
1084       }
1085     }
1086     
1087     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1088       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
1089       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
1090     }
1091
1092     //=== nSigma and signal ===
1093     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1094       TH2 *h=(TH2*)fListQAtrdNsig->At(ispecie);
1095       TH2 *hTPCTOF=(TH2*)fListQAtrdNsigTPCTOF->At(ispecie);
1096       if (!h || !hTPCTOF) continue;
1097       Float_t nSigmaTPC=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, (AliPID::EParticleType)ispecie);
1098       Float_t nSigmaTRD=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, (AliPID::EParticleType)ispecie);
1099       Float_t nSigmaTOF=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, (AliPID::EParticleType)ispecie);
1100       h->Fill(momentum,nSigmaTRD);
1101
1102       if (TMath::Abs(nSigmaTPC)<3 && TMath::Abs(nSigmaTOF)<3) {
1103         hTPCTOF->Fill(momentum,nSigmaTRD);
1104       }
1105     }
1106
1107     TH2 *h=(TH2*)fListQAtrdNsig->Last();
1108     
1109     if (h) {
1110       Double_t sig=track->GetTRDsignal();
1111       h->Fill(momentum,sig);
1112     }
1113     
1114   }
1115 }
1116
1117 //______________________________________________________________________________
1118 void AliAnalysisTaskPIDqa::FillTOFqa()
1119 {
1120   //
1121   // Fill TOF information
1122   //
1123   AliVEvent *event=InputEvent();
1124
1125   Int_t ntracks=event->GetNumberOfTracks();
1126   Int_t tracksAtTof = 0;
1127   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1128     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1129
1130     //
1131     //basic track cuts
1132     //
1133     ULong_t status=track->GetStatus();
1134     // TPC refit + ITS refit +
1135     // TOF out + kTIME
1136     // kTIME
1137     // (we don't use kTOFmismatch because it depends on TPC and kTOFpid because it prevents light nuclei
1138     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1139         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1140         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1141         //        !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1142         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1143
1144     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1145     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1146     if (track->GetTPCNclsF()>0) {
1147       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1148     }
1149
1150     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1151
1152     tracksAtTof++;
1153
1154     Double_t mom=track->P();
1155
1156     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1157       TH2 *h=(TH2*)fListQAtof->At(ispecie);
1158       if (!h) continue;
1159       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1160       h->Fill(mom,nSigma);
1161     }
1162
1163     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
1164     if (h) {
1165       Double_t sig=track->GetTOFsignal()/1000.;
1166       h->Fill(mom,sig);
1167     }
1168
1169     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
1170     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
1171
1172     if (mom >= 0.75 && mom <= 1.25 ) {
1173       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
1174       if (mask == 0) {
1175         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
1176       } else if (mask == 1) {
1177         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
1178       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
1179         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
1180       } else {
1181         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
1182       }
1183       if (mask & 0x1) { //at least TOF-T0 present
1184         Double_t delta=0;
1185         (void)fPIDResponse->GetSignalDelta((AliPIDResponse::EDetector)AliPIDResponse::kTOF,track,(AliPID::EParticleType)AliPID::kPion,delta);
1186         ((TH1F*)fListQAtof->FindObject("hDelta_TOF_Pion"))->Fill(delta);
1187       }
1188     }
1189
1190     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
1191     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
1192
1193     Double_t startTimeT0 = event->GetT0TOF(0);
1194     if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1195     else {
1196       startTimeT0 = event->GetT0TOF(1);
1197       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1198       startTimeT0 = event->GetT0TOF(2);
1199       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1200     }
1201   }
1202   if (tracksAtTof > 0) {
1203     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
1204     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
1205     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
1206   }
1207 }
1208
1209 //______________________________________________________________________________
1210 void AliAnalysisTaskPIDqa::FillT0qa()
1211 {
1212   //
1213   // Fill TOF information
1214   //
1215   AliVEvent *event=InputEvent();
1216
1217   Int_t ntracks=event->GetNumberOfTracks();
1218
1219   Int_t tracksAtT0 = 0;
1220
1221   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1222     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1223
1224     //
1225     //basic track cuts
1226     //
1227     ULong_t status=track->GetStatus();
1228     // TPC refit + ITS refit +
1229     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1230         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1231     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1232     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1233     if (track->GetTPCNclsF()>0) {
1234       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1235     }
1236     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1237
1238     tracksAtT0++;
1239   }
1240
1241   Bool_t t0A = kFALSE;
1242   Bool_t t0C = kFALSE;
1243   Bool_t t0And = kFALSE;
1244   Double_t startTimeT0 = event->GetT0TOF(0);     // AND
1245   if (startTimeT0 < 90000) {
1246     t0And = kTRUE;
1247     ((TH1F*)fListQAt0->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1248     }
1249   startTimeT0 = event->GetT0TOF(1);             // T0A 
1250   if (startTimeT0 < 90000) {
1251     t0A = kTRUE;
1252     ((TH1F*)fListQAt0->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1253     
1254   }
1255   startTimeT0 = event->GetT0TOF(2);             // T0C 
1256   if (startTimeT0 < 90000) {
1257     t0C = kTRUE;
1258     ((TH1F*)fListQAt0->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1259   }
1260   
1261   ((TH1F* )fListQAt0->FindObject("hnTracksAt_T0"))->Fill(tracksAtT0);
1262   if (t0A) ((TH1F*)fListQAt0->FindObject("hT0AEff"))->Fill(tracksAtT0);
1263   if (t0C) ((TH1F*)fListQAt0->FindObject("hT0CEff"))->Fill(tracksAtT0);
1264   if (t0And) ((TH1F*)fListQAt0->FindObject("hT0AndEff"))->Fill(tracksAtT0);
1265   if (t0A || t0C) ((TH1F*)fListQAt0->FindObject("hT0OrEff"))->Fill(tracksAtT0);
1266 }
1267
1268
1269 //______________________________________________________________________________
1270 void AliAnalysisTaskPIDqa::FillEMCALqa()
1271 {
1272   //
1273   // Fill PID qa histograms for the EMCAL
1274   //
1275
1276   AliVEvent *event=InputEvent();
1277   
1278   Int_t ntracks=event->GetNumberOfTracks();
1279   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1280     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1281     
1282     //
1283     //basic track cuts
1284     //
1285     ULong_t status=track->GetStatus();
1286     // not that nice. status bits not in virtual interface
1287     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1288
1289     Double_t pt=track->Pt();
1290    
1291     //EMCAL nSigma (only for electrons at the moment)
1292     TH2 *h=(TH2*)fListQAemcal->At(0);
1293     if (!h) continue;
1294     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
1295     h->Fill(pt,nSigma);
1296     
1297   }
1298
1299    //EMCAL signal (E/p vs. pT) for electrons from V0
1300   for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
1301     AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
1302
1303     //
1304     //basic track cuts
1305     //
1306     ULong_t status=track->GetStatus();
1307     // not that nice. status bits not in virtual interface
1308     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1309
1310     Double_t pt=track->Pt();
1311
1312     TH2 *h=(TH2*)fListQAemcal->At(1);
1313     if (h) {
1314
1315       Int_t nMatchClus = track->GetEMCALcluster();
1316       Double_t mom     = track->P();
1317       Double_t eop     = -1.;
1318
1319       if(nMatchClus > -1){
1320     
1321         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1322
1323         if(matchedClus){
1324
1325           // matched cluster is EMCAL
1326           if(matchedClus->IsEMCAL()){
1327
1328             Double_t fClsE       = matchedClus->E();
1329             eop                  = fClsE/mom;
1330
1331             h->Fill(pt,eop);
1332
1333           }
1334         }
1335       }
1336     }
1337   }
1338
1339    //EMCAL signal (E/p vs. pT) for pions from V0
1340   for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
1341     AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
1342
1343     //
1344     //basic track cuts
1345     //
1346     ULong_t status=track->GetStatus();
1347     // not that nice. status bits not in virtual interface
1348     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1349
1350     Double_t pt=track->Pt();
1351
1352     TH2 *h=(TH2*)fListQAemcal->At(2);
1353     if (h) {
1354
1355       Int_t nMatchClus = track->GetEMCALcluster();
1356       Double_t mom     = track->P();
1357       Double_t eop     = -1.;
1358
1359       if(nMatchClus > -1){
1360     
1361         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1362
1363         if(matchedClus){
1364
1365           // matched cluster is EMCAL
1366           if(matchedClus->IsEMCAL()){
1367
1368             Double_t fClsE       = matchedClus->E();
1369             eop                  = fClsE/mom;
1370
1371             h->Fill(pt,eop);
1372
1373           }
1374         }
1375       }
1376     }
1377   }
1378
1379    //EMCAL signal (E/p vs. pT) for protons from V0
1380   for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1381     AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1382
1383     //
1384     //basic track cuts
1385     //
1386     ULong_t status=track->GetStatus();
1387     // not that nice. status bits not in virtual interface
1388     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1389
1390     Double_t pt=track->Pt();
1391
1392     TH2 *hP=(TH2*)fListQAemcal->At(3);
1393     TH2 *hAP=(TH2*)fListQAemcal->At(4);
1394     if (hP && hAP) {
1395
1396       Int_t nMatchClus = track->GetEMCALcluster();
1397       Double_t mom     = track->P();
1398       Int_t charge     = track->Charge();             
1399       Double_t eop     = -1.;
1400
1401       if(nMatchClus > -1){
1402     
1403         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1404
1405         if(matchedClus){
1406
1407           // matched cluster is EMCAL
1408           if(matchedClus->IsEMCAL()){
1409
1410             Double_t fClsE       = matchedClus->E();
1411             eop                  = fClsE/mom;
1412
1413             if(charge > 0)      hP->Fill(pt,eop);
1414             else if(charge < 0) hAP->Fill(pt,eop);
1415
1416           }
1417         }
1418       }
1419     }
1420   }
1421
1422 }
1423
1424
1425 //______________________________________________________________________________
1426 void AliAnalysisTaskPIDqa::FillHMPIDqa()
1427 {
1428   //
1429   // Fill PID qa histograms for the HMPID
1430   //
1431   
1432   AliVEvent *event=InputEvent();
1433   
1434   Int_t ntracks=event->GetNumberOfTracks();
1435   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1436     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1437     
1438     //
1439     //basic track cuts
1440     //
1441     const ULong_t status=track->GetStatus();
1442     // not that nice. status bits not in virtual interface
1443     // TPC refit + ITS refit +
1444     // TOF out + TOFpid +
1445     // kTIME
1446     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1447         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1448
1449     const Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1450     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1451     if (track->GetTPCNclsF()>0) {
1452       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1453     }
1454
1455     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1456     
1457     const Double_t mom = track->P();
1458     const Double_t ckovAngle = track->GetHMPIDsignal();
1459
1460     Int_t nhists=0;
1461     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1462       if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
1463       TH2 *h=(TH2*)fListQAhmpid->At(nhists);
1464       if (!h) {++nhists; continue;}
1465       const Double_t nSigma=fPIDResponse->NumberOfSigmasHMPID(track, (AliPID::EParticleType)ispecie);
1466       h->Fill(mom,nSigma);
1467       ++nhists;
1468     }
1469     
1470     TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(AliPID::kSPECIESC);
1471     
1472     if (hThetavsMom) hThetavsMom->Fill(mom,ckovAngle);
1473   
1474   }
1475 }
1476 //______________________________________________________________________________
1477 void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
1478 {
1479   //
1480   // Fill PID qa histograms for the HMPID
1481   //
1482   
1483   AliVEvent *event=InputEvent();
1484   
1485   Int_t ntracks=event->GetNumberOfTracks();
1486   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1487     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1488     
1489     //
1490     //basic track cuts
1491     //
1492     ULong_t status=track->GetStatus();
1493     // not that nice. status bits not in virtual interface
1494     // TPC refit + ITS refit +
1495     // TOF out + TOFpid +
1496     // kTIME
1497     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1498         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1499         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1500         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1501         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1502
1503     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1504     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1505     if (track->GetTPCNclsF()>0) {
1506       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1507     }
1508
1509     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1510     
1511     Double_t mom = track->P();
1512     Double_t ckovAngle = track->GetHMPIDsignal();
1513     
1514     Double_t nSigmaTOF[3]; 
1515     TH1F *h[3];
1516     
1517     for (Int_t ispecie=2; ispecie<5; ++ispecie){
1518       //TOF nSigma
1519       nSigmaTOF[ispecie-2]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1520       h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
1521       
1522     if(TMath::Abs(nSigmaTOF[0])<2)                                                              h[0]->Fill(mom,ckovAngle);
1523     
1524     if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3)                                h[1]->Fill(mom,ckovAngle);
1525
1526     if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3)  h[2]->Fill(mom,ckovAngle);
1527       
1528   }
1529   
1530 }
1531
1532 //______________________________________________________________________________
1533 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
1534 {
1535   //
1536   // Fill PID qa histograms for the TOF
1537   //   Here also the TPC histograms after TOF selection are filled
1538   //
1539
1540   AliVEvent *event=InputEvent();
1541
1542   Int_t ntracks=event->GetNumberOfTracks();
1543   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1544     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1545
1546     //
1547     //basic track cuts
1548     //
1549     ULong_t status=track->GetStatus();
1550     // not that nice. status bits not in virtual interface
1551     // TPC refit + ITS refit +
1552     // TOF out + TOFpid +
1553     // kTIME
1554     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1555         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1556 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
1557         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1558         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1559         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1560
1561     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1562     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1563     if (track->GetTPCNclsF()>0) {
1564       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1565     }
1566
1567     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1568
1569
1570     Double_t mom=track->P();
1571     Double_t momTPC=track->GetTPCmomentum();
1572
1573     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1574       //TOF nSigma
1575       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1576       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
1577
1578       //TPC after TOF cut
1579       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
1580       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
1581
1582       //TOF after TPC cut
1583       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
1584       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
1585
1586       //EMCAL after TOF and TPC cut
1587       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
1588       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
1589
1590         Int_t nMatchClus = track->GetEMCALcluster();
1591         Double_t pt      = track->Pt();
1592         Double_t eop     = -1.;
1593         
1594         if(nMatchClus > -1){
1595           
1596           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1597           
1598           if(matchedClus){
1599             
1600             // matched cluster is EMCAL
1601             if(matchedClus->IsEMCAL()){
1602               
1603               Double_t fClsE       = matchedClus->E();
1604               eop                  = fClsE/mom;
1605
1606               h->Fill(pt,eop);
1607  
1608               
1609             }
1610           }
1611         }
1612       }
1613     }
1614   }
1615 }
1616
1617 //_____________________________________________________________________________
1618 void AliAnalysisTaskPIDqa::FillQAinfo()
1619 {
1620   //
1621   // Fill the QA information
1622   //
1623
1624
1625   //TPC QA info
1626   TObjArray *arrTPC=static_cast<TObjArray*>(fListQAinfo->At(0));
1627   if (fPIDResponse && arrTPC){
1628     AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
1629     // fill spline names
1630     if (!arrTPC->UncheckedAt(0)){
1631       
1632       TObjArray *arrTPCsplineNames=new TObjArray(AliPID::kSPECIESC);
1633       arrTPCsplineNames->SetOwner();
1634       arrTPCsplineNames->SetName("TPC_spline_names");
1635       arrTPC->AddAt(arrTPCsplineNames,0);
1636       
1637       for (Int_t iresp=0; iresp<AliPID::kSPECIESC; ++iresp){
1638         const TObject *o=tpcResp.GetResponseFunction((AliPID::EParticleType)iresp);
1639         if (!o) continue;
1640         arrTPCsplineNames->Add(new TObjString(Form("%02d: %s",iresp, o->GetName())));
1641       }
1642     }
1643
1644     // tpc response config
1645     if (!arrTPC->UncheckedAt(1)){
1646       
1647       TObjArray *arrTPCconfigInfo=new TObjArray;
1648       arrTPCconfigInfo->SetOwner();
1649       arrTPCconfigInfo->SetName("TPC_config_info");
1650       arrTPC->AddAt(arrTPCconfigInfo,1);
1651
1652       TObjString *ostr=0x0;
1653       ostr=new TObjString;
1654       ostr->String().Form("Eta Corr map: %s", tpcResp.GetEtaCorrMap()?tpcResp.GetEtaCorrMap()->GetName():"none");
1655       arrTPCconfigInfo->Add(ostr);
1656
1657       ostr=new TObjString;
1658       ostr->String().Form("Sigma Par map: %s", tpcResp.GetSigmaPar1Map()?tpcResp.GetSigmaPar1Map()->GetName():"none");
1659       arrTPCconfigInfo->Add(ostr);
1660
1661       ostr=new TObjString;
1662       ostr->String().Form("MIP: %.2f", tpcResp.GetMIP());
1663       arrTPCconfigInfo->Add(ostr);
1664       
1665       ostr=new TObjString;
1666       ostr->String().Form("Res: Def %.3g (%.3g) : AllHigh %.3g (%.3g) : OROC high %.3g (%.3g)",
1667                           tpcResp.GetRes0(AliTPCPIDResponse::kDefault), tpcResp.GetResN2(AliTPCPIDResponse::kDefault),
1668                           tpcResp.GetRes0(AliTPCPIDResponse::kALLhigh), tpcResp.GetResN2(AliTPCPIDResponse::kALLhigh),
1669                           tpcResp.GetRes0(AliTPCPIDResponse::kOROChigh), tpcResp.GetResN2(AliTPCPIDResponse::kOROChigh)
1670                          );
1671       arrTPCconfigInfo->Add(ostr);
1672     }
1673   }
1674 }
1675
1676 //______________________________________________________________________________
1677 void AliAnalysisTaskPIDqa::SetupITSqa()
1678 {
1679   //
1680   // Create the ITS qa objects
1681   //
1682   
1683   TVectorD *vX=MakeLogBinning(200,.1,30);
1684   
1685   //ITS+TPC tracks
1686   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1687     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
1688                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1689                               vX->GetNrows()-1,vX->GetMatrixArray(),
1690                               200,-10,10);
1691     fListQAits->Add(hNsigmaP);
1692   }
1693   TH2F *hSig = new TH2F("hSigP_ITS",
1694                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1695                         vX->GetNrows()-1,vX->GetMatrixArray(),
1696                         300,0,300);
1697   fListQAits->Add(hSig);
1698
1699   //ITS Standalone tracks
1700   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1701     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
1702                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1703                                 vX->GetNrows()-1,vX->GetMatrixArray(),
1704                                 200,-10,10);
1705     fListQAitsSA->Add(hNsigmaPSA);
1706   }
1707   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
1708                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1709                           vX->GetNrows()-1,vX->GetMatrixArray(),
1710                           300,0,300);
1711   fListQAitsSA->Add(hSigSA);
1712   
1713   //ITS Pure Standalone tracks
1714   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1715     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
1716                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1717                                     vX->GetNrows()-1,vX->GetMatrixArray(),
1718                                     200,-10,10);
1719     fListQAitsPureSA->Add(hNsigmaPPureSA);
1720   }
1721   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
1722                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1723                               vX->GetNrows()-1,vX->GetMatrixArray(),
1724                               300,0,300);
1725   fListQAitsPureSA->Add(hSigPureSA);
1726   
1727   delete vX;  
1728 }
1729
1730 //_____________________________________________________________________________
1731 void AliAnalysisTaskPIDqa::AddTPCHistogramsSignal(TList *sublist, const char *scenario)
1732 {
1733   //
1734   // Create the TPC qa objects: create histograms for the TPC signal for different settings
1735   //
1736
1737   TVectorD *vX=MakeLogBinning(200,.1,30);
1738   Int_t nBinsMult = 38;
1739   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
1740                            120, 140, 160, 180, 200, 
1741                            300, 400, 500, 600, 700, 800, 900, 1000, 
1742                            1200, 1400, 1600, 1800, 2000, 
1743                            2200, 2400, 2600, 2800, 3000, 
1744                            3200, 3400, 3600, 3800, 4000
1745                            };
1746   const Int_t binsEta=110;
1747   Float_t etaMin=-1.1;
1748   Float_t etaMax=1.1;
1749
1750   char signal[4][12]={"std","IROC","OROCmedium","OROClong"};
1751
1752
1753   // TPC signal vs. p for all particles (standard, IROC, OROCmedium, OROClong)
1754   for (Int_t iSig=0; iSig<4; iSig++) {
1755     TH2F *hSigP = new TH2F(Form("hSigP_TPC_%s_%s",signal[iSig],scenario),
1756                              Form("TPC_%s signal (%s) vs. p;p [GeV]; TPC signal [arb. units]",scenario,signal[iSig]),
1757                              vX->GetNrows()-1,vX->GetMatrixArray(),
1758                              300,0,300);
1759     sublist->Add(hSigP);
1760   }
1761
1762   // MIP pions: TPC signal vs. eta
1763   for (Int_t iSig=0; iSig<4; iSig++) {
1764     TH2F *hSigEtaMIPpi = new TH2F(Form("hSigEta_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1765                                   Form("TPC_%s signal (%s) MIPpi vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1766                                   binsEta,etaMin,etaMax,
1767                                   300,0,300);
1768     sublist->Add(hSigEtaMIPpi);
1769   }
1770
1771   // MIP pions: TPC signal vs. multiplicity
1772   for (Int_t iSig=0; iSig<4; iSig++) {
1773     TH2F *hSigMultMPIpi = new TH2F(Form("hSigMult_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1774                                    Form("TPC_%s signal (%s) MIPpi vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1775                                    nBinsMult,xBinsMult,
1776                                    300,0,300);
1777     sublist->Add(hSigMultMPIpi);
1778   }
1779  
1780   // Electrons: TPC signal vs. eta
1781   for (Int_t iSig=0; iSig<4; iSig++) {
1782     TH2F *hSigEtaEle = new TH2F(Form("hSigEta_TPC_%s_%s_Ele",signal[iSig],scenario),
1783                                 Form("TPC_%s signal (%s) electrons vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1784                                 binsEta,etaMin,etaMax,
1785                                 300,0,300);
1786     sublist->Add(hSigEtaEle);
1787   }
1788
1789   // Electrons: TPC signal vs. multiplicity
1790   for (Int_t iSig=0; iSig<4; iSig++) {
1791     TH2F *hSigMultEle = new TH2F(Form("hSigMult_TPC_%s_%s_Ele",signal[iSig],scenario),
1792                                  Form("TPC_%s signal (%s) electrons vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1793                                  nBinsMult,xBinsMult,
1794                                  300,0,300);
1795     sublist->Add(hSigMultEle);
1796   }
1797
1798   delete vX;
1799
1800 }
1801
1802 //_____________________________________________________________________________
1803 void AliAnalysisTaskPIDqa::AddTPCHistogramsNsigma(TList *sublist, const char *scenario, Int_t scnumber)
1804 {
1805   //
1806   // Create the TPC qa objects: create histograms for TPC Nsigma for different settings
1807   //
1808
1809   TVectorD *vX=MakeLogBinning(200,.1,30.);
1810   Int_t nBinsMult = 38;
1811   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
1812                            120, 140, 160, 180, 200, 
1813                            300, 400, 500, 600, 700, 800, 900, 1000, 
1814                            1200, 1400, 1600, 1800, 2000, 
1815                            2200, 2400, 2600, 2800, 3000, 
1816                            3200, 3400, 3600, 3800, 4000
1817                            };
1818   const Int_t binsEta=110;
1819   Float_t etaMin=-1.1;
1820   Float_t etaMax=1.1;
1821
1822   Int_t nSpecies=0;
1823   
1824   if (scnumber == 4) nSpecies=(Int_t)AliPID::kSPECIES;
1825   else nSpecies=(Int_t)AliPID::kSPECIESC; 
1826
1827   // Nsigma vs. p for different particle species
1828   for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1829     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1830                               Form("TPC_%s n#sigma %s vs. p;p [GeV]; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1831                               vX->GetNrows()-1,vX->GetMatrixArray(),
1832                               200,-10,10);
1833     sublist->Add(hNsigmaP);
1834   }
1835
1836   // Nsigma vs. eta for different particle species (only for some scenarios)
1837   if ( scnumber == 1 || scnumber == 4 ) {
1838     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1839       TH2F *hNsigmaEta = new TH2F(Form("hNsigmaEta_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1840                               Form("TPC_%s n#sigma %s vs. eta;#eta; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1841                               binsEta,etaMin,etaMax,
1842                               200,-10,10);
1843       sublist->Add(hNsigmaEta);
1844     }
1845   }
1846
1847   // Nsigma vs. multiplicity for different particle species (only for some scenarios)
1848   if ( scnumber == 1 || scnumber == 4 ) {
1849     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1850       TH2F *hNsigmaMult = new TH2F(Form("hNsigmaMult_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1851                               Form("TPC_%s n#sigma %s vs. mult;multiplicity; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1852                               nBinsMult,xBinsMult,
1853                               200,-10,10);
1854       sublist->Add(hNsigmaMult);
1855     }
1856   }
1857
1858   // - Beginn: Adding histograms for MIP pions and electrons (only for some scenarios) -
1859   if ( scnumber == 0 || scnumber == 2 || scnumber == 3 ) {
1860
1861     // MIP pions: Nsigma vs. eta 
1862     TH2F *hNsigmaEtaMIPpi = new TH2F(Form("hNsigmaEta_TPC_%s_MIPpi",scenario),
1863                               Form("TPC_%s n#sigma MIPpi vs. eta;#eta; n#sigma",scenario),
1864                               binsEta,etaMin,etaMax,
1865                               200,-10,10);
1866     sublist->Add(hNsigmaEtaMIPpi);
1867
1868     // MIP pions: Nsigma vs. multiplicity
1869     TH2F *hNsigmaMultMIPpi = new TH2F(Form("hNsigmaMult_TPC_%s_MIPpi",scenario),
1870                                Form("TPC_%s n#sigma MIPpi vs. mult;multiplicity; n#sigma",scenario),
1871                                nBinsMult,xBinsMult,
1872                                200,-10,10);
1873     sublist->Add(hNsigmaMultMIPpi);
1874
1875     // Electrons: Nsigma vs. eta
1876     TH2F *hNsigmaEtaEle = new TH2F(Form("hNsigmaEta_TPC_%s_Ele",scenario),
1877                               Form("TPC_%s n#sigma electrons vs. eta;#eta; n#sigma",scenario),
1878                               binsEta,etaMin,etaMax,
1879                               200,-10,10);
1880     sublist->Add(hNsigmaEtaEle);
1881
1882     // Electrons: Nsigma vs. multiplicity
1883     TH2F *hNsigmaMultEle = new TH2F(Form("hNsigmaMult_TPC_%s_Ele",scenario),
1884                                Form("TPC_%s n#sigma electrons vs. mult;multiplicity; n#sigma",scenario),
1885                                nBinsMult,xBinsMult,
1886                                200,-10,10);
1887     sublist->Add(hNsigmaMultEle);
1888   } // - End: Adding histograms for MIP pions and electrons
1889
1890   delete vX;
1891
1892 }
1893
1894 //______________________________________________________________________________
1895 void AliAnalysisTaskPIDqa::SetupTPCqa(Bool_t fillMC, Bool_t fill11h, Bool_t fillV0)
1896 {
1897   //
1898   // Create the TPC qa objects
1899   //
1900   
1901   // Set up the multiplicity binning
1902   Int_t nBinsMult = 38;
1903   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
1904                            120, 140, 160, 180, 200, 
1905                            300, 400, 500, 600, 700, 800, 900, 1000, 
1906                            1200, 1400, 1600, 1800, 2000, 
1907                            2200, 2400, 2600, 2800, 3000, 
1908                            3200, 3400, 3600, 3800, 4000
1909                            };
1910
1911
1912   // Create TPC sublists for different scenarios 
1913   // corresponding to available information, 
1914   // e.g. MC or not, special settings for LHC11h
1915
1916   // basic/default scenario, used always
1917   fListQAtpcBasic=new TList;
1918   fListQAtpcBasic->SetOwner();
1919   fListQAtpcBasic->SetName("TPCBasic");
1920   fListQAtpc->Add(fListQAtpcBasic);
1921  
1922   // MC truth scenario: use only MC truth identified particles
1923   // only available for MC
1924   if (fillMC == kTRUE) {
1925     fListQAtpcMCtruth=new TList;
1926     fListQAtpcMCtruth->SetOwner();
1927     fListQAtpcMCtruth->SetName("TPCMCtruth");
1928     fListQAtpc->Add(fListQAtpcMCtruth);
1929   }
1930   
1931   // Hybrid and OROChigh scenarios, 
1932   // special settings only available for PbPb LHC11h data
1933   if (fill11h == kTRUE) {
1934     fListQAtpcHybrid=new TList;
1935     fListQAtpcHybrid->SetOwner();
1936     fListQAtpcHybrid->SetName("TPCHybrid");
1937     fListQAtpc->Add(fListQAtpcHybrid);
1938   
1939     fListQAtpcOROChigh=new TList;
1940     fListQAtpcOROChigh->SetOwner();
1941     fListQAtpcOROChigh->SetName("TPCOROChigh");
1942     fListQAtpc->Add(fListQAtpcOROChigh);
1943   }
1944
1945   // scenario only for V0s, 
1946   // only available for ESDs
1947   if (fillV0 == kTRUE) {
1948     fListQAtpcV0=new TList;
1949     fListQAtpcV0->SetOwner();
1950     fListQAtpcV0->SetName("TPCV0");
1951     fListQAtpc->Add(fListQAtpcV0);
1952   }
1953
1954
1955   // the default ("basic") scenario
1956   AddTPCHistogramsNsigma(fListQAtpcBasic,"Basic",0);
1957   AddTPCHistogramsSignal(fListQAtpcBasic,"Basic");
1958
1959   // only MC truth identified particles
1960   if (fillMC) {
1961     AddTPCHistogramsNsigma(fListQAtpcMCtruth,"MCtruth",1);
1962   }
1963
1964   // the "hybrid" scenario (only for period LHC11h)
1965   if (fill11h) {
1966     AddTPCHistogramsNsigma(fListQAtpcHybrid,"Hybrid",2);
1967   }
1968
1969   // the "OROC high" scenario (only for period LHC11h)
1970   if (fill11h) {
1971     AddTPCHistogramsNsigma(fListQAtpcOROChigh,"OROChigh",3);
1972   }
1973
1974   // only for V0s
1975   if (fillV0) {
1976     AddTPCHistogramsNsigma(fListQAtpcV0,"V0",4);
1977   }
1978  
1979
1980   // Multiplicity distribution --- as check
1981   TH1F *hMult = new TH1F("hMult_TPC",
1982                          "Multiplicity distribution;multiplicity;counts",
1983                          nBinsMult,xBinsMult);
1984   fListQAtpc->Add(hMult);
1985
1986 }
1987
1988 //______________________________________________________________________________
1989 void AliAnalysisTaskPIDqa::SetupTRDqa()
1990 {
1991   //
1992   // Create the TRD qa objects
1993   //
1994   TVectorD *vX=MakeLogBinning(200,.1,30);
1995   for(Int_t itl = 0; itl < 6; ++itl){
1996     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1997       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
1998                               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)),
1999                               vX->GetNrows()-1, vX->GetMatrixArray(),
2000                               100, 0., 1.);
2001       fListQAtrd->Add(hLikeP);
2002     }
2003   }
2004
2005   // === nSigma Values and signal ===
2006   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2007     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_%s",AliPID::ParticleName(ispecie)),
2008                               Form("TRD n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2009                               vX->GetNrows()-1,vX->GetMatrixArray(),
2010                               100,-10,10);
2011     fListQAtrdNsig->Add(hNsigmaP);
2012   }
2013
2014   TH2F *hSig = new TH2F("hSigP_TRD",
2015                         "TRD signal vs. p;p [GeV]; TRD signal [arb. units]",
2016                         vX->GetNrows()-1,vX->GetMatrixArray(),
2017                         100,0,100);
2018   fListQAtrdNsig->Add(hSig);
2019
2020   fListQAtrd->Add(fListQAtrdNsig);
2021
2022   // === Same after 3 sigma in TPC and TOF
2023   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2024     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_TPCTOF_%s",AliPID::ParticleName(ispecie)),
2025                               Form("TRD n#sigma %s vs. p after 3#sigma cut in TPC&TOF;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2026                               vX->GetNrows()-1,vX->GetMatrixArray(),
2027                               100,-10,10);
2028     fListQAtrdNsigTPCTOF->Add(hNsigmaP);
2029   }
2030   
2031   fListQAtrd->Add(fListQAtrdNsigTPCTOF);
2032   
2033   delete vX;
2034 }
2035
2036 //______________________________________________________________________________
2037 void AliAnalysisTaskPIDqa::SetupTOFqa()
2038 {
2039   //
2040   // Create the TOF qa objects
2041   //
2042   
2043   TVectorD *vX=MakeLogBinning(200,.1,30);
2044
2045   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2046     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
2047                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2048                               vX->GetNrows()-1,vX->GetMatrixArray(),
2049                               200,-10,10);
2050     fListQAtof->Add(hNsigmaP);
2051   }
2052
2053   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
2054   fListQAtof->Add(hnSigT0Fill);
2055   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
2056   fListQAtof->Add(hnSigT0T0);
2057   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
2058   fListQAtof->Add(hnSigT0TOF);
2059   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
2060   fListQAtof->Add(hnSigT0Best);
2061   TH1F *hnDeltaPi = new TH1F("hDelta_TOF_Pion","DeltaT (Pion) [0.75-1.25 GeV/c]",50,-500,500);
2062   fListQAtof->Add(hnDeltaPi);
2063   
2064   TH2F *hSig = new TH2F("hSigP_TOF",
2065                         "TOF signal vs. p;p [GeV]; TOF signal [ns]",
2066                         vX->GetNrows()-1,vX->GetMatrixArray(),
2067                         300,0,30);
2068
2069   delete vX;
2070   
2071   fListQAtof->Add(hSig);
2072
2073   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
2074   fListQAtof->Add(hStartTimeMaskTOF);
2075   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
2076   fListQAtof->Add(hStartTimeResTOF);
2077
2078   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",100,0,100);
2079   fListQAtof->Add(hnTracksAtTOF);
2080   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",100,0,100);
2081   fListQAtof->Add(hT0MakerEff);
2082
2083   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
2084   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2085   fListQAtof->Add(hStartTimeAT0);
2086   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2087   fListQAtof->Add(hStartTimeCT0);
2088   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2089   fListQAtof->Add(hStartTimeACT0);
2090 }
2091
2092
2093 //______________________________________________________________________________
2094 void AliAnalysisTaskPIDqa::SetupT0qa()
2095 {
2096   //
2097   // Create the T0 qa objects
2098   //
2099   
2100   // these are similar to plots inside TOFqa, but these are for all events
2101   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2102   fListQAt0->Add(hStartTimeAT0);
2103   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2104   fListQAt0->Add(hStartTimeCT0);
2105   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2106   fListQAt0->Add(hStartTimeACT0);
2107
2108   TH1F *hnTracksAtT0 = new TH1F("hnTracksAt_T0","Tracks for events selected for T0",100,0,100);
2109   fListQAt0->Add(hnTracksAtT0);
2110   TH1F *hT0AEff = new TH1F("hT0AEff","Events with T0A vs nTracks",100,0,100);
2111   fListQAt0->Add(hT0AEff);
2112   TH1F *hT0CEff = new TH1F("hT0CEff","Events with T0C vs nTracks",100,0,100);
2113   fListQAt0->Add(hT0CEff);
2114   TH1F *hT0AndEff = new TH1F("hT0AndEff","Events with T0AC (AND) vs nTracks",100,0,100);
2115   fListQAt0->Add(hT0AndEff);
2116   TH1F *hT0OrEff = new TH1F("hT0OrEff","Events with T0AC (OR) vs nTracks",100,0,100);
2117   fListQAt0->Add(hT0OrEff);
2118
2119
2120 }
2121
2122 //______________________________________________________________________________
2123 void AliAnalysisTaskPIDqa::SetupEMCALqa()
2124 {
2125   //
2126   // Create the EMCAL qa objects
2127   //
2128
2129   TVectorD *vX=MakeLogBinning(200,.1,30);
2130   
2131   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
2132                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
2133                              vX->GetNrows()-1,vX->GetMatrixArray(),
2134                              200,-10,10);
2135   fListQAemcal->Add(hNsigmaPt);  
2136   
2137
2138   TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele",
2139                         "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2140                         vX->GetNrows()-1,vX->GetMatrixArray(),
2141                         200,0,2);
2142   fListQAemcal->Add(hSigPtEle);
2143
2144   TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions",
2145                         "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2146                         vX->GetNrows()-1,vX->GetMatrixArray(),
2147                         200,0,2);
2148   fListQAemcal->Add(hSigPtPions);
2149
2150   TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons",
2151                         "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2152                         vX->GetNrows()-1,vX->GetMatrixArray(),
2153                         200,0,2);
2154   fListQAemcal->Add(hSigPtProtons);
2155
2156   TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons",
2157                         "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2158                         vX->GetNrows()-1,vX->GetMatrixArray(),
2159                         200,0,2);
2160   fListQAemcal->Add(hSigPtAntiProtons);
2161
2162   delete vX;  
2163 }
2164
2165 //______________________________________________________________________________
2166 void AliAnalysisTaskPIDqa::SetupHMPIDqa()
2167 {
2168   //
2169   // Create the HMPID qa objects
2170   //
2171
2172   TVectorD *vX=MakeLogBinning(200,.1,30);
2173
2174   // nSigmas
2175   Int_t nhists=0;
2176   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2177     if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
2178     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_HMPID_%s",AliPID::ParticleName(ispecie)),
2179                               Form("HMPID n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2180                               vX->GetNrows()-1,vX->GetMatrixArray(),
2181                               200,-10,10);
2182     fListQAhmpid->AddAt(hNsigmaP, nhists);
2183     ++nhists;
2184   }
2185   
2186   // cherenkov angle
2187   TH2F *hCkovAnglevsMom   = new TH2F("hCkovAnglevsMom",  "Cherenkov angle vs momentum",
2188                                      vX->GetNrows()-1,vX->GetMatrixArray(),
2189                                      500,0,1);
2190   fListQAhmpid->AddAt(hCkovAnglevsMom,nhists);
2191   
2192   delete vX;
2193 }
2194
2195 //______________________________________________________________________________
2196 void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
2197 {
2198   //
2199   // Create the HMPID qa objects
2200   //
2201   
2202   TH2F *hCkovAnglevsMomPion   = new TH2F("hCkovAnglevsMom_pion",  "Cherenkov angle vs momentum for pions",500,0,5.,500,0,1);
2203   fListQAtofhmpid->Add(hCkovAnglevsMomPion);
2204   
2205   TH2F *hCkovAnglevsMomKaon   = new TH2F("hCkovAnglevsMom_kaon",  "Cherenkov angle vs momentum for kaons",500,0,5.,500,0,1);
2206   fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
2207   
2208   TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momentum for protons",500,0,5.,500,0,1);
2209   fListQAtofhmpid->Add(hCkovAnglevsMomProton);
2210   
2211   
2212 }  
2213
2214 //______________________________________________________________________________
2215 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
2216 {
2217   //
2218   // Create the qa objects for TPC + TOF combination
2219   //
2220   
2221   TVectorD *vX=MakeLogBinning(200,.1,30);
2222
2223   //TPC signals after TOF cut
2224   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2225     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
2226                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2227                               vX->GetNrows()-1,vX->GetMatrixArray(),
2228                               200,-10,10);
2229     fListQAtpctof->Add(hNsigmaP);
2230   }
2231
2232   //TOF signals after TPC cut
2233   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2234     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2235                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2236                               vX->GetNrows()-1,vX->GetMatrixArray(),
2237                               200,-10,10);
2238     fListQAtpctof->Add(hNsigmaP);
2239   }
2240
2241   //EMCAL signal after TOF and TPC cut
2242   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2243     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2244                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
2245                             vX->GetNrows()-1,vX->GetMatrixArray(),
2246                             200,0,2);
2247     fListQAtpctof->Add(heopPt);
2248   }
2249
2250   delete vX;
2251 }
2252 //______________________________________________________________________________
2253 void AliAnalysisTaskPIDqa::SetupV0qa()
2254 {
2255   //
2256   // Create the qa objects for V0 Kine cuts
2257   //
2258   
2259   TH2F *hArmenteros  = new TH2F("hArmenteros",  "Armenteros plot",200,-1.,1.,200,0.,0.4);
2260   fListQAV0->Add(hArmenteros);
2261  
2262 }
2263
2264 //_____________________________________________________________________________
2265 void AliAnalysisTaskPIDqa::SetupQAinfo(){
2266   //
2267   // Setup the info of QA objects
2268   //
2269
2270   TObjArray *arr=new TObjArray;
2271   arr->SetName("TPC_info");
2272   fListQAinfo->Add(arr);
2273 }
2274
2275 //______________________________________________________________________________
2276 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2277 {
2278   //
2279   // Make logarithmic binning
2280   // the user has to delete the array afterwards!!!
2281   //
2282   
2283   //check limits
2284   if (xmin<1e-20 || xmax<1e-20){
2285     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
2286     return MakeLinBinning(nbinsX, xmin, xmax);
2287   }
2288   if (xmax<xmin){
2289     Double_t tmp=xmin;
2290     xmin=xmax;
2291     xmax=tmp;
2292   }
2293   TVectorD *binLim=new TVectorD(nbinsX+1);
2294   Double_t first=xmin;
2295   Double_t last=xmax;
2296   Double_t expMax=TMath::Log(last/first);
2297   for (Int_t i=0; i<nbinsX+1; ++i){
2298     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
2299   }
2300   return binLim;
2301 }
2302
2303 //______________________________________________________________________________
2304 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2305 {
2306   //
2307   // Make linear binning
2308   // the user has to delete the array afterwards!!!
2309   //
2310   if (xmax<xmin){
2311     Double_t tmp=xmin;
2312     xmin=xmax;
2313     xmax=tmp;
2314   }
2315   TVectorD *binLim=new TVectorD(nbinsX+1);
2316   Double_t first=xmin;
2317   Double_t last=xmax;
2318   Double_t binWidth=(last-first)/nbinsX;
2319   for (Int_t i=0; i<nbinsX+1; ++i){
2320     (*binLim)[i]=first+binWidth*(Double_t)i;
2321   }
2322   return binLim;
2323 }
2324
2325 //_____________________________________________________________________________
2326 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
2327 {
2328   //
2329   // Make arbitrary binning, bins separated by a ','
2330   //
2331   TString limits(bins);
2332   if (limits.IsNull()){
2333     AliError("Bin Limit string is empty, cannot add the variable");
2334     return 0x0;
2335   }
2336   
2337   TObjArray *arr=limits.Tokenize(",");
2338   Int_t nLimits=arr->GetEntries();
2339   if (nLimits<2){
2340     AliError("Need at leas 2 bin limits, cannot add the variable");
2341     delete arr;
2342     return 0x0;
2343   }
2344   
2345   TVectorD *binLimits=new TVectorD(nLimits);
2346   for (Int_t iLim=0; iLim<nLimits; ++iLim){
2347     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
2348   }
2349   
2350   delete arr;
2351   return binLimits;
2352 }
2353