1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //Task for analysis if TPC dEdx and PID information
20 //-----------------------------------------------------------------------
21 // Author : M.Ivanov marian.ivanov@cern.ch -
22 //-----------------------------------------------------------------------
25 // 3 6D histograms - THnSparse created in the task:
27 // TPC normalized dEdx (dEdx_rec/dNdx_mc)
28 // TPC PID probabilities
30 // The values are binned in following variables:
31 // Some of them are correlated - but THnSpase handle it
32 // ~ 14 MBy per object needed
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]
39 // 5 - measurement - dEdx, dEdx/BB resp. PID probability
40 // 6 - BB resp. rec particle
50 #include <TParticle.h>
53 #include <AliAnalysisManager.h>
54 #include <AliESDInputHandler.h>
56 #include "AliMCEvent.h"
57 #include "AliMCEventHandler.h"
58 #include "AliMathBase.h"
61 #include "AliExternalTrackParam.h"
62 #include "AliTracker.h"
63 #include "AliTPCseed.h"
65 #include "AliTPCtaskPID.h"
67 #include <THnSparse.h>
77 ClassImp(AliTPCtaskPID)
79 //________________________________________________________________________
80 AliTPCtaskPID::AliTPCtaskPID() :
82 fMCinfo(0), //! MC event handler
90 // Default constructor (should not be used)
94 AliTPCtaskPID::AliTPCtaskPID(const AliTPCtaskPID& info) :
95 AliAnalysisTask(info),
96 fMCinfo(info.fMCinfo), //! MC event handler
104 // Dummy Copy constructor - no copy constructor for THnSparse
106 fList = (TObjArray*)(info.fList->Clone());
111 //________________________________________________________________________
112 AliTPCtaskPID::AliTPCtaskPID(const char *name) :
113 AliAnalysisTask(name, "AliTPCtaskPID"),
114 fMCinfo(0), //! MC event handler
122 // Normal constructor
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());
133 void AliTPCtaskPID::Init(){
135 // Init dEdx histogram
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]
142 // 5 - measurement - dEdx or dEdx/BB
145 Double_t xmin[7], xmax[7];
149 xmin[0]=0; xmax[0]=10;
152 xmin[1]=0.1; xmax[1]=100;
155 xmin[2]=-1.4; xmax[2]=1.4;
160 xmin[3]=0.1; xmax[3]=1000;
163 xmin[4] =50; xmax[4]=160;
167 xmin[6]=0.5; xmax[6]=6;
170 xmin[5]=20; xmax[5]=400;
171 fTPCsignal = new THnSparseF("TPC signal","TPC signal",7,nbins,xmin,xmax);
173 xmin[5]=25; xmax[5]=75;
174 fTPCsignalNorm = new THnSparseF("TPC signal Norm","TPC signal Norm",7,nbins,xmin,xmax);
177 xmin[5]=-0.1; xmax[5]=1.1;
179 xmin[6]=0; xmax[6]=10;
180 fTPCr = new THnSparseF("TPC pid probability ","TPC pid probability",7,nbins,xmin,xmax);
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));
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"};
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]);
205 fList = new TObjArray(3);
206 fList->AddAt(fTPCsignal,0);
207 fList->AddAt(fTPCsignalNorm,1);
208 fList->AddAt(fTPCr,2);
214 AliTPCtaskPID::~AliTPCtaskPID(){
219 delete fTPCsignalNorm;
224 //________________________________________________________________________
225 void AliTPCtaskPID::ConnectInputData(Option_t *)
228 // Connect the input data
231 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
233 //Printf("ERROR: Could not read chain from input slot 0");
236 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
238 //Printf("ERROR: Could not get ESDInputHandler");
241 fESD = esdH->GetEvent();
242 //Printf("*** CONNECTED NEW EVENT ****");
245 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
246 mcinfo->SetReadTR(kTRUE);
248 fMCinfo = mcinfo->MCEvent();
258 //________________________________________________________________________
259 void AliTPCtaskPID::CreateOutputObjects()
262 // Connect the output objects
268 //________________________________________________________________________
269 void AliTPCtaskPID::Exec(Option_t *) {
271 // Execute analysis for current event
275 // If MC has been connected
278 cout << "Not MC info\n" << endl;
291 //________________________________________________________________________
292 void AliTPCtaskPID::Terminate(Option_t *) {
306 void AliTPCtaskPID::ProcessMCInfo(){
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;
319 TParticle * particle= new TParticle;
320 TClonesArray * trefs = new TClonesArray("AliTrackReference");
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();
327 Int_t ipart = TMath::Abs(track->GetLabel());
329 Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
330 if (status<0) continue;
331 if (!particle) continue;
332 if (!trefs) continue;
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;
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);
348 Int_t pdg = particle->GetPdgCode();
349 for (Int_t iType=0;iType<5;iType++) {
350 if (AliPID::ParticleCode(iType)==TMath::Abs(pdg)){
352 if (TDatabasePDG::Instance()->GetParticle(pdg)->Charge()<0) x[0]+=5;
356 x[2]= -0.5*TMath::Log((track->P()+track->Pz())/(track->P()-track->Pz()));
358 x[4]=track->GetTPCsignalN();
364 x[5]=dedx/betheBloch;
365 fTPCsignalNorm->Fill(x);
366 for (Int_t ipart2=0;ipart2<5;ipart2++){
368 Double_t tpcpid[AliPID::kSPECIES];
369 track->GetTPCpid(tpcpid);
379 void AliTPCtaskPID::FinishTaskOutput()
382 // According description in AliAnalisysTask this method is call
383 // on the slaves before sending data
386 gSystem->Exec("pwd");
392 void AliTPCtaskPID::BinLogX(TAxis *axis) {
396 Int_t bins = axis->GetNbins();
398 Double_t from = axis->GetXmin();
399 Double_t to = axis->GetXmax();
400 Double_t *new_bins = new Double_t[bins + 1];
403 Double_t factor = pow(to/from, 1./bins);
405 for (Int_t i = 1; i <= bins; i++) {
406 new_bins[i] = factor * new_bins[i-1];
408 axis->Set(bins, new_bins);