don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliTPCtaskPID.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 //
17 // 
18 //Task for analysis if TPC dEdx and PID information 
19 //
20 //-----------------------------------------------------------------------
21 // Author : M.Ivanov  marian.ivanov@cern.ch - 
22 //-----------------------------------------------------------------------
23
24
25 // 3 6D histograms  - THnSparse created in the task:
26 // TPC raw dEdx
27 // TPC normalized dEdx (dEdx_rec/dNdx_mc)
28 // TPC PID probabilities
29 //
30 // The values are binned in following variables:
31 // Some of them are correlated - but THnSpase handle it  
32 //                               ~ 14 MBy per object needed
33 //
34 // 0 - MC particle species as defined in the AliPID - negatives+5 <0,9>
35 // 1 - momenta - at the entrance of the TPC
36 // 2 - tan lambda- fP[3]
37 // 3 - betagamma
38 // 4 - npoints
39 // 5 - measurement - dEdx, dEdx/BB resp.  PID probability
40 // 6 - BB resp. rec particle
41
42
43 //
44 // ROOT includes
45 #include <TChain.h>
46 #include <TMath.h>
47 #include <TVectorD.h>
48 #include <TSystem.h>
49 #include <TFile.h>
50 #include <TParticle.h>
51
52 // ALIROOT includes
53 #include <AliAnalysisManager.h>
54 #include <AliESDInputHandler.h>
55 #include "AliStack.h"
56 #include "AliMCEvent.h"
57 #include "AliMCEventHandler.h"
58 #include "AliMathBase.h"
59
60 #include <AliESD.h>
61 #include "AliExternalTrackParam.h"
62 #include "AliTracker.h"
63 #include "AliTPCseed.h"
64 //
65 #include "AliTPCtaskPID.h"
66 //
67 #include <THnSparse.h>
68
69 //
70
71 // STL includes
72 #include <iostream>
73
74 using namespace std;
75
76
77 ClassImp(AliTPCtaskPID)
78
79 //________________________________________________________________________
80 AliTPCtaskPID::AliTPCtaskPID() : 
81   AliAnalysisTask(), 
82   fMCinfo(0),     //! MC event handler
83   fESD(0),
84   fList(0),
85   fTPCsignal(0),
86   fTPCsignalNorm(0),
87   fTPCr(0)
88 {
89   //
90   // Default constructor (should not be used)
91   //
92 }
93
94 AliTPCtaskPID::AliTPCtaskPID(const AliTPCtaskPID& info) : 
95   AliAnalysisTask(info), 
96   fMCinfo(info.fMCinfo),     //! MC event handler
97   fESD(info.fESD),        //!
98   fList(0),
99   fTPCsignal(0),
100   fTPCsignalNorm(0),
101   fTPCr(0)
102 {
103   //
104   // Dummy Copy  constructor - no copy constructor for THnSparse 
105   //
106   fList = (TObjArray*)(info.fList->Clone());
107 }
108
109
110
111 //________________________________________________________________________
112 AliTPCtaskPID::AliTPCtaskPID(const char *name) : 
113   AliAnalysisTask(name, "AliTPCtaskPID"), 
114   fMCinfo(0),     //! MC event handler
115   fESD(0),
116   fList(0),
117   fTPCsignal(0),
118   fTPCsignalNorm(0),
119   fTPCr(0)
120 {
121   //
122   // Normal constructor
123   //
124   // Input slot #0 works with a TChain
125   DefineInput(0, TChain::Class());
126   // Output slot #0 writes into a TList
127   DefineOutput(0, TObjArray::Class());
128   //
129   //make histos
130   Init(); 
131 }
132
133 void AliTPCtaskPID::Init(){
134   //
135   // Init dEdx histogram
136   // Dimensions
137   // 0 - particle specie as defined in the AliPID - negatives+5 <0,9>
138   // 1 - momenta - at the entrance of the TPC
139   // 2 - tan lambda- fP[3]
140   // 3 - betagamma
141   // 4 - npoints
142   // 5 - measurement - dEdx or dEdx/BB
143   // 6 - BB
144   //
145   Double_t xmin[7],  xmax[7];
146   Int_t    nbins[7];
147   // pid
148   nbins[0]=10;
149   xmin[0]=0; xmax[0]=10;
150   // momenta
151   nbins[1]=50;
152   xmin[1]=0.1; xmax[1]=100;
153   //pseudorapidity
154   nbins[2]=40;
155   xmin[2]=-1.4; xmax[2]=1.4;
156   //
157   // betagamma
158   //
159   nbins[3]=100;
160   xmin[3]=0.1; xmax[3]=1000;
161   //
162   nbins[4]=11;
163   xmin[4] =50; xmax[4]=160;
164   //
165   // 
166   nbins[6]=60;
167   xmin[6]=0.5; xmax[6]=6;
168
169   nbins[5]=400;
170   xmin[5]=20; xmax[5]=400;
171   fTPCsignal = new THnSparseF("TPC signal","TPC signal",7,nbins,xmin,xmax);
172   nbins[5]=100;
173   xmin[5]=25; xmax[5]=75;
174   fTPCsignalNorm = new THnSparseF("TPC signal Norm","TPC signal Norm",7,nbins,xmin,xmax);
175   //
176   nbins[5]=256;
177   xmin[5]=-0.1; xmax[5]=1.1;
178   nbins[6]=10;
179   xmin[6]=0; xmax[6]=10;
180   fTPCr = new THnSparseF("TPC pid probability ","TPC pid probability",7,nbins,xmin,xmax);
181   //
182   //
183   BinLogX(fTPCsignal->GetAxis(1));
184   BinLogX(fTPCsignal->GetAxis(3));
185   BinLogX(fTPCsignal->GetAxis(6));
186   BinLogX(fTPCsignalNorm->GetAxis(1));
187   BinLogX(fTPCsignalNorm->GetAxis(3));
188   BinLogX(fTPCsignalNorm->GetAxis(6));
189   BinLogX(fTPCr->GetAxis(1));
190   BinLogX(fTPCr->GetAxis(3));
191   
192   const char *hisAxisName[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","dEdx_{rec},dedx_{mc}"};
193   const char *hisAxisNameNorm[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","dEdx_{rec}/dEdx_{mc},dedx_{mc}"};
194   const char *hisAxisNameR[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","TPCr","pid2"};
195   //  
196   for (Int_t i=0;i<7;i++) {
197     fTPCsignal->GetAxis(i)->SetTitle(hisAxisName[i]);
198     fTPCsignal->GetAxis(i)->SetName(hisAxisName[i]);
199     fTPCsignalNorm->GetAxis(i)->SetTitle(hisAxisNameNorm[i]);
200     fTPCsignalNorm->GetAxis(i)->SetName(hisAxisNameNorm[i]);
201     fTPCr->GetAxis(i)->SetTitle(hisAxisNameR[i]);
202     fTPCr->GetAxis(i)->SetName(hisAxisNameR[i]);
203   }
204
205   fList = new TObjArray(3);
206   fList->AddAt(fTPCsignal,0);
207   fList->AddAt(fTPCsignalNorm,1);
208   fList->AddAt(fTPCr,2);
209 }
210
211
212
213
214 AliTPCtaskPID::~AliTPCtaskPID(){
215   //
216   //
217   //
218   delete fTPCsignal;
219   delete fTPCsignalNorm;
220   delete fTPCr;
221 }
222
223
224 //________________________________________________________________________
225 void AliTPCtaskPID::ConnectInputData(Option_t *) 
226 {
227   //
228   // Connect the input data
229   //
230
231   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
232   if (!tree) {
233     //Printf("ERROR: Could not read chain from input slot 0");
234   }
235   else {
236     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
237     if (!esdH) {
238       //Printf("ERROR: Could not get ESDInputHandler");
239     }
240     else {
241       fESD = esdH->GetEvent();
242       //Printf("*** CONNECTED NEW EVENT ****");
243     }  
244   }
245   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  
246   mcinfo->SetReadTR(kTRUE);
247   
248   fMCinfo = mcinfo->MCEvent();
249
250
251 }
252
253
254
255
256
257
258 //________________________________________________________________________
259 void AliTPCtaskPID::CreateOutputObjects() 
260 {
261   //
262   // Connect the output objects
263   //
264
265 }
266
267
268 //________________________________________________________________________
269 void AliTPCtaskPID::Exec(Option_t *) {
270   //
271   // Execute analysis for current event 
272   //
273
274
275   // If MC has been connected   
276
277   if (!fMCinfo){
278     cout << "Not MC info\n" << endl;
279   }else{
280     ProcessMCInfo();
281     //mcinfo->Print();
282     //DumpInfo();
283   }
284   //
285   PostData(0, fList);
286 }      
287
288
289
290
291 //________________________________________________________________________
292 void AliTPCtaskPID::Terminate(Option_t *) {
293     //
294     // Terminate loop
295     //
296   //
297   return;
298 }
299
300
301
302
303
304
305
306 void  AliTPCtaskPID::ProcessMCInfo(){
307   //
308   //
309   //
310   //
311
312   if (!fTPCsignal) Init();
313   Int_t npart   = fMCinfo->GetNumberOfTracks();
314   Int_t ntracks = fESD->GetNumberOfTracks(); 
315   if (npart<=0) return;
316   if (ntracks<=0) return;
317   //
318   //
319   TParticle * particle= new TParticle;
320   TClonesArray * trefs = new TClonesArray("AliTrackReference");
321   
322   for (Int_t itrack=0;itrack<ntracks;itrack++){
323     AliESDtrack *track = fESD->GetTrack(itrack);
324     const AliExternalTrackParam *in=track->GetInnerParam();
325     const AliExternalTrackParam *out=track->GetOuterParam();
326     if (!in) continue;
327     Int_t ipart = TMath::Abs(track->GetLabel());
328     //
329     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
330     if (status<0) continue;
331     if (!particle) continue;
332     if (!trefs) continue;
333     //
334     //
335     Double_t mom = in->GetP();
336     if (track->GetP()>5)  mom= track->GetP();
337     if (out&&out->GetX()<260&&TMath::Abs(out->GetZ())<250)  mom= (in->GetP()+out->GetP())*0.5;
338     //
339     Double_t dedx=track->GetTPCsignal();
340     Double_t mass = particle->GetMass();
341     Double_t bg  =mom/mass;
342     Double_t betheBloch = AliMathBase::BetheBlochAleph(bg);
343     //
344     // Fill histos
345     //
346     Double_t x[7];
347     //PID
348     Int_t pdg = particle->GetPdgCode();
349     for (Int_t iType=0;iType<5;iType++) {
350       if (AliPID::ParticleCode(iType)==TMath::Abs(pdg)){
351         x[0]=iType;
352         if (TDatabasePDG::Instance()->GetParticle(pdg)->Charge()<0) x[0]+=5;
353       }
354     }
355     x[1]= mom;
356     x[2]= -0.5*TMath::Log((track->P()+track->Pz())/(track->P()-track->Pz()));
357     x[3]= bg;
358     x[4]=track->GetTPCsignalN();
359     x[5]= dedx;
360     x[6]= betheBloch;
361     //
362     fTPCsignal->Fill(x);
363     //
364     x[5]=dedx/betheBloch;
365     fTPCsignalNorm->Fill(x);
366     for (Int_t ipart2=0;ipart2<5;ipart2++){
367       x[6]=ipart2;
368       Double_t tpcpid[AliPID::kSPECIES];
369       track->GetTPCpid(tpcpid);
370       x[5]=tpcpid[ipart2];
371       fTPCr->Fill(x);
372     }
373   }
374 }
375
376
377
378
379 void AliTPCtaskPID::FinishTaskOutput()
380 {
381   //
382   // According description in AliAnalisysTask this method is call
383   // on the slaves before sending data
384   //
385   Terminate("slave");
386   gSystem->Exec("pwd");
387 }
388
389
390
391
392 void AliTPCtaskPID::BinLogX(TAxis *axis) {
393   //
394   //
395   //
396   Int_t bins = axis->GetNbins();
397   
398   Double_t from = axis->GetXmin();
399   Double_t to   = axis->GetXmax();
400   Double_t *new_bins = new Double_t[bins + 1];
401   
402   new_bins[0] = from;
403   Double_t factor = pow(to/from, 1./bins);
404   
405   for (Int_t i = 1; i <= bins; i++) {
406     new_bins[i] = factor * new_bins[i-1];
407   }
408   axis->Set(bins, new_bins);
409   delete [] new_bins;
410 }