ownership
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliAnalysisTaskPIDflowQA.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 pures 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 pure. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // AliAnalysisTaskPIDflowQA:
19 // QA for pid
20 //
21 //origin: Marek Chojnacki, Marek.Chojnacki@cern.ch
22 //modified: Mikolaj Krzewicki, Mikolaj.Krzewicki@cern.ch
23
24 #include "AliAnalysisTaskPIDflowQA.h"
25 #include "AliAnalysisManager.h"
26 #include "AliESDEvent.h"
27 #include "AliInputEventHandler.h"
28 #include "AliStack.h"
29 #include "AliMCEvent.h"
30 #include "TH2F.h"
31 #include "TProfile.h"
32 #include "TMath.h"
33 #include "AliVEvent.h"
34 #include "AliESDtrackCuts.h"
35 #include "AliPID.h"
36 #include "AliCDBManager.h"
37 #include "AliFlowEventCuts.h"
38 #include "AliFlowTrackCuts.h"
39
40 ClassImp( AliAnalysisTaskPIDflowQA)
41
42 //________________________________________________________________________
43 AliAnalysisTaskPIDflowQA:: AliAnalysisTaskPIDflowQA():
44   AliAnalysisTaskSE("AliAnalysisTaskPIDflowQA"),
45   fESD(NULL),
46   fCuts(NULL),
47   fEventCuts(NULL),
48   fESDpid(NULL),
49   fMC(kFALSE),
50   fITSsignal(NULL),
51   fTPCsignal(NULL),
52   fTOFsignal(NULL),
53   fITSsignalpip(NULL),
54   fTPCsignalpip(NULL),
55   fTOFsignalpip(NULL),
56   fITSsignalKp(NULL),
57   fTPCsignalKp(NULL),
58   fTOFsignalKp(NULL),
59   fITSsignalpp(NULL),
60   fTPCsignalpp(NULL),
61   fTOFsignalpp(NULL),
62   fITSsignalpiMCp(NULL),
63   fTPCsignalpiMCp(NULL),
64   fTOFsignalpiMCp(NULL),
65   fITSsignalKMCp(NULL),
66   fTPCsignalKMCp(NULL),
67   fTOFsignalKMCp(NULL),
68   fITSsignalpMCp(NULL),
69   fTPCsignalpMCp(NULL),
70   fTOFsignalpMCp(NULL),
71   fTOFsignalBeta(NULL),
72   fTOFsignalPiBeta(NULL),
73   fTOFsignalKBeta(NULL),
74   fTOFsignalPBeta(NULL),
75   fPvsPt(NULL),
76   fMeanPvsP(NULL),
77   fTPCvsGlobalMult(NULL),
78   fStandardGlobalCuts(NULL),
79   fStandardTPCCuts(NULL),
80   fOutputList(NULL)
81 {
82   //def ctor
83 }
84
85 //________________________________________________________________________
86 AliAnalysisTaskPIDflowQA:: AliAnalysisTaskPIDflowQA(const char *name):
87   AliAnalysisTaskSE(name),
88   fESD(NULL),
89   fCuts(NULL),
90   fEventCuts(NULL),
91   fESDpid(NULL),
92   fMC(kFALSE),
93   fITSsignal(NULL),
94   fTPCsignal(NULL),
95   fTOFsignal(NULL),
96   fITSsignalpip(NULL),
97   fTPCsignalpip(NULL),
98   fTOFsignalpip(NULL),
99   fITSsignalKp(NULL),
100   fTPCsignalKp(NULL),
101   fTOFsignalKp(NULL),
102   fITSsignalpp(NULL),
103   fTPCsignalpp(NULL),
104   fTOFsignalpp(NULL),
105   fITSsignalpiMCp(NULL),
106   fTPCsignalpiMCp(NULL),
107   fTOFsignalpiMCp(NULL),
108   fITSsignalKMCp(NULL),
109   fTPCsignalKMCp(NULL),
110   fTOFsignalKMCp(NULL),
111   fITSsignalpMCp(NULL),
112   fTPCsignalpMCp(NULL),
113   fTOFsignalpMCp(NULL),
114   fTOFsignalBeta(NULL),
115   fTOFsignalPiBeta(NULL),
116   fTOFsignalKBeta(NULL),
117   fTOFsignalPBeta(NULL),
118   fPvsPt(NULL),
119   fMeanPvsP(NULL),
120   fTPCvsGlobalMult(NULL),
121   fStandardGlobalCuts(NULL),
122   fStandardTPCCuts(NULL),
123   fOutputList(NULL)
124 {
125   //Constructor
126   fESDpid=new AliESDpid();
127
128   //old
129   fESDpid->GetTPCResponse().SetBetheBlochParameters(0.0283086,
130       2.63394e+01,
131       5.04114e-11,
132       2.12543e+00,
133       4.88663e+00 );
134   //new
135   //fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28949/50.,
136   //                                                  2.74095e+01,
137   //                                                  TMath::Exp(-3.21763e+01),
138   //                                                  2.44026,
139   //                                                  6.58800);
140
141   DefineOutput(1, TList::Class());
142 }
143
144 //________________________________________________________________________
145 void  AliAnalysisTaskPIDflowQA::UserCreateOutputObjects()
146 {
147   //UserCreateOutputObject
148   if (fOutputList) fOutputList->Delete();
149   delete fOutputList;
150   fOutputList=new TList();
151   fOutputList->SetOwner(kTRUE);
152
153   const  Int_t ndec=2;
154   Int_t startvalue=-1;
155   const  Int_t npredec=50;
156   Double_t tabx[ndec*npredec+1];
157   for (Int_t i=0; i<ndec; i++)
158   {
159     for (Int_t j=0; j<npredec; j++)
160     {
161       tabx[npredec*i+j]=TMath::Power(10,((Double_t)i)+((Double_t)startvalue)+((Double_t)j)/((Double_t)npredec));
162     }
163   }
164   tabx[ndec*npredec]=TMath::Power(10,ndec+startvalue);
165
166   fITSsignal=new TH2F("fITSsignal",";p [GeV/c];dEdx",ndec*npredec,tabx,900,0,900);
167   fOutputList->Add(fITSsignal);
168   fTPCsignal=new TH2F("fTPCsignal",";p [GeV/c];dEdx",ndec*npredec,tabx,500,0,500);
169   fOutputList->Add(fTPCsignal);
170   fTOFsignal=new TH2F("fTOFsignal",";p [GeV/c];t-t_{#pi}",ndec*npredec,tabx,1000,-2000,10000);
171   fOutputList->Add(fTOFsignal);
172
173   Int_t kPtBins=60;
174   Double_t binsPtDummy[kPtBins+1];
175   binsPtDummy[0]=0.0;
176   for(int i=1; i<=kPtBins+1; i++)
177   {
178     if(binsPtDummy[i-1]+0.05<1.01)
179       binsPtDummy[i]=binsPtDummy[i-1]+0.05;
180     else
181       binsPtDummy[i]=binsPtDummy[i-1]+0.1;
182   }
183
184   Int_t kPBins=60;
185   Double_t binsPDummy[kPBins+1];
186   binsPDummy[0]=0.0;
187   for(int i=1; i<=kPBins+1; i++)
188   {
189     if(binsPDummy[i-1]+0.05<1.01)
190       binsPDummy[i]=binsPDummy[i-1]+0.05;
191     else
192       binsPDummy[i]=binsPDummy[i-1]+0.1;
193   }
194
195   fITSsignalpip=new TH2F("fITSsignalPions",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of p for pi+
196   fTPCsignalpip=new TH2F("fTPCsignalPions",";p [GeV/c];signal",kPBins,binsPDummy,300,-1,1);//TPC PID signal as function of p for pi+
197   fTOFsignalpip=new TH2F("fTOFsignalPions",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of p for pi+
198   fOutputList->Add(fITSsignalpip);
199   fOutputList->Add(fTPCsignalpip);
200   fOutputList->Add(fTOFsignalpip);
201
202   fITSsignalKp=new TH2F("fITSsignalKaons",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of p for K+
203   fTPCsignalKp=new TH2F("fTPCsignalKaons",";p [GeV/c];signal",kPBins,binsPDummy,300,-1,1);//TPC PID signal as function of p for K+
204   fTOFsignalKp=new TH2F("fTOFsignalKaons",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of p for K+
205   fOutputList->Add(fITSsignalKp);
206   fOutputList->Add(fTPCsignalKp);
207   fOutputList->Add(fTOFsignalKp);
208
209   fITSsignalpp=new TH2F("fITSsignalProtons",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of p for p
210   fTPCsignalpp=new TH2F("fTPCsignalProtons",";p [GeV/c];signal",kPBins,binsPDummy,300,-1,1);//TPC PID signal as function of p for p
211   fTOFsignalpp=new TH2F("fTOFsignalProtons",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of p for p
212   fOutputList->Add(fITSsignalpp);
213   fOutputList->Add(fTPCsignalpp);
214   fOutputList->Add(fTOFsignalpp);
215
216   fTOFsignalBeta=new TH2F("fTOFbeta",";p[GeV/c];#beta",kPBins,binsPDummy,1000, 0.4, 1.1);//
217   fTOFsignalPiBeta=new TH2F("fTOFbetaPions",";p [GeV/c];#beta-#beta_{#pi}",kPBins,binsPDummy,500, -0.25, 0.25);//
218   fTOFsignalKBeta=new TH2F("fTOFbetaKaons",";p [GeV/c];#beta-#beta_{K}",kPBins,binsPDummy,500, -0.25, 0.2);//
219   fTOFsignalPBeta=new TH2F("fTOFbetaProtons",";p [GeV/c];#beta-#beta_{p}",kPBins,binsPDummy,500, -0.25, 0.25);//
220   fOutputList->Add(fTOFsignalBeta);
221   fOutputList->Add(fTOFsignalPiBeta);
222   fOutputList->Add(fTOFsignalKBeta);
223   fOutputList->Add(fTOFsignalPBeta);
224
225   if(fMC)
226   {
227     fITSsignalpiMCp=new TH2F("fITSsignalPionsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of pt for pi+
228     fTPCsignalpiMCp=new TH2F("fTPCsignalPionsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-1,1);//TPC PID signal as function of pt for pi+
229     fTOFsignalpiMCp=new TH2F("fTOFsignalPionsMC",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of pt for pi+
230     fOutputList->Add(fITSsignalpiMCp);
231     fOutputList->Add(fTPCsignalpiMCp);
232     fOutputList->Add(fTOFsignalpiMCp);
233
234     fITSsignalKMCp=new TH2F("fITSsignalKaonsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of pt for K+
235     fTPCsignalKMCp=new TH2F("fTPCsignalKaonsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-1,1);//TPC PID signal as function of pt for K+
236     fTOFsignalKMCp=new TH2F("fTOFsignalKaonsMC",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of pt for K+
237     fOutputList->Add(fITSsignalKMCp);
238     fOutputList->Add(fTPCsignalKMCp);
239     fOutputList->Add(fTOFsignalKMCp);
240
241     fITSsignalpMCp=new TH2F("fITSsignalProtonsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-4,4);//ITS PID signal as function of pt for p
242     fTPCsignalpMCp=new TH2F("fTPCsignalProtonsMC",";p [GeV/c];signal",kPBins,binsPDummy,600,-1,1);//TPC PID signal as function of pt for p
243     fTOFsignalpMCp=new TH2F("fTOFsignalProtonsMC",";p [GeV/c];signal",kPBins,binsPDummy,1000,-8000,8000);//TOF PID signal as function of pt for p
244     fOutputList->Add(fITSsignalpMCp);
245     fOutputList->Add(fTPCsignalpMCp);
246     fOutputList->Add(fTOFsignalpMCp);
247   }
248
249   fPvsPt=new TH2F("fPvsPt","p vs p_{t};p [GeV/c];p_{t} [GeV/c]",kPBins,binsPDummy,kPtBins,binsPtDummy);
250   fOutputList->Add(fPvsPt);
251
252   fMeanPvsP = new TProfile("fMeanPvsP","Mean P vs P;p [Gev/c];<p> [GeV/c]",kPBins,binsPDummy);
253   fOutputList->Add(fMeanPvsP);
254
255   fTPCvsGlobalMult = new TH2F("fTPCvsGlobalMult","TPC only vs Global track multiplicity;global;TPC only",500,0,2500,500,0,3500);
256   fOutputList->Add(fTPCvsGlobalMult);
257
258   fStandardGlobalCuts = AliFlowTrackCuts::GetStandardGlobalTrackCuts2010();
259   fStandardTPCCuts = AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts2010();
260
261   //fOutputList->Add(fESDpid);
262
263   PostData(1,  fOutputList);
264 }
265
266 //________________________________________________________________________
267 void  AliAnalysisTaskPIDflowQA::UserExec(Option_t *)
268 {
269   fESD = dynamic_cast<AliESDEvent*> (InputEvent());
270   if (!fESD) return;
271
272   //do the calibration bit
273   fESDpid->SetTOFResponse(fESD,AliESDpid::kTOF_T0); // to use T0-TOF 
274   fESDpid->MakePID(fESD,kFALSE);
275
276   if(!fCuts || !fEventCuts)
277   {
278     Printf("No CUTS Defined.........\n");
279     PostData(1,  fOutputList);
280     return;
281   }
282
283   if (!(fEventCuts->IsSelected(fESD)))
284   {
285     return;
286   }
287
288   AliStack* stack=0x0;
289   if(fMC)
290   {
291     AliMCEvent* mcEvent  = (AliMCEvent*) MCEvent();
292     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
293     stack = mcEvent->Stack();
294   }
295
296   Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
297   Int_t nTracks=fESD->GetNumberOfTracks();
298
299   AliESDtrack *trackESD=0;
300
301   for(int tr1=0; tr1<nTracks; tr1++)
302   {
303     trackESD=fESD->GetTrack(tr1);
304     if (!trackESD) continue;
305
306     if(!(fCuts->AcceptTrack(trackESD))) continue;
307
308     Int_t label=-1;
309     if(fMC) label=trackESD->GetLabel();
310
311     Int_t pdgcode=0;
312     if(stack&&fMC)
313     {
314       TParticle* particle2 = stack->Particle(TMath::Abs(label));
315       pdgcode=particle2->GetPdgCode();
316     }
317
318     Double_t p=trackESD->GetP();
319     Double_t pt=trackESD->Pt();
320     fPvsPt->Fill(p,pt);
321     fMeanPvsP->Fill(p,p);
322
323     pidITS(trackESD,pdgcode);
324     pidTPC(trackESD,pdgcode);
325     pidTOF(trackESD,pdgcode);
326   }
327
328   //check the correlation between the global and TPConly number of tracks
329   fStandardGlobalCuts->SetEvent(fESD);
330   fStandardTPCCuts->SetEvent(fESD);
331   fTPCvsGlobalMult->Fill(fStandardGlobalCuts->Count(),fStandardTPCCuts->Count());
332
333   // Post output data.
334   PostData(1,  fOutputList);
335 }
336
337 //________________________________________________________________________
338 void  AliAnalysisTaskPIDflowQA::Terminate(Option_t *)
339 {
340   //Terminate
341   if(fCuts)
342     fCuts->Dump();
343   if(fMC)
344     Printf("MC On\n");
345
346   Printf("AliAnalysisTaskPIDflowQA: end of Terminate");
347 }
348
349 //________________________________________________________________________
350 void AliAnalysisTaskPIDflowQA::pidITS(AliESDtrack* t, Int_t pdgcode)
351 {
352   Int_t ngoodSDDSSD=0;
353   Double_t sample[4]= {0.0,0.0,0.0,0.0};
354   t->GetITSdEdxSamples(sample);
355   for(int i=0; i<4; i++)
356   {
357     if(sample[i]>50.0)
358       ngoodSDDSSD++;
359   }
360   if(ngoodSDDSSD<3) return;
361
362   Float_t dedx=(Float_t)t->GetITSsignal();
363   if(dedx<1.0) return;
364
365   Bool_t ifSA=!(t->IsOn(AliESDtrack::kTPCin));
366   Float_t p=t->GetP();
367
368   Float_t signalpi=fESDpid->GetITSResponse().Bethe(p,AliPID::ParticleMass(2),ifSA);
369   Float_t signalK=fESDpid->GetITSResponse().Bethe(p,AliPID::ParticleMass(3),ifSA);
370   Float_t signalp=fESDpid->GetITSResponse().Bethe(p,AliPID::ParticleMass(4),ifSA);
371   if(signalpi<1.0||signalK<1.0||signalp<1.0) return;
372
373   fITSsignal->Fill(p,dedx);
374
375   fITSsignalpip->Fill(p,TMath::Log(dedx)-TMath::Log(signalpi));
376   fITSsignalKp->Fill(p,TMath::Log(dedx)-TMath::Log(signalK));
377   fITSsignalpp->Fill(p,TMath::Log(dedx)-TMath::Log(signalp));
378
379   if(fMC)
380   {
381     if(TMath::Abs(pdgcode)==211)
382       fITSsignalpiMCp->Fill(p,TMath::Log(dedx)-TMath::Log(signalpi));
383     else if(TMath::Abs(pdgcode)==321)
384       fITSsignalKMCp->Fill(p,TMath::Log(dedx)-TMath::Log(signalK));
385     else if (TMath::Abs(pdgcode)==2212)
386       fITSsignalpMCp->Fill(p,TMath::Log(dedx)-TMath::Log(signalp));
387   }
388 }
389
390 //________________________________________________________________________
391 void AliAnalysisTaskPIDflowQA::pidTPC(AliESDtrack* t, Int_t pdgcode)
392 {
393   Double_t pinTPCglobal=t->GetInnerParam()->GetP();
394   Float_t sigPion     = fESDpid->GetTPCResponse().GetExpectedSignal(pinTPCglobal, AliPID::kPion);
395   Float_t sigKaon     = fESDpid->GetTPCResponse().GetExpectedSignal(pinTPCglobal, AliPID::kKaon);
396   Float_t sigProton   = fESDpid->GetTPCResponse().GetExpectedSignal(pinTPCglobal, AliPID::kProton);
397   Float_t p=t->GetP();
398   Double_t tpcSignal =t ->GetTPCsignal();
399   if(!(sigPion>0.0&&sigKaon>0.0&&sigProton>0.0))
400     return;
401
402   fTPCsignal->Fill(pinTPCglobal,tpcSignal);
403
404   fTPCsignalpip->Fill(p,(tpcSignal-sigPion)/sigPion);
405   fTPCsignalKp->Fill(p,(tpcSignal-sigKaon)/sigKaon);
406   fTPCsignalpp->Fill(p,(tpcSignal-sigProton)/sigProton);
407
408   if(fMC)
409   {
410     if(TMath::Abs(pdgcode)==211)
411       fTPCsignalpiMCp->Fill(p,(tpcSignal-sigPion)/sigPion);
412     else if(TMath::Abs(pdgcode)==321)
413       fTPCsignalKMCp->Fill(p,(tpcSignal-sigKaon)/sigKaon);
414     else if (TMath::Abs(pdgcode)==2212)
415       fTPCsignalpMCp->Fill(p,(tpcSignal-sigProton)/sigProton);
416   }
417 }
418
419 //________________________________________________________________________
420 void AliAnalysisTaskPIDflowQA::pidTOF(AliESDtrack* t, Int_t pdgcode)
421 {
422   Bool_t goodtrack = (t) &&
423                      (t->GetStatus() & AliESDtrack::kTOFpid) &&
424                      (t->GetTOFsignal() > 12000) &&
425                      (t->GetTOFsignal() < 100000) &&
426                      (t->GetIntegratedLength() > 365) &&
427                      !(t->GetStatus() & AliESDtrack::kTOFmismatch);
428
429   if (!goodtrack) return;
430
431   const Float_t c = 2.99792457999999984e-02;
432   Float_t p=t->GetP();
433   Float_t L = t->GetIntegratedLength();
434   Float_t fT0track=fESDpid->GetTOFResponse().GetStartTime(p);
435   Float_t timeTOF=t->GetTOFsignal()- fT0track;
436
437   //calculate beta for the track
438   Float_t beta = L/timeTOF/c;
439   
440   //2=pion 3=kaon 4=protons
441   Double_t inttimes[5]= {-1.0,-1.0,-1.0,-1.0,-1.0};
442   t->GetIntegratedTimes(inttimes);
443   Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
444   for (Int_t i=0;i<5;i++)
445   {
446     betaHypothesis[i] = L/inttimes[i]/c;
447   }
448
449   fTOFsignal->Fill(p,timeTOF-inttimes[2]);
450
451   //beta part
452   fTOFsignalBeta->Fill(p,beta);
453   fTOFsignalPiBeta->Fill(p,beta-betaHypothesis[2]);
454   fTOFsignalKBeta->Fill(p,beta-betaHypothesis[3]);
455   fTOFsignalPBeta->Fill(p,beta-betaHypothesis[4]);
456
457   //P part
458   fTOFsignalpip->Fill(p,timeTOF-inttimes[2]);
459   fTOFsignalKp->Fill(p,timeTOF-inttimes[3]);
460   fTOFsignalpp->Fill(p,timeTOF-inttimes[4]);
461
462   if(fMC)
463   {
464     if(TMath::Abs(pdgcode)==211)
465       fTOFsignalpiMCp->Fill(p,timeTOF-inttimes[2]);
466     else if(TMath::Abs(pdgcode)==321)
467       fTOFsignalKMCp->Fill(p,timeTOF-inttimes[3]);
468     else if (TMath::Abs(pdgcode)==2212)
469       fTOFsignalpMCp->Fill(p,timeTOF-inttimes[4]);
470   }
471 }
472
473 Float_t AliAnalysisTaskPIDflowQA::Beta(Float_t m, Float_t p) 
474 {
475   //get theoretical beta
476   return TMath::Sqrt(1. / (1. + m * m / (p * p)));
477 }
478