]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx
introducing standard contructor and some access functions for data members
[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAnalysisTaskITSsaTracks.cxx
CommitLineData
045571a6 1#include "AliAnalysisTaskSE.h"
2#include "AliTrackPointArray.h"
3#include "AliAnalysisManager.h"
4#include "AliAnalysisDataContainer.h"
5#include "AliESDEvent.h"
6#include "AliESDfriend.h"
7#include "AliStack.h"
8#include "AliPID.h"
9#include "AliITSPIDResponse.h"
10#include "AliMCEventHandler.h"
11#include "AliMCEvent.h"
12#include <TParticle.h>
13#include <TSystem.h>
14#include <TTree.h>
2afdea69 15#include <TNtuple.h>
045571a6 16#include <TH1F.h>
17#include <TH2F.h>
18#include <TChain.h>
19#include "AliESDInputHandlerRP.h"
20#include "AliAnalysisTaskITSsaTracks.h"
21
22/**************************************************************************
23 * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
24 * *
25 * Author: The ALICE Off-line Project. *
26 * Contributors are mentioned in the code where appropriate. *
27 * *
28 * Permission to use, copy, modify and distribute this software and its *
29 * documentation strictly for non-commercial purposes is hereby granted *
30 * without fee, provided that the above copyright notice appears in all *
31 * copies and that both the copyright notice and this permission notice *
32 * appear in the supporting documentation. The authors make no claims *
33 * about the suitability of this software for any purpose. It is *
34 * provided "as is" without express or implied warranty. *
35 **************************************************************************/
36
37//*************************************************************************
38// Implementation of class AliAnalysisTaskITSsaTracks
39// AliAnalysisTaskSE to extract QA and performance histos for ITS standalone tracks
40//
41//
42// Authors: L. Milano, milano@to.infn.it
43// F. Prino, prino@to.infn.it
44//
45//*************************************************************************
46
47ClassImp(AliAnalysisTaskITSsaTracks)
48//______________________________________________________________________________
49AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("ITSsa resolution"),
50 fOutput(0),
51 fHistNEvents(0),
2afdea69 52 fHistMCPhiResid(0),
53 fHistPhiResid(0),
54 fNtupleTracks(0),
045571a6 55 fNPtBins(kMaxPtBins),
56 fMinITSpts(4),
57 fMinSPDpts(1),
2afdea69 58 fMinPtsforPid(3),
045571a6 59 fMinTPCpts(50),
2afdea69 60 fMaxITSChi2Clu(100.),
61 fFillNtuple(kFALSE),
045571a6 62 fReadMC(kFALSE),
63 fUseMCId(kFALSE)
64{
65 //
66 for(Int_t ilay=0; ilay<6;ilay++) fRequirePoint[ilay]=kFALSE;
67 DefineInput(0, TChain::Class());
68 DefineOutput(1, TList::Class());
69 const Int_t nbins = 29;
70 Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
71 0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
72 0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
73 SetPtBins(nbins,xbins);
74}
75
76
77//___________________________________________________________________________
78AliAnalysisTaskITSsaTracks::~AliAnalysisTaskITSsaTracks(){
79 //
2afdea69 80 delete fHistNEvents;
81 for(Int_t iType=0; iType<kNtrackTypes; iType++){
82 delete fHistPt[iType];
83 delete fHistPtGood[iType];
84 delete fHistPtFake[iType];
85 delete fHistEtaPhi[iType];
86 delete fHistEtaPhiGood[iType];
87 delete fHistEtaPhiFake[iType];
88 delete fHistChi2[iType];
89 delete fHistChi2Good[iType];
90 delete fHistChi2Fake[iType];
91 delete fHistNclu[iType];
92 delete fHistNcluGood[iType];
93 delete fHistNcluFake[iType];
94 delete fHistdedxvsPt3cls[iType];
95 delete fHistdedxvsPt4cls[iType];
96 }
97 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
98 delete fHistPtTPCITS[iSpec];
99 delete fHistPtITSsa[iSpec];
100 delete fHistPtITSpureSA[iSpec];
101 delete fHistEtaPhiTPCITS[iSpec];
102 delete fHistEtaPhiITSsa[iSpec];
103 delete fHistEtaPhiITSpureSA[iSpec];
104 delete fHistNcluTPCITS[iSpec];
105 delete fHistNcluITSsa[iSpec];
106 delete fHistNcluITSpureSA[iSpec];
107 delete fHistd0rphiITSpureSA[iSpec];
108 delete fHistd0zITSpureSA[iSpec];
109 delete fHistCluInLayTPCITS[iSpec];
110 delete fHistCluInLayITSsa[iSpec];
111 delete fHistCluInLayITSpureSA[iSpec];
112 delete fHistOuterLayITSpureSA[iSpec];
113 delete fHistPtResid[iSpec];
114 delete fHistPtRelResid[iSpec];
115 delete fHistInvPtResid[iSpec];
116 delete fHistInvPtRelResid[iSpec];
117 delete fHistMCPtResid[iSpec];
118 delete fHistMCPtRelResid[iSpec];
119 delete fHistMCInvPtResid[iSpec];
120 delete fHistMCInvPtRelResid[iSpec];
121 }
122
123 delete fHistMCPhiResid;
124 delete fHistPhiResid;
125 delete fNtupleTracks;
045571a6 126 if (fOutput) {
127 delete fOutput;
128 fOutput = 0;
129 }
130}
131
132//________________________________________________________________________
133void AliAnalysisTaskITSsaTracks::SetPtBins(Int_t n, Double_t* lim){
134 // define pt bins for analysis
135 if(n>kMaxPtBins){
136 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
137 return;
138 }else{
139 fNPtBins=n;
140 for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
141 for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
142 }
143}
144//___________________________________________________________________________
145void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
146 // create output histos
147
148 fOutput = new TList();
149 fOutput->SetOwner();
150 fOutput->SetName("OutputHistos");
151
2afdea69 152 fNtupleTracks= new TNtuple("ntupleTracks","ntupleTracks","pt:px:py:pz:eta:phi:d0xy:d0z:dedx3:dedx4:dedx5:dedx6:truncmean:chi2:status:ITSrefit:TPCin:TPCrefit:ITSpureSA:nITSclu:nTPCclu:clumap:spdin:spdout:sddin:sddout:ssdin:ssdout:label:ptgen:pxgen:pygen:pzgen:etagen:phigen:ntracks");
153 // kinematics: pt, p eta,phi -> 6 variables
154 // impact parameters: d0xy, d0z -> 2 variables
155 // dE/dx: 4 Layers + trunc mean -> 5 variables
156 // chi2 and track status: -> 6 variables
157 // cluster infos: -> 9 variables
158 // MC info: -> 7 variables
159 // Multiplicity -> 1 variable
160 // Total: 36
161 fOutput->Add(fNtupleTracks);
162
045571a6 163 fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
164 fHistNEvents->Sumw2();
165 fHistNEvents->SetMinimum(0);
166 fOutput->Add(fHistNEvents);
167
168 //binning for the dedx histogram
169 const Int_t hnbinsdedx=400;
170 Double_t hxmindedx = 0.01;
171 Double_t hxmaxdedx = 10;
172 Double_t hlogxmindedx = TMath::Log10(hxmindedx);
173 Double_t hlogxmaxdedx = TMath::Log10(hxmaxdedx);
174 Double_t hbinwidthdedx = (hlogxmaxdedx-hlogxmindedx)/hnbinsdedx;
175 Double_t hxbinsdedx[hnbinsdedx+1];
176 hxbinsdedx[0] = 0.01;
177 for (Int_t i=1;i<=hnbinsdedx;i++) {
178 hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
179 }
180
2afdea69 181 TString tit[kNtrackTypes]={"TPCITS","ITSsa","ITSpureSA"};
182 for(Int_t iType=0; iType<kNtrackTypes; iType++){
183 fHistPt[iType]=new TH1F(Form("hPt%s",tit[iType].Data()),"",100,0.,2.);
184 fHistPt[iType]->Sumw2();
185 fOutput->Add(fHistPt[iType]);
186 fHistPtGood[iType]=new TH1F(Form("hPtGood%s",tit[iType].Data()),"",100,0.,2.);
187 fHistPtGood[iType]->Sumw2();
188 fOutput->Add(fHistPtGood[iType]);
189 fHistPtFake[iType]=new TH1F(Form("hPtFake%s",tit[iType].Data()),"",100,0.,2.);
190 fHistPtFake[iType]->Sumw2();
191 fOutput->Add(fHistPtFake[iType]);
192
193 fHistEtaPhi[iType] = new TH2F(Form("hEtaPhi%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
194 fHistEtaPhi[iType]->Sumw2();
195 fOutput->Add(fHistEtaPhi[iType]);
196 fHistEtaPhiGood[iType] = new TH2F(Form("hEtaPhiGood%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
197 fHistEtaPhiGood[iType]->Sumw2();
198 fOutput->Add(fHistEtaPhiGood[iType]);
199 fHistEtaPhiFake[iType] = new TH2F(Form("hEtaPhiFake%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
200 fHistEtaPhiFake[iType]->Sumw2();
201 fOutput->Add(fHistEtaPhiFake[iType]);
202
203 fHistChi2[iType]=new TH1F(Form("hChi2%s",tit[iType].Data()),"",100,0.,10.);
204 fHistChi2[iType]->Sumw2();
205 fOutput->Add(fHistChi2[iType]);
206 fHistChi2Good[iType]=new TH1F(Form("hChi2Good%s",tit[iType].Data()),"",100,0.,10.);
207 fHistChi2Good[iType]->Sumw2();
208 fOutput->Add(fHistChi2Good[iType]);
209 fHistChi2Fake[iType]=new TH1F(Form("hChi2Fake%s",tit[iType].Data()),"",100,0.,10.);
210 fHistChi2Fake[iType]->Sumw2();
211 fOutput->Add(fHistChi2Fake[iType]);
212
213 fHistNclu[iType]=new TH1F(Form("hNclu%s",tit[iType].Data()),"",7,-0.5,6.5);
214 fHistNclu[iType]->Sumw2();
215 fOutput->Add(fHistNclu[iType]);
216 fHistNcluGood[iType]=new TH1F(Form("hNcluGood%s",tit[iType].Data()),"",7,-0.5,6.5);
217 fHistNcluGood[iType]->Sumw2();
218 fOutput->Add(fHistNcluGood[iType]);
219 fHistNcluFake[iType]=new TH1F(Form("hNcluFake%s",tit[iType].Data()),"",7,-0.5,6.5);
220 fHistNcluFake[iType]->Sumw2();
221 fOutput->Add(fHistNcluFake[iType]);
222
223 fHistdedxvsPt3cls[iType] = new TH2F(Form("hdedxvsPt3cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
224 fHistdedxvsPt3cls[iType]->Sumw2();
225 fOutput->Add(fHistdedxvsPt3cls[iType]);
226
227 fHistdedxvsPt4cls[iType] = new TH2F(Form("hdedxvsPt4cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
228 fHistdedxvsPt4cls[iType]->Sumw2();
229 fOutput->Add(fHistdedxvsPt4cls[iType]);
230 }
231
232 //-----------------------------------------------------------
233
045571a6 234 TString spname[3]={"Pion","Kaon","Proton"};
235 TString hisname;
236 const Int_t nbins = fNPtBins;
237 Double_t xbins[nbins+1];
238 for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
239
045571a6 240 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
241
242 hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
243 fHistPtTPCITS[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
244 fHistPtTPCITS[iSpec]->Sumw2();
245 fOutput->Add(fHistPtTPCITS[iSpec]);
246
247 hisname.Form("hPtITSsa%s",spname[iSpec].Data());
248 fHistPtITSsa[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
249 fHistPtITSsa[iSpec]->Sumw2();
250 fOutput->Add(fHistPtITSsa[iSpec]);
251
252 hisname.Form("hPtITSpureSA%s",spname[iSpec].Data());
253 fHistPtITSpureSA[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
254 fHistPtITSpureSA[iSpec]->Sumw2();
255 fOutput->Add(fHistPtITSpureSA[iSpec]);
256
257 //---
258
259 hisname.Form("hEtaPhiTPCITS%s",spname[iSpec].Data());
260 fHistEtaPhiTPCITS[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
261 fHistEtaPhiTPCITS[iSpec]->Sumw2();
262 fOutput->Add(fHistEtaPhiTPCITS[iSpec]);
263
264 hisname.Form("hEtaPhiITSsa%s",spname[iSpec].Data());
265 fHistEtaPhiITSsa[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
266 fHistEtaPhiITSsa[iSpec]->Sumw2();
267 fOutput->Add(fHistEtaPhiITSsa[iSpec]);
268
269 hisname.Form("hEtaPhiITSpureSA%s",spname[iSpec].Data());
270 fHistEtaPhiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
271 fHistEtaPhiITSpureSA[iSpec]->Sumw2();
272 fOutput->Add(fHistEtaPhiITSpureSA[iSpec]);
273
274 //---
275
276 hisname.Form("hNcluTPCITS%s",spname[iSpec].Data());
277 fHistNcluTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
278 fHistNcluTPCITS[iSpec]->Sumw2();
279 fOutput->Add(fHistNcluTPCITS[iSpec]);
280
281 hisname.Form("hNcluITSsa%s",spname[iSpec].Data());
282 fHistNcluITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
283 fHistNcluITSsa[iSpec]->Sumw2();
284 fOutput->Add(fHistNcluITSsa[iSpec]);
285
286 hisname.Form("hNcluITSpureSA%s",spname[iSpec].Data());
287 fHistNcluITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
288 fHistNcluITSpureSA[iSpec]->Sumw2();
289 fOutput->Add(fHistNcluITSpureSA[iSpec]);
290
291 //---
292
293 hisname.Form("hd0rphiITSpureSA%s",spname[iSpec].Data());
294 fHistd0rphiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
295 fHistd0rphiITSpureSA[iSpec]->Sumw2();
296 fOutput->Add(fHistd0rphiITSpureSA[iSpec]);
297
298 hisname.Form("hd0zITSpureSA%s",spname[iSpec].Data());
299 fHistd0zITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
300 fHistd0zITSpureSA[iSpec]->Sumw2();
301 fOutput->Add(fHistd0zITSpureSA[iSpec]);
302
303 //---
304
305 hisname.Form("hCluInLayTPCITS%s",spname[iSpec].Data());
306 fHistCluInLayTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
307 fHistCluInLayTPCITS[iSpec]->Sumw2();
308 fOutput->Add(fHistCluInLayTPCITS[iSpec]);
309
310 hisname.Form("hCluInLayITSsa%s",spname[iSpec].Data());
311 fHistCluInLayITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
312 fHistCluInLayITSsa[iSpec]->Sumw2();
313 fOutput->Add(fHistCluInLayITSsa[iSpec]);
314
315 hisname.Form("hCluInLayITSpureSA%s",spname[iSpec].Data());
316 fHistCluInLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
317 fHistCluInLayITSpureSA[iSpec]->Sumw2();
318 fOutput->Add(fHistCluInLayITSpureSA[iSpec]);
319
320 hisname.Form("hOuterLayITSpureSA%s",spname[iSpec].Data());
321 fHistOuterLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
322 fHistOuterLayITSpureSA[iSpec]->Sumw2();
323 fOutput->Add(fHistOuterLayITSpureSA[iSpec]);
324
325 //---
326
327 hisname.Form("hPtResid%s",spname[iSpec].Data());
328 fHistPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
329 fHistPtResid[iSpec]->Sumw2();
330 fOutput->Add(fHistPtResid[iSpec]);
331 hisname.Form("hPtRelResid%s",spname[iSpec].Data());
332 fHistPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
333 fHistPtRelResid[iSpec]->Sumw2();
334 fOutput->Add(fHistPtRelResid[iSpec]);
335
336 hisname.Form("hInvPtResid%s",spname[iSpec].Data());
337 fHistInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
338 fHistInvPtResid[iSpec]->Sumw2();
339 fOutput->Add(fHistInvPtResid[iSpec]);
340 hisname.Form("hInvPtRelResid%s",spname[iSpec].Data());
341 fHistInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
342 fHistInvPtRelResid[iSpec]->Sumw2();
343 fOutput->Add(fHistInvPtRelResid[iSpec]);
344
345 hisname.Form("hMCPtResid%s",spname[iSpec].Data());
346 fHistMCPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
347 fHistMCPtResid[iSpec]->Sumw2();
348 fOutput->Add(fHistMCPtResid[iSpec]);
349 hisname.Form("hMCPtRelResid%s",spname[iSpec].Data());
350 fHistMCPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
351 fHistMCPtRelResid[iSpec]->Sumw2();
352 fOutput->Add(fHistMCPtRelResid[iSpec]);
353
354 hisname.Form("hMCInvPtResid%s",spname[iSpec].Data());
355 fHistMCInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
356 fHistMCInvPtResid[iSpec]->Sumw2();
357 fOutput->Add(fHistMCInvPtResid[iSpec]);
358 hisname.Form("hMCInvPtRelResid%s",spname[iSpec].Data());
359 fHistMCInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
360 fHistMCInvPtRelResid[iSpec]->Sumw2();
361 fOutput->Add(fHistMCInvPtRelResid[iSpec]);
362
363 }
364
2afdea69 365 fHistMCPhiResid=new TH2F("hMCPhiResid","",nbins,xbins,100,-0.5,0.5);
366 fHistMCPhiResid->Sumw2();
367 fOutput->Add(fHistMCPhiResid);
368 fHistPhiResid=new TH2F("hPhiResid","",nbins,xbins,100,-0.5,0.5);
369 fHistPhiResid->Sumw2();
370 fOutput->Add(fHistPhiResid);
371
372 PostData(1,fOutput);
373
045571a6 374}
375//______________________________________________________________________________
376void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
377{
378 //
379
380 AliESDEvent *esd = (AliESDEvent*) (InputEvent());
381
382
383 if(!esd) {
384 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
385 return;
386 }
387
388
389 if(!ESDfriend()) {
390 printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
391 return;
392 }
393
394 AliStack* stack=0;
395
396 if(fReadMC){
397 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
398 if (!eventHandler) {
399 Printf("ERROR: Could not retrieve MC event handler");
400 return;
401 }
402 AliMCEvent* mcEvent = eventHandler->MCEvent();
403 if (!mcEvent) {
404 Printf("ERROR: Could not retrieve MC event");
405 return;
406 }
407 stack = mcEvent->Stack();
408 if (!stack) {
409 Printf("ERROR: stack not available");
410 return;
411 }
412 }
413
045571a6 414
415 fHistNEvents->Fill(0);
416 const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
417 if(spdv->GetNContributors()<=0) return;
418 fHistNEvents->Fill(1);
2afdea69 419
420 const Int_t ntSize=36;
421 Float_t xnt[ntSize];
422
045571a6 423 Int_t ntracks = esd->GetNumberOfTracks();
424 for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
425 AliESDtrack * track = esd->GetTrack(iTrack);
426 if (!track) continue;
2afdea69 427 Float_t pttrack=track->Pt();
428 Float_t ptrack=track->P();
429 Float_t pxtrack=track->Px();
430 Float_t pytrack=track->Py();
431 Float_t pztrack=track->Pz();
432 Float_t impactXY=-999, impactZ=-999;
433 track->GetImpactParameters(impactXY, impactZ);
434 Double_t dedxlay[4];
435 track->GetITSdEdxSamples(dedxlay);
436 Float_t dedx=track->GetITSsignal();
437 Float_t chi2=track->GetITSchi2();
045571a6 438 Int_t status=track->GetStatus();
2afdea69 439 Bool_t isITSrefit=kFALSE;
440 Bool_t isTPCin=kFALSE;
441 Bool_t isTPCrefit=kFALSE;
442 Bool_t isPureSA=kFALSE;
443 if(status&AliESDtrack::kITSrefit) isITSrefit=kTRUE;
444 if(status&AliESDtrack::kTPCin) isTPCin=kTRUE;
445 if(status&AliESDtrack::kTPCrefit) isTPCrefit=kTRUE;
446 if(status&AliESDtrack::kITSpureSA) isPureSA=kTRUE;
045571a6 447 Bool_t isSA=kTRUE;
448 if(status&AliESDtrack::kTPCin) isSA=kFALSE;
449 Int_t nTPCclus=track->GetNcls(1);
450 Int_t nITSclus=track->GetNcls(0);
2afdea69 451 UChar_t clumap=track->GetITSClusterMap();
452 Int_t idet;
453 Float_t xloc,zloc;
454 Int_t statusLay[6];
455 for(Int_t iLay=0; iLay<6;iLay++) track->GetITSModuleIndexInfo(iLay,idet,statusLay[iLay],xloc,zloc);
456 Int_t trlabel=track->GetLabel();
457 Float_t ptgen=-999.;
458 Float_t pgen=-999.;
459 Float_t pxgen=-999.;
460 Float_t pygen=-999.;
461 Float_t pzgen=-999.;
462 Float_t etagen=-999.;
463 Float_t phigen=-999.;
464 if(fReadMC){
465 TParticle* part = stack->Particle(TMath::Abs(trlabel));
466 ptgen=part->Pt();
467 pgen=part->P();
468 pxgen=part->Px();
469 pygen=part->Py();
470 pzgen=part->Pz();
471 etagen=part->Eta();
472 phigen=part->Phi();
473 }
474
475 Int_t indexn=0;
476 xnt[indexn++]=pttrack;
477 xnt[indexn++]=pxtrack;
478 xnt[indexn++]=pytrack;
479 xnt[indexn++]=pztrack;
480 xnt[indexn++]=track->Eta();
481 xnt[indexn++]=track->Phi();
482 xnt[indexn++]=impactXY;
483 xnt[indexn++]=impactZ;
484 for(Int_t iLay=0; iLay<4; iLay++) xnt[indexn++]=dedxlay[iLay];
485 xnt[indexn++]=dedx;
486 xnt[indexn++]=chi2;
487 xnt[indexn++]=status;
488 xnt[indexn++]=isITSrefit;
489 xnt[indexn++]=isTPCin;
490 xnt[indexn++]=isTPCrefit;
491 xnt[indexn++]=isPureSA;
492 xnt[indexn++]=nITSclus;
493 xnt[indexn++]=nTPCclus;
494 xnt[indexn++]=clumap;
495 for(Int_t iLay=0; iLay<6; iLay++) xnt[indexn++]=(Float_t)statusLay[iLay];
496 xnt[indexn++]=trlabel;
497 xnt[indexn++]=ptgen;
498 xnt[indexn++]=pxgen;
499 xnt[indexn++]=pygen;
500 xnt[indexn++]=pzgen;
501 xnt[indexn++]=etagen;
502 xnt[indexn++]=phigen;
503 xnt[indexn++]=ntracks;
504
505 if(indexn>ntSize) printf("AliAnalysisTaskITSsaTracks: ERROR ntuple insexout of range\n");
506 if(fFillNtuple) fNtupleTracks->Fill(xnt);
507
508 if(!(status&AliESDtrack::kITSrefit)) continue;
045571a6 509 if(nITSclus<fMinITSpts)continue;
2afdea69 510 if((chi2/nITSclus) > fMaxITSChi2Clu) continue;
045571a6 511
045571a6 512 Int_t nSPDPoints=0;
513 for(Int_t i=0; i<2; i++){
514 if(clumap&(1<<i)) ++nSPDPoints;
515 }
516 if(nSPDPoints<fMinSPDpts) continue;
517
518 Int_t nPointsForPid=0;
519 for(Int_t i=2; i<6; i++){
520 if(clumap&(1<<i)) ++nPointsForPid;
521 }
2afdea69 522
523 Int_t iTrackType=-1;
524 if(status&AliESDtrack::kTPCin){
525 iTrackType=kTypeTPCITS;
045571a6 526 }else{
2afdea69 527 if(status&AliESDtrack::kITSpureSA) iTrackType=kTypeITSpureSA;
528 else iTrackType=kTypeITSsa;
529 }
530 if(iTrackType<0 || iTrackType>=kNtrackTypes) continue;
531 fHistPt[iTrackType]->Fill(pttrack);
532 fHistEtaPhi[iTrackType]->Fill(track->Eta(),track->Phi());
533 fHistChi2[iTrackType]->Fill(chi2/nITSclus);
534 fHistNclu[iTrackType]->Fill(nITSclus);
535 if(nPointsForPid==3){
536 fHistdedxvsPt3cls[iTrackType]->Fill(pttrack,dedx);
537 }
538 if(nPointsForPid==4){
539 fHistdedxvsPt4cls[iTrackType]->Fill(pttrack,dedx);
540 }
541 if(fReadMC){
542 if(trlabel<0){
543 fHistPtFake[iTrackType]->Fill(pttrack);
544 fHistEtaPhiFake[iTrackType]->Fill(track->Eta(),track->Phi());
545 fHistChi2Fake[iTrackType]->Fill(chi2/nITSclus);
546 fHistNcluFake[iTrackType]->Fill(nITSclus);
045571a6 547 }else{
2afdea69 548 fHistPtGood[iTrackType]->Fill(pttrack);
549 fHistEtaPhiGood[iTrackType]->Fill(track->Eta(),track->Phi());
550 fHistChi2Good[iTrackType]->Fill(chi2/nITSclus);
551 fHistNcluGood[iTrackType]->Fill(nITSclus);
045571a6 552 }
553 }
554
2afdea69 555 Int_t hadronSpecie=-1;
045571a6 556 if(fReadMC && fUseMCId){
045571a6 557 if(trlabel>=0){
2afdea69 558 TParticle* part = stack->Particle(trlabel);
559 Int_t pdg=TMath::Abs(part->GetPdgCode());
560 if(pdg==211) hadronSpecie=kPion;
561 else if(pdg==321) hadronSpecie=kKaon;
562 else if(pdg==2212) hadronSpecie=kProton;
045571a6 563 }
045571a6 564 }else{
2afdea69 565 if(nPointsForPid<fMinPtsforPid)continue;
045571a6 566 AliITSPIDResponse pidits(fReadMC);
567 Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
568 if(nSigmaPion>-2. && nSigmaPion<2.){
569 hadronSpecie=kPion;
570 }else{
571 Float_t nSigmaKaon=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kKaon,nPointsForPid,isSA);
572 if(nSigmaKaon>-2. && nSigmaKaon<2.){
573 hadronSpecie=kKaon;
574 }else{
575 Float_t nSigmaProton=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kProton,nPointsForPid,isSA);
576 if(nSigmaProton>-2. && nSigmaProton<2.){
577 hadronSpecie=kProton;
578 }
579 }
580 }
581 }
582 if(hadronSpecie<0) continue;
2afdea69 583 if(iTrackType==kTypeTPCITS){ // TPC+ITS tracks
045571a6 584 fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
585 fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
586 fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
587 fHistCluInLayTPCITS[hadronSpecie]->Fill(-1.);
588 for(Int_t iBit=0; iBit<6; iBit++){
589 if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
590 }
2afdea69 591 }else if(iTrackType==kTypeITSpureSA){ // TPC+ITS tracks
592 fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
593 fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
594 fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
595 fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
596 fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
597 fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
598 Int_t outerLay=-1;
599 for(Int_t iBit=0; iBit<6; iBit++){
600 if(clumap&(1<<iBit)){
601 fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
602 if(iBit>outerLay) outerLay=iBit;
045571a6 603 }
045571a6 604 }
2afdea69 605 fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);
606
045571a6 607
2afdea69 608 if(fReadMC){
609 if(trlabel>=0){
045571a6 610 TParticle* part = stack->Particle(trlabel);
2afdea69 611 ptgen=part->Pt();
045571a6 612 Float_t invpttrack=track->OneOverPt();
613 Float_t invptgen=0.;
614 if(ptgen>0.) invptgen=1./ptgen;
615 fHistMCPtResid[hadronSpecie]->Fill(pttrack,pttrack-ptgen);
616 fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
617 fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
618 fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);
2afdea69 619 Float_t deltaphi=track->Phi()-part->Phi();
620 if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
621 if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
622 fHistMCPhiResid->Fill(pttrack,deltaphi);
045571a6 623 }
2afdea69 624 }
625 }else if(iTrackType==kTypeITSsa){ // TPC+ITS tracks
626 fHistPtITSsa[hadronSpecie]->Fill(pttrack);
627 fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
628 fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
629 fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
630 for(Int_t iBit=0; iBit<6; iBit++){
631 if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
045571a6 632 }
633 }
2afdea69 634
635 // match TPCITS with ITSpureSA
636 if(nITSclus<fMinITSpts || nTPCclus<fMinTPCpts) continue;
045571a6 637 Int_t idxMI[12],idxSA[12];
638 for(Int_t icl=0; icl<12; icl++){
639 idxMI[icl]=-1;
640 idxSA[icl]=-1;
641 }
642 Int_t ncls=track->GetClusters(0,idxMI);
643 if(fMinITSpts<6){
644 Bool_t accept=kTRUE;
645 for(Int_t ilay=0; ilay<6; ilay++){
646 if(fRequirePoint[ilay]){
647 Int_t mask = 1<<ilay;
648 if(!(clumap & mask)){
649 accept=kFALSE;
650 break;
651 }
652 }
653 }
654 if(!accept) continue;
655 }
656 // Sort
657 for(Int_t i=0;i<12;i++){
658 for(Int_t j=i+1;j<12;j++){
659 if(idxMI[j]>idxMI[i]){
660 Int_t tmp=idxMI[j];
661 idxMI[j]=idxMI[i];
662 idxMI[i]=tmp;
663 }
664 }
665 }
666 // for(Int_t i=0; i<12; i++) printf("%d ",idxMI[i]);
667 // printf("\n");
668 if(idxMI[0]<0 && idxMI[0]==idxMI[1]) continue;
669 Bool_t matched=kFALSE;
670 for (Int_t iTrack2 = 0; iTrack2 < ntracks; iTrack2++) {
671 if(matched) break;
672 if(iTrack2==iTrack) continue;
673 AliESDtrack* track2 = esd->GetTrack(iTrack2);
674 Int_t status2=track2->GetStatus();
675 if(!(status2&AliESDtrack::kITSrefit)) continue;
676 if(!(status2&AliESDtrack::kITSpureSA)) continue;
677 Int_t clumap2=track2->GetITSClusterMap();
678 Int_t nITSclus2=track2->GetNcls(0);
679 Int_t nTPCclus2=track2->GetNcls(1);
680 if(nITSclus2<fMinITSpts || nTPCclus2>0) continue;
681 Int_t ncls2=track2->GetClusters(0,idxSA);
682 if(ncls2!=ncls) continue;
683 if(fMinITSpts<6){
684 Bool_t accept=kTRUE;
685 for(Int_t ilay=0; ilay<6; ilay++){
686 if(fRequirePoint[ilay]){
687 Int_t mask = 1<<ilay;
688 if(!(clumap2 & mask)){
689 accept=kFALSE;
690 break;
691 }
692 }
693 }
694 if(!accept) continue;
695 }
696 // Sort
697 for(Int_t i=0;i<12;i++){
698 for(Int_t j=i+1;j<12;j++){
699 if(idxSA[j]>idxSA[i]){
700 Int_t tmp=idxSA[j];
701 idxSA[j]=idxSA[i];
702 idxSA[i]=tmp;
703 }
704 }
705 }
706 Int_t match=0;
707 for(Int_t icl=0; icl<ncls; icl++){
708 if(idxSA[icl]!=idxMI[icl]){
709 match=0;
710 break;
711 }
712 else match++;
713 }
714 if(match==ncls && match>0){
715 matched=kTRUE;
716 Float_t pt1=track->Pt();
717 Float_t pt2=track2->Pt();
718 Float_t ptm1=track->OneOverPt();
719 Float_t ptm2=track2->OneOverPt();
720 fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
721 fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
722 fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
2afdea69 723 fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
724 Float_t deltaphi=track->Phi()-track2->Phi();
725 if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
726 if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
727 fHistPhiResid->Fill(pt1,deltaphi);
045571a6 728 }
729 }
730 }
731
732 PostData(1,fOutput);
733
734}
735//______________________________________________________________________________
736void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
737{
738 // Terminate analysis
739 fOutput = dynamic_cast<TList*> (GetOutputData(1));
740 if (!fOutput) {
741 printf("ERROR: fOutput not available\n");
742 return;
743 }
2afdea69 744 fNtupleTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("ntupleTracks"));
045571a6 745 fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
045571a6 746
2afdea69 747 TString spname[3]={"Pion","Kaon","Proton"};
045571a6 748
749 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
750 fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
751 fHistPtITSsa[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSsa%s",spname[iSpec].Data())));
752 fHistPtITSpureSA[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSpureSA%s",spname[iSpec].Data())));
753
754 fHistEtaPhiTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiTPCITS%s",spname[iSpec].Data())));
755 fHistEtaPhiITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSsa%s",spname[iSpec].Data())));
756 fHistEtaPhiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSpureSA%s",spname[iSpec].Data())));
757
758 fHistNcluTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluTPCITS%s",spname[iSpec].Data())));
759 fHistNcluITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSsa%s",spname[iSpec].Data())));
760 fHistNcluITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSpureSA%s",spname[iSpec].Data())));
761
762 fHistd0rphiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0rphiITSpureSA%s",spname[iSpec].Data())));
763 fHistd0zITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0zITSpureSA%s",spname[iSpec].Data())));
764
765 fHistCluInLayTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayTPCITS%s",spname[iSpec].Data())));
766 fHistCluInLayITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSsa%s",spname[iSpec].Data())));
767 fHistCluInLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSpureSA%s",spname[iSpec].Data())));
768
769 fHistOuterLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hOuterLayITSpureSA%s",spname[iSpec].Data())));
770
771 fHistPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtResid%s",spname[iSpec].Data())));
772 fHistPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtRelResid%s",spname[iSpec].Data())));
773 fHistInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtResid%s",spname[iSpec].Data())));
774 fHistInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtRelResid%s",spname[iSpec].Data())));
775 fHistMCPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtResid%s",spname[iSpec].Data())));
776 fHistMCPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtRelResid%s",spname[iSpec].Data())));
777 fHistMCInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtResid%s",spname[iSpec].Data())));
778 fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
779
780 }
2afdea69 781 fHistMCPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hMCPhiResid"));
782 fHistPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hPhiResid"));
045571a6 783 return;
784}
785
786
787
788
789