Fix in the last caall to CleanOwnPrimaryVertex
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskPIDqa.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 /* $Id: AliAnalysisTaskPIDqa.cxx 43811 2010-09-23 14:13:31Z wiechula $ */\r
17 #include <TList.h>\r
18 #include <TVectorD.h>\r
19 #include <TObjArray.h>\r
20 #include <TH2.h>\r
21 #include <TFile.h>\r
22 #include <TPRegexp.h>\r
23 #include <TChain.h>\r
24 #include <TF1.h>\r
25 #include <TSpline.h>\r
26 \r
27 #include <AliAnalysisManager.h>\r
28 #include <AliInputEventHandler.h>\r
29 #include <AliVEventHandler.h>\r
30 #include <AliVEvent.h>\r
31 #include <AliVParticle.h>\r
32 #include <AliVTrack.h>\r
33 #include <AliLog.h>\r
34 #include <AliPID.h>\r
35 #include <AliPIDResponse.h>\r
36 #include <AliITSPIDResponse.h>\r
37 #include <AliTPCPIDResponse.h>\r
38 #include <AliTRDPIDResponse.h>\r
39 #include <AliTOFPIDResponse.h>\r
40 \r
41 #include "AliAnalysisTaskPIDqa.h"\r
42 \r
43 \r
44 ClassImp(AliAnalysisTaskPIDqa)\r
45 \r
46 //______________________________________________________________________________\r
47 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():\r
48 AliAnalysisTaskSE(),\r
49 fPIDResponse(0x0),\r
50 fListQA(0x0),\r
51 fListQAits(0x0),\r
52 fListQAitsSA(0x0),\r
53 fListQAitsPureSA(0x0),\r
54 fListQAtpc(0x0),\r
55 fListQAtrd(0x0),\r
56 fListQAtof(0x0),\r
57 fListQAemcal(0x0),\r
58 fListQAtpctof(0x0)\r
59 {\r
60   //\r
61   // Dummy constructor\r
62   //\r
63 }\r
64 \r
65 //______________________________________________________________________________\r
66 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):\r
67 AliAnalysisTaskSE(name),\r
68 fPIDResponse(0x0),\r
69 fListQA(0x0),\r
70 fListQAits(0x0),\r
71 fListQAitsSA(0x0),\r
72 fListQAitsPureSA(0x0),\r
73 fListQAtpc(0x0),\r
74 fListQAtrd(0x0),\r
75 fListQAtof(0x0),\r
76 fListQAemcal(0x0),\r
77 fListQAtpctof(0x0)\r
78 {\r
79   //\r
80   // Default constructor\r
81   //\r
82   DefineInput(0,TChain::Class());\r
83   DefineOutput(1,TList::Class());\r
84 }\r
85 \r
86 //______________________________________________________________________________\r
87 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()\r
88 {\r
89   //\r
90   // Destructor\r
91   //\r
92 \r
93 }\r
94 \r
95 //______________________________________________________________________________\r
96 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()\r
97 {\r
98   //\r
99   // Create the output QA objects\r
100   //\r
101 \r
102   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);\r
103 \r
104   //input hander\r
105   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();\r
106   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());\r
107   if (!inputHandler) AliFatal("Input handler needed");\r
108 \r
109   //pid response object\r
110   fPIDResponse=inputHandler->GetPIDResponse();\r
111   if (!fPIDResponse) AliError("PIDResponse object was not created");\r
112   \r
113   //\r
114   fListQA=new TList;\r
115   fListQA->SetOwner();\r
116   \r
117   fListQAits=new TList;\r
118   fListQAits->SetOwner();\r
119   fListQAits->SetName("ITS");\r
120 \r
121   fListQAitsSA=new TList;\r
122   fListQAitsSA->SetOwner();\r
123   fListQAitsSA->SetName("ITS");\r
124 \r
125   fListQAitsPureSA=new TList;\r
126   fListQAitsPureSA->SetOwner();\r
127   fListQAitsPureSA->SetName("ITS");\r
128 \r
129   fListQAtpc=new TList;\r
130   fListQAtpc->SetOwner();\r
131   fListQAtpc->SetName("TPC");\r
132   \r
133   fListQAtrd=new TList;\r
134   fListQAtrd->SetOwner();\r
135   fListQAtrd->SetName("TRD");\r
136   \r
137   fListQAtof=new TList;\r
138   fListQAtof->SetOwner();\r
139   fListQAtof->SetName("TOF");\r
140   \r
141   fListQAemcal=new TList;\r
142   fListQAemcal->SetOwner();\r
143   fListQAemcal->SetName("EMCAL");\r
144 \r
145   fListQAtpctof=new TList;\r
146   fListQAtpctof->SetOwner();\r
147   fListQAtpctof->SetName("TPC_TOF");\r
148 \r
149   fListQA->Add(fListQAits);\r
150   fListQA->Add(fListQAitsSA);\r
151   fListQA->Add(fListQAitsPureSA);\r
152   fListQA->Add(fListQAtpc);\r
153   fListQA->Add(fListQAtrd);\r
154   fListQA->Add(fListQAtof);\r
155   fListQA->Add(fListQAemcal);\r
156   fListQA->Add(fListQAtpctof);\r
157 \r
158   SetupITSqa();\r
159   SetupTPCqa();\r
160   SetupTRDqa();\r
161   SetupTOFqa();\r
162   SetupEMCALqa();\r
163   SetupTPCTOFqa();\r
164 \r
165   PostData(1,fListQA);\r
166 }\r
167 \r
168 \r
169 //______________________________________________________________________________\r
170 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)\r
171 {\r
172   //\r
173   // Setup the PID response functions and fill the QA histograms\r
174   //\r
175 \r
176   AliVEvent *event=InputEvent();\r
177   if (!event||!fPIDResponse) return;\r
178 \r
179   \r
180   FillITSqa();\r
181   FillTPCqa();\r
182   FillTRDqa();\r
183   FillTOFqa();\r
184   FillEMCALqa();\r
185   FillTPCTOFqa();\r
186 \r
187   PostData(1,fListQA);\r
188 }\r
189 \r
190 //______________________________________________________________________________\r
191 void AliAnalysisTaskPIDqa::FillITSqa()\r
192 {\r
193   //\r
194   // Fill PID qa histograms for the ITS\r
195   //\r
196 \r
197   AliVEvent *event=InputEvent();\r
198   \r
199   Int_t ntracks=event->GetNumberOfTracks();\r
200   for(Int_t itrack = 0; itrack < ntracks; itrack++){\r
201     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);\r
202     ULong_t status=track->GetStatus();\r
203     // not that nice. status bits not in virtual interface\r
204     // ITS refit + ITS pid selection
205     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
206            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
207     Double_t mom=track->P();\r
208     \r
209     TList *theList = 0x0;
210     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
211       //ITS+TPC tracks
212       theList=fListQAits;
213     }else{
214       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){ 
215         //ITS Standalone tracks
216         theList=fListQAitsSA;
217       }else{
218         //ITS Pure Standalone tracks
219         theList=fListQAitsPureSA;
220       }
221     }
222     
223     
224     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
225       TH2 *h=(TH2*)theList->At(ispecie);
226       if (!h) continue;\r
227       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);\r
228       h->Fill(mom,nSigma);\r
229     }\r
230     TH2 *h=(TH2*)theList->At(AliPID::kSPECIES);
231     if (h) {\r
232       Double_t sig=track->GetITSsignal();\r
233       h->Fill(mom,sig);\r
234     }\r
235   }\r
236 }\r
237 \r
238 //______________________________________________________________________________\r
239 void AliAnalysisTaskPIDqa::FillTPCqa()\r
240 {\r
241   //\r
242   // Fill PID qa histograms for the TPC\r
243   //\r
244   \r
245   AliVEvent *event=InputEvent();\r
246   \r
247   Int_t ntracks=event->GetNumberOfTracks();\r
248   for(Int_t itrack = 0; itrack < ntracks; itrack++){\r
249     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);\r
250     \r
251     //\r
252     //basic track cuts\r
253     //\r
254     ULong_t status=track->GetStatus();\r
255     // not that nice. status bits not in virtual interface\r
256     // TPC refit + ITS refit + TPC pid\r
257     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||\r
258         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||\r
259         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) ) continue;\r
260     \r
261     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);\r
262     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;\r
263     if (track->GetTPCNclsF()>0) {\r
264       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();\r
265     }\r
266     \r
267     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;\r
268     \r
269     Double_t mom=track->GetTPCmomentum();\r
270     \r
271     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
272       TH2 *h=(TH2*)fListQAtpc->At(ispecie);\r
273       if (!h) continue;\r
274       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);\r
275       h->Fill(mom,nSigma);\r
276     }\r
277     \r
278     TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);\r
279     if (h) {\r
280       Double_t sig=track->GetTPCsignal();\r
281       h->Fill(mom,sig);\r
282     }\r
283   }\r
284 }\r
285 \r
286 //______________________________________________________________________________\r
287 void AliAnalysisTaskPIDqa::FillTRDqa()\r
288 {\r
289   //\r
290   // Fill PID qa histograms for the TRD\r
291   //\r
292 \r
293 }\r
294 \r
295 //______________________________________________________________________________\r
296 void AliAnalysisTaskPIDqa::FillTOFqa()\r
297 {\r
298   //\r
299   // Fill PID qa histograms for the TOF\r
300   //\r
301 \r
302   AliVEvent *event=InputEvent();\r
303   \r
304   Int_t ntracks=event->GetNumberOfTracks();\r
305   for(Int_t itrack = 0; itrack < ntracks; itrack++){\r
306     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);\r
307     \r
308     //\r
309     //basic track cuts\r
310     //\r
311     ULong_t status=track->GetStatus();\r
312     // not that nice. status bits not in virtual interface\r
313     // TPC refit + ITS refit +\r
314     // TOF out + TOFpid +\r
315     // kTIME\r
316     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||\r
317         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||\r
318         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||\r
319         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||\r
320         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;\r
321     \r
322     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);\r
323     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;\r
324     if (track->GetTPCNclsF()>0) {\r
325       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();\r
326     }\r
327     \r
328     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;\r
329     \r
330     \r
331     Double_t mom=track->P();\r
332     \r
333     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
334       //TOF nSigma\r
335       TH2 *h=(TH2*)fListQAtof->At(ispecie);\r
336       if (!h) continue;\r
337       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);\r
338       h->Fill(mom,nSigma);\r
339     }\r
340     \r
341     TH2 *h=(TH2*)fListQAtof->At(AliPID::kSPECIES);\r
342     if (h) {\r
343       Double_t sig=track->GetTOFsignal();\r
344       h->Fill(mom,sig);\r
345     }\r
346     \r
347   }\r
348 }\r
349 \r
350 //______________________________________________________________________________\r
351 void AliAnalysisTaskPIDqa::FillEMCALqa()\r
352 {\r
353   //\r
354   // Fill PID qa histograms for the EMCAL\r
355   //\r
356 \r
357   AliVEvent *event=InputEvent();\r
358   \r
359   Int_t ntracks=event->GetNumberOfTracks();\r
360   for(Int_t itrack = 0; itrack < ntracks; itrack++){\r
361     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);\r
362     \r
363     //\r
364     //basic track cuts\r
365     //\r
366     ULong_t status=track->GetStatus();\r
367     // not that nice. status bits not in virtual interface\r
368     // TPC refit + ITS refit +\r
369     // TOF out + TOFpid +\r
370     // kTIME\r
371     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;\r
372 \r
373     Double_t pt=track->Pt();\r
374    \r
375     //EMCAL nSigma (only for electrons at the moment)\r
376     TH2 *h=(TH2*)fListQAemcal->At(0);\r
377     if (!h) continue;\r
378     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);\r
379     h->Fill(pt,nSigma);\r
380     \r
381     //EMCAL signal (E/p vs. pT)\r
382     h=(TH2*)fListQAemcal->At(1);\r
383     if (h) {\r
384 \r
385       Int_t nMatchClus = track->GetEMCALcluster();\r
386       Double_t mom     = track->P();\r
387       Double_t eop     = -1.;\r
388 \r
389       if(nMatchClus > -1){\r
390     \r
391         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);\r
392     \r
393         if(matchedClus){\r
394 \r
395           // matched cluster is EMCAL\r
396           if(matchedClus->IsEMCAL()){\r
397       \r
398             Double_t fClsE       = matchedClus->E();\r
399             eop                  = fClsE/mom;\r
400 \r
401             h->Fill(pt,eop);\r
402             \r
403           }\r
404         }\r
405       }\r
406       else{\r
407         Printf("status status = AliVTrack::kEMCALmatch, BUT no matched cluster!");\r
408       }\r
409     }\r
410   }\r
411 }\r
412 \r
413 //______________________________________________________________________________\r
414 void AliAnalysisTaskPIDqa::FillTPCTOFqa()\r
415 {\r
416   //\r
417   // Fill PID qa histograms for the TOF\r
418   //   Here also the TPC histograms after TOF selection are filled\r
419   //\r
420 \r
421   AliVEvent *event=InputEvent();\r
422 \r
423   Int_t ntracks=event->GetNumberOfTracks();\r
424   for(Int_t itrack = 0; itrack < ntracks; itrack++){\r
425     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);\r
426 \r
427     //\r
428     //basic track cuts\r
429     //\r
430     ULong_t status=track->GetStatus();\r
431     // not that nice. status bits not in virtual interface\r
432     // TPC refit + ITS refit +\r
433     // TOF out + TOFpid +\r
434     // kTIME\r
435     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||\r
436         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||\r
437         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) ||\r
438         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||\r
439         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||\r
440         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;\r
441 \r
442     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);\r
443     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;\r
444     if (track->GetTPCNclsF()>0) {\r
445       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();\r
446     }\r
447 \r
448     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;\r
449 \r
450 \r
451     Double_t mom=track->P();\r
452     Double_t momTPC=track->GetTPCmomentum();\r
453 \r
454     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
455       //TOF nSigma\r
456       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);\r
457       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);\r
458 \r
459       //TPC after TOF cut\r
460       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);\r
461       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);\r
462 \r
463       //TOF after TPC cut\r
464       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);\r
465       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);\r
466 \r
467       //EMCAL after TOF and TPC cut\r
468       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);\r
469       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){\r
470 \r
471         Int_t nMatchClus = track->GetEMCALcluster();\r
472         Double_t pt      = track->Pt();\r
473         Double_t eop     = -1.;\r
474         \r
475         if(nMatchClus > -1){\r
476           \r
477           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);\r
478           \r
479           if(matchedClus){\r
480             \r
481             // matched cluster is EMCAL\r
482             if(matchedClus->IsEMCAL()){\r
483               \r
484               Double_t fClsE       = matchedClus->E();\r
485               eop                  = fClsE/mom;\r
486 \r
487               h->Fill(pt,eop);\r
488  \r
489               \r
490             }\r
491           }\r
492         }\r
493       }\r
494     }\r
495   }\r
496 }\r
497 \r
498 //______________________________________________________________________________\r
499 void AliAnalysisTaskPIDqa::SetupITSqa()\r
500 {\r
501   //\r
502   // Create the ITS qa objects\r
503   //\r
504   \r
505   TVectorD *vX=MakeLogBinning(200,.1,30);\r
506   \r
507   //ITS+TPC tracks
508   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
509     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),\r
510                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
511                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
512                               200,-10,10);\r
513     fListQAits->Add(hNsigmaP);\r
514   }\r
515   TH2F *hSig = new TH2F("hSigP_ITS",\r
516                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",\r
517                         vX->GetNrows()-1,vX->GetMatrixArray(),\r
518                         300,0,300);\r
519   fListQAits->Add(hSig);\r
520 \r
521   //ITS Standalone tracks
522   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
523     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
524                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
525                                 vX->GetNrows()-1,vX->GetMatrixArray(),
526                                 200,-10,10);
527     fListQAitsSA->Add(hNsigmaPSA);
528   }
529   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
530                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
531                           vX->GetNrows()-1,vX->GetMatrixArray(),
532                           300,0,300);
533   fListQAitsSA->Add(hSigSA);
534   
535   //ITS Pure Standalone tracks
536   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
537     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
538                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
539                                     vX->GetNrows()-1,vX->GetMatrixArray(),
540                                     200,-10,10);
541     fListQAitsPureSA->Add(hNsigmaPPureSA);
542   }
543   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
544                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
545                               vX->GetNrows()-1,vX->GetMatrixArray(),
546                               300,0,300);
547   fListQAitsPureSA->Add(hSigPureSA);
548   
549   delete vX;  \r
550 }\r
551 \r
552 //______________________________________________________________________________\r
553 void AliAnalysisTaskPIDqa::SetupTPCqa()\r
554 {\r
555   //\r
556   // Create the TPC qa objects\r
557   //\r
558   \r
559   TVectorD *vX=MakeLogBinning(200,.1,30);\r
560   \r
561   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
562     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),\r
563                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
564                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
565                               200,-10,10);\r
566     fListQAtpc->Add(hNsigmaP);\r
567   }\r
568   \r
569   \r
570   TH2F *hSig = new TH2F("hSigP_TPC",\r
571                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",\r
572                         vX->GetNrows()-1,vX->GetMatrixArray(),\r
573                         300,0,300);\r
574   fListQAtpc->Add(hSig);\r
575 \r
576   delete vX;  \r
577 }\r
578 \r
579 //______________________________________________________________________________\r
580 void AliAnalysisTaskPIDqa::SetupTRDqa()\r
581 {\r
582   //\r
583   // Create the TRD qa objects\r
584   //\r
585   \r
586 }\r
587 \r
588 //______________________________________________________________________________\r
589 void AliAnalysisTaskPIDqa::SetupTOFqa()\r
590 {\r
591   //\r
592   // Create the TOF qa objects\r
593   //\r
594   \r
595   TVectorD *vX=MakeLogBinning(200,.1,30);\r
596   \r
597   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
598     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),\r
599                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
600                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
601                               200,-10,10);\r
602     fListQAtof->Add(hNsigmaP);\r
603   }\r
604   \r
605   \r
606   TH2F *hSig = new TH2F("hSigP_TOF",\r
607                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",\r
608                         vX->GetNrows()-1,vX->GetMatrixArray(),\r
609                         300,0,300);\r
610   fListQAtof->Add(hSig);\r
611 \r
612   delete vX;  \r
613 }\r
614 \r
615 //______________________________________________________________________________\r
616 void AliAnalysisTaskPIDqa::SetupEMCALqa()\r
617 {\r
618   //\r
619   // Create the EMCAL qa objects\r
620   //\r
621 \r
622   TVectorD *vX=MakeLogBinning(200,.1,30);\r
623   \r
624   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),\r
625                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),\r
626                              vX->GetNrows()-1,vX->GetMatrixArray(),\r
627                              200,-10,10);\r
628   fListQAemcal->Add(hNsigmaPt);  \r
629   \r
630   TH2F *hSigPt = new TH2F("hSigPt_EMCAL",\r
631                         "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",\r
632                         vX->GetNrows()-1,vX->GetMatrixArray(),\r
633                         200,0,2);\r
634   fListQAemcal->Add(hSigPt);\r
635 \r
636   delete vX;  \r
637 }\r
638 \r
639 //______________________________________________________________________________\r
640 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()\r
641 {\r
642   //\r
643   // Create the qa objects for TPC + TOF combination\r
644   //\r
645   \r
646   TVectorD *vX=MakeLogBinning(200,.1,30);\r
647 \r
648   //TPC signals after TOF cut\r
649   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
650     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),\r
651                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
652                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
653                               200,-10,10);\r
654     fListQAtpctof->Add(hNsigmaP);\r
655   }\r
656 \r
657   //TOF signals after TPC cut\r
658   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
659     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),\r
660                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
661                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
662                               200,-10,10);\r
663     fListQAtpctof->Add(hNsigmaP);\r
664   }\r
665 \r
666   //EMCAL signal after TOF and TPC cut\r
667   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
668     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),\r
669                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),\r
670                             vX->GetNrows()-1,vX->GetMatrixArray(),\r
671                             200,0,2);\r
672     fListQAtpctof->Add(heopPt);\r
673   }\r
674 \r
675   delete vX;\r
676 }\r
677 \r
678 //______________________________________________________________________________\r
679 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)\r
680 {\r
681   //\r
682   // Make logarithmic binning\r
683   // the user has to delete the array afterwards!!!\r
684   //\r
685   \r
686   //check limits\r
687   if (xmin<1e-20 || xmax<1e-20){\r
688     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");\r
689     return MakeLinBinning(nbinsX, xmin, xmax);\r
690   }\r
691   if (xmax<xmin){\r
692     Double_t tmp=xmin;\r
693     xmin=xmax;\r
694     xmax=tmp;\r
695   }\r
696   TVectorD *binLim=new TVectorD(nbinsX+1);\r
697   Double_t first=xmin;\r
698   Double_t last=xmax;\r
699   Double_t expMax=TMath::Log(last/first);\r
700   for (Int_t i=0; i<nbinsX+1; ++i){\r
701     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);\r
702   }\r
703   return binLim;\r
704 }\r
705 \r
706 //______________________________________________________________________________\r
707 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)\r
708 {\r
709   //\r
710   // Make linear binning\r
711   // the user has to delete the array afterwards!!!\r
712   //\r
713   if (xmax<xmin){\r
714     Double_t tmp=xmin;\r
715     xmin=xmax;\r
716     xmax=tmp;\r
717   }\r
718   TVectorD *binLim=new TVectorD(nbinsX+1);\r
719   Double_t first=xmin;\r
720   Double_t last=xmax;\r
721   Double_t binWidth=(last-first)/nbinsX;\r
722   for (Int_t i=0; i<nbinsX+1; ++i){\r
723     (*binLim)[i]=first+binWidth*(Double_t)i;\r
724   }\r
725   return binLim;\r
726 }\r
727 \r
728 //_____________________________________________________________________________\r
729 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)\r
730 {\r
731   //\r
732   // Make arbitrary binning, bins separated by a ','\r
733   //\r
734   TString limits(bins);\r
735   if (limits.IsNull()){\r
736     AliError("Bin Limit string is empty, cannot add the variable");\r
737     return 0x0;\r
738   }\r
739   \r
740   TObjArray *arr=limits.Tokenize(",");\r
741   Int_t nLimits=arr->GetEntries();\r
742   if (nLimits<2){\r
743     AliError("Need at leas 2 bin limits, cannot add the variable");\r
744     delete arr;\r
745     return 0x0;\r
746   }\r
747   \r
748   TVectorD *binLimits=new TVectorD(nLimits);\r
749   for (Int_t iLim=0; iLim<nLimits; ++iLim){\r
750     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();\r
751   }\r
752   \r
753   delete arr;\r
754   return binLimits;\r
755 }\r
756