]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/ITS/AliAnalysisTaskITSsaTracks.cxx
fine tune split parameters, adding some control histograms
[u/mrichter/AliRoot.git] / PWGPP / 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 //
9c12af5f 80 if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
81 // RS: why do we delete all histos? they are owned by the output!!!
82 /*
2afdea69 83 delete fHistNEvents;
9c12af5f 84 for(Int_t iType=0; iType<kNtrackTypes; iType++){
2afdea69 85 delete fHistPt[iType];
86 delete fHistPtGood[iType];
87 delete fHistPtFake[iType];
88 delete fHistEtaPhi[iType];
89 delete fHistEtaPhiGood[iType];
90 delete fHistEtaPhiFake[iType];
519eed9b 91 delete fHistEtaPhiAny[iType];
92 delete fHistEtaPhi1SPD[iType];
93 delete fHistEtaPhi4Clu[iType];
94 delete fHistEtaPhi6Clu[iType];
2afdea69 95 delete fHistChi2[iType];
96 delete fHistChi2Good[iType];
97 delete fHistChi2Fake[iType];
98 delete fHistNclu[iType];
99 delete fHistNcluGood[iType];
100 delete fHistNcluFake[iType];
519eed9b 101 delete fHistdedxvsP2cls[iType];
102 delete fHistdedxvsP3cls[iType];
103 delete fHistdedxvsP4cls[iType];
2afdea69 104 }
105 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
106 delete fHistPtTPCITS[iSpec];
107 delete fHistPtITSsa[iSpec];
108 delete fHistPtITSpureSA[iSpec];
109 delete fHistEtaPhiTPCITS[iSpec];
110 delete fHistEtaPhiITSsa[iSpec];
111 delete fHistEtaPhiITSpureSA[iSpec];
112 delete fHistNcluTPCITS[iSpec];
113 delete fHistNcluITSsa[iSpec];
114 delete fHistNcluITSpureSA[iSpec];
115 delete fHistd0rphiITSpureSA[iSpec];
116 delete fHistd0zITSpureSA[iSpec];
117 delete fHistCluInLayTPCITS[iSpec];
118 delete fHistCluInLayITSsa[iSpec];
119 delete fHistCluInLayITSpureSA[iSpec];
120 delete fHistOuterLayITSpureSA[iSpec];
121 delete fHistPtResid[iSpec];
122 delete fHistPtRelResid[iSpec];
123 delete fHistInvPtResid[iSpec];
124 delete fHistInvPtRelResid[iSpec];
125 delete fHistMCPtResid[iSpec];
126 delete fHistMCPtRelResid[iSpec];
127 delete fHistMCInvPtResid[iSpec];
128 delete fHistMCInvPtRelResid[iSpec];
129 }
9c12af5f 130 */
2afdea69 131 delete fHistMCPhiResid;
132 delete fHistPhiResid;
133 delete fNtupleTracks;
045571a6 134 if (fOutput) {
135 delete fOutput;
136 fOutput = 0;
137 }
138}
139
140//________________________________________________________________________
141void AliAnalysisTaskITSsaTracks::SetPtBins(Int_t n, Double_t* lim){
142 // define pt bins for analysis
143 if(n>kMaxPtBins){
144 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
145 return;
146 }else{
147 fNPtBins=n;
148 for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
149 for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
150 }
151}
152//___________________________________________________________________________
153void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
154 // create output histos
155
156 fOutput = new TList();
157 fOutput->SetOwner();
158 fOutput->SetName("OutputHistos");
159
2afdea69 160 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");
161 // kinematics: pt, p eta,phi -> 6 variables
162 // impact parameters: d0xy, d0z -> 2 variables
163 // dE/dx: 4 Layers + trunc mean -> 5 variables
164 // chi2 and track status: -> 6 variables
165 // cluster infos: -> 9 variables
166 // MC info: -> 7 variables
167 // Multiplicity -> 1 variable
168 // Total: 36
169 fOutput->Add(fNtupleTracks);
170
045571a6 171 fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
172 fHistNEvents->Sumw2();
173 fHistNEvents->SetMinimum(0);
174 fOutput->Add(fHistNEvents);
175
176 //binning for the dedx histogram
177 const Int_t hnbinsdedx=400;
178 Double_t hxmindedx = 0.01;
179 Double_t hxmaxdedx = 10;
180 Double_t hlogxmindedx = TMath::Log10(hxmindedx);
181 Double_t hlogxmaxdedx = TMath::Log10(hxmaxdedx);
182 Double_t hbinwidthdedx = (hlogxmaxdedx-hlogxmindedx)/hnbinsdedx;
183 Double_t hxbinsdedx[hnbinsdedx+1];
184 hxbinsdedx[0] = 0.01;
185 for (Int_t i=1;i<=hnbinsdedx;i++) {
186 hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
187 }
188
2afdea69 189 TString tit[kNtrackTypes]={"TPCITS","ITSsa","ITSpureSA"};
190 for(Int_t iType=0; iType<kNtrackTypes; iType++){
191 fHistPt[iType]=new TH1F(Form("hPt%s",tit[iType].Data()),"",100,0.,2.);
192 fHistPt[iType]->Sumw2();
193 fOutput->Add(fHistPt[iType]);
194 fHistPtGood[iType]=new TH1F(Form("hPtGood%s",tit[iType].Data()),"",100,0.,2.);
195 fHistPtGood[iType]->Sumw2();
196 fOutput->Add(fHistPtGood[iType]);
197 fHistPtFake[iType]=new TH1F(Form("hPtFake%s",tit[iType].Data()),"",100,0.,2.);
198 fHistPtFake[iType]->Sumw2();
199 fOutput->Add(fHistPtFake[iType]);
200
201 fHistEtaPhi[iType] = new TH2F(Form("hEtaPhi%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
202 fHistEtaPhi[iType]->Sumw2();
203 fOutput->Add(fHistEtaPhi[iType]);
204 fHistEtaPhiGood[iType] = new TH2F(Form("hEtaPhiGood%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
205 fHistEtaPhiGood[iType]->Sumw2();
206 fOutput->Add(fHistEtaPhiGood[iType]);
207 fHistEtaPhiFake[iType] = new TH2F(Form("hEtaPhiFake%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
208 fHistEtaPhiFake[iType]->Sumw2();
209 fOutput->Add(fHistEtaPhiFake[iType]);
210
519eed9b 211 fHistEtaPhiAny[iType] = new TH2F(Form("hEtaPhiAny%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
212 fHistEtaPhiAny[iType]->Sumw2();
213 fOutput->Add(fHistEtaPhiAny[iType]);
214 fHistEtaPhi1SPD[iType] = new TH2F(Form("hEtaPhi1SPD%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
215 fHistEtaPhi1SPD[iType]->Sumw2();
216 fOutput->Add(fHistEtaPhi1SPD[iType]);
217 fHistEtaPhi4Clu[iType] = new TH2F(Form("hEtaPhi4Clu%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
218 fHistEtaPhi4Clu[iType]->Sumw2();
219 fOutput->Add(fHistEtaPhi4Clu[iType]);
220 fHistEtaPhi6Clu[iType] = new TH2F(Form("hEtaPhi6Clu%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
221 fHistEtaPhi6Clu[iType]->Sumw2();
222 fOutput->Add(fHistEtaPhi6Clu[iType]);
223
2afdea69 224 fHistChi2[iType]=new TH1F(Form("hChi2%s",tit[iType].Data()),"",100,0.,10.);
225 fHistChi2[iType]->Sumw2();
226 fOutput->Add(fHistChi2[iType]);
227 fHistChi2Good[iType]=new TH1F(Form("hChi2Good%s",tit[iType].Data()),"",100,0.,10.);
228 fHistChi2Good[iType]->Sumw2();
229 fOutput->Add(fHistChi2Good[iType]);
230 fHistChi2Fake[iType]=new TH1F(Form("hChi2Fake%s",tit[iType].Data()),"",100,0.,10.);
231 fHistChi2Fake[iType]->Sumw2();
232 fOutput->Add(fHistChi2Fake[iType]);
233
234 fHistNclu[iType]=new TH1F(Form("hNclu%s",tit[iType].Data()),"",7,-0.5,6.5);
235 fHistNclu[iType]->Sumw2();
236 fOutput->Add(fHistNclu[iType]);
237 fHistNcluGood[iType]=new TH1F(Form("hNcluGood%s",tit[iType].Data()),"",7,-0.5,6.5);
238 fHistNcluGood[iType]->Sumw2();
239 fOutput->Add(fHistNcluGood[iType]);
240 fHistNcluFake[iType]=new TH1F(Form("hNcluFake%s",tit[iType].Data()),"",7,-0.5,6.5);
241 fHistNcluFake[iType]->Sumw2();
242 fOutput->Add(fHistNcluFake[iType]);
243
519eed9b 244 fHistdedxvsP2cls[iType] = new TH2F(Form("hdedxvsP2cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
245 fHistdedxvsP2cls[iType]->Sumw2();
246 fOutput->Add(fHistdedxvsP2cls[iType]);
247
248 fHistdedxvsP3cls[iType] = new TH2F(Form("hdedxvsP3cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
249 fHistdedxvsP3cls[iType]->Sumw2();
250 fOutput->Add(fHistdedxvsP3cls[iType]);
2afdea69 251
519eed9b 252 fHistdedxvsP4cls[iType] = new TH2F(Form("hdedxvsP4cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
253 fHistdedxvsP4cls[iType]->Sumw2();
254 fOutput->Add(fHistdedxvsP4cls[iType]);
2afdea69 255 }
256
257 //-----------------------------------------------------------
258
045571a6 259 TString spname[3]={"Pion","Kaon","Proton"};
260 TString hisname;
261 const Int_t nbins = fNPtBins;
262 Double_t xbins[nbins+1];
263 for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
264
045571a6 265 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
266
267 hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
268 fHistPtTPCITS[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
269 fHistPtTPCITS[iSpec]->Sumw2();
270 fOutput->Add(fHistPtTPCITS[iSpec]);
271
272 hisname.Form("hPtITSsa%s",spname[iSpec].Data());
273 fHistPtITSsa[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
274 fHistPtITSsa[iSpec]->Sumw2();
275 fOutput->Add(fHistPtITSsa[iSpec]);
276
277 hisname.Form("hPtITSpureSA%s",spname[iSpec].Data());
278 fHistPtITSpureSA[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
279 fHistPtITSpureSA[iSpec]->Sumw2();
280 fOutput->Add(fHistPtITSpureSA[iSpec]);
281
282 //---
283
284 hisname.Form("hEtaPhiTPCITS%s",spname[iSpec].Data());
285 fHistEtaPhiTPCITS[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
286 fHistEtaPhiTPCITS[iSpec]->Sumw2();
287 fOutput->Add(fHistEtaPhiTPCITS[iSpec]);
288
289 hisname.Form("hEtaPhiITSsa%s",spname[iSpec].Data());
290 fHistEtaPhiITSsa[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
291 fHistEtaPhiITSsa[iSpec]->Sumw2();
292 fOutput->Add(fHistEtaPhiITSsa[iSpec]);
293
294 hisname.Form("hEtaPhiITSpureSA%s",spname[iSpec].Data());
295 fHistEtaPhiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
296 fHistEtaPhiITSpureSA[iSpec]->Sumw2();
297 fOutput->Add(fHistEtaPhiITSpureSA[iSpec]);
298
299 //---
300
301 hisname.Form("hNcluTPCITS%s",spname[iSpec].Data());
302 fHistNcluTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
303 fHistNcluTPCITS[iSpec]->Sumw2();
304 fOutput->Add(fHistNcluTPCITS[iSpec]);
305
306 hisname.Form("hNcluITSsa%s",spname[iSpec].Data());
307 fHistNcluITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
308 fHistNcluITSsa[iSpec]->Sumw2();
309 fOutput->Add(fHistNcluITSsa[iSpec]);
310
311 hisname.Form("hNcluITSpureSA%s",spname[iSpec].Data());
312 fHistNcluITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
313 fHistNcluITSpureSA[iSpec]->Sumw2();
314 fOutput->Add(fHistNcluITSpureSA[iSpec]);
315
316 //---
317
318 hisname.Form("hd0rphiITSpureSA%s",spname[iSpec].Data());
319 fHistd0rphiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
320 fHistd0rphiITSpureSA[iSpec]->Sumw2();
321 fOutput->Add(fHistd0rphiITSpureSA[iSpec]);
322
323 hisname.Form("hd0zITSpureSA%s",spname[iSpec].Data());
324 fHistd0zITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
325 fHistd0zITSpureSA[iSpec]->Sumw2();
326 fOutput->Add(fHistd0zITSpureSA[iSpec]);
327
328 //---
329
330 hisname.Form("hCluInLayTPCITS%s",spname[iSpec].Data());
331 fHistCluInLayTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
332 fHistCluInLayTPCITS[iSpec]->Sumw2();
333 fOutput->Add(fHistCluInLayTPCITS[iSpec]);
334
335 hisname.Form("hCluInLayITSsa%s",spname[iSpec].Data());
336 fHistCluInLayITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
337 fHistCluInLayITSsa[iSpec]->Sumw2();
338 fOutput->Add(fHistCluInLayITSsa[iSpec]);
339
340 hisname.Form("hCluInLayITSpureSA%s",spname[iSpec].Data());
341 fHistCluInLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
342 fHistCluInLayITSpureSA[iSpec]->Sumw2();
343 fOutput->Add(fHistCluInLayITSpureSA[iSpec]);
344
345 hisname.Form("hOuterLayITSpureSA%s",spname[iSpec].Data());
346 fHistOuterLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
347 fHistOuterLayITSpureSA[iSpec]->Sumw2();
348 fOutput->Add(fHistOuterLayITSpureSA[iSpec]);
349
350 //---
351
352 hisname.Form("hPtResid%s",spname[iSpec].Data());
353 fHistPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
354 fHistPtResid[iSpec]->Sumw2();
355 fOutput->Add(fHistPtResid[iSpec]);
356 hisname.Form("hPtRelResid%s",spname[iSpec].Data());
357 fHistPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
358 fHistPtRelResid[iSpec]->Sumw2();
359 fOutput->Add(fHistPtRelResid[iSpec]);
360
361 hisname.Form("hInvPtResid%s",spname[iSpec].Data());
362 fHistInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
363 fHistInvPtResid[iSpec]->Sumw2();
364 fOutput->Add(fHistInvPtResid[iSpec]);
365 hisname.Form("hInvPtRelResid%s",spname[iSpec].Data());
366 fHistInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
367 fHistInvPtRelResid[iSpec]->Sumw2();
368 fOutput->Add(fHistInvPtRelResid[iSpec]);
369
370 hisname.Form("hMCPtResid%s",spname[iSpec].Data());
371 fHistMCPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
372 fHistMCPtResid[iSpec]->Sumw2();
373 fOutput->Add(fHistMCPtResid[iSpec]);
374 hisname.Form("hMCPtRelResid%s",spname[iSpec].Data());
375 fHistMCPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
376 fHistMCPtRelResid[iSpec]->Sumw2();
377 fOutput->Add(fHistMCPtRelResid[iSpec]);
378
379 hisname.Form("hMCInvPtResid%s",spname[iSpec].Data());
380 fHistMCInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
381 fHistMCInvPtResid[iSpec]->Sumw2();
382 fOutput->Add(fHistMCInvPtResid[iSpec]);
383 hisname.Form("hMCInvPtRelResid%s",spname[iSpec].Data());
384 fHistMCInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
385 fHistMCInvPtRelResid[iSpec]->Sumw2();
386 fOutput->Add(fHistMCInvPtRelResid[iSpec]);
387
388 }
389
2afdea69 390 fHistMCPhiResid=new TH2F("hMCPhiResid","",nbins,xbins,100,-0.5,0.5);
391 fHistMCPhiResid->Sumw2();
392 fOutput->Add(fHistMCPhiResid);
393 fHistPhiResid=new TH2F("hPhiResid","",nbins,xbins,100,-0.5,0.5);
394 fHistPhiResid->Sumw2();
395 fOutput->Add(fHistPhiResid);
396
397 PostData(1,fOutput);
398
045571a6 399}
400//______________________________________________________________________________
401void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
402{
403 //
404
405 AliESDEvent *esd = (AliESDEvent*) (InputEvent());
406
407
408 if(!esd) {
409 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
410 return;
411 }
412
413
414 if(!ESDfriend()) {
415 printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
416 return;
417 }
418
419 AliStack* stack=0;
420
421 if(fReadMC){
422 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
423 if (!eventHandler) {
424 Printf("ERROR: Could not retrieve MC event handler");
425 return;
426 }
427 AliMCEvent* mcEvent = eventHandler->MCEvent();
428 if (!mcEvent) {
429 Printf("ERROR: Could not retrieve MC event");
430 return;
431 }
432 stack = mcEvent->Stack();
433 if (!stack) {
434 Printf("ERROR: stack not available");
435 return;
436 }
437 }
438
045571a6 439
440 fHistNEvents->Fill(0);
441 const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
442 if(spdv->GetNContributors()<=0) return;
443 fHistNEvents->Fill(1);
2afdea69 444
445 const Int_t ntSize=36;
446 Float_t xnt[ntSize];
447
045571a6 448 Int_t ntracks = esd->GetNumberOfTracks();
449 for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
450 AliESDtrack * track = esd->GetTrack(iTrack);
451 if (!track) continue;
2afdea69 452 Float_t pttrack=track->Pt();
453 Float_t ptrack=track->P();
454 Float_t pxtrack=track->Px();
455 Float_t pytrack=track->Py();
456 Float_t pztrack=track->Pz();
457 Float_t impactXY=-999, impactZ=-999;
458 track->GetImpactParameters(impactXY, impactZ);
459 Double_t dedxlay[4];
460 track->GetITSdEdxSamples(dedxlay);
461 Float_t dedx=track->GetITSsignal();
462 Float_t chi2=track->GetITSchi2();
045571a6 463 Int_t status=track->GetStatus();
2afdea69 464 Bool_t isITSrefit=kFALSE;
465 Bool_t isTPCin=kFALSE;
466 Bool_t isTPCrefit=kFALSE;
467 Bool_t isPureSA=kFALSE;
468 if(status&AliESDtrack::kITSrefit) isITSrefit=kTRUE;
469 if(status&AliESDtrack::kTPCin) isTPCin=kTRUE;
470 if(status&AliESDtrack::kTPCrefit) isTPCrefit=kTRUE;
471 if(status&AliESDtrack::kITSpureSA) isPureSA=kTRUE;
045571a6 472 Bool_t isSA=kTRUE;
473 if(status&AliESDtrack::kTPCin) isSA=kFALSE;
474 Int_t nTPCclus=track->GetNcls(1);
475 Int_t nITSclus=track->GetNcls(0);
2afdea69 476 UChar_t clumap=track->GetITSClusterMap();
477 Int_t idet;
478 Float_t xloc,zloc;
479 Int_t statusLay[6];
480 for(Int_t iLay=0; iLay<6;iLay++) track->GetITSModuleIndexInfo(iLay,idet,statusLay[iLay],xloc,zloc);
481 Int_t trlabel=track->GetLabel();
482 Float_t ptgen=-999.;
483 Float_t pgen=-999.;
484 Float_t pxgen=-999.;
485 Float_t pygen=-999.;
486 Float_t pzgen=-999.;
487 Float_t etagen=-999.;
488 Float_t phigen=-999.;
489 if(fReadMC){
490 TParticle* part = stack->Particle(TMath::Abs(trlabel));
491 ptgen=part->Pt();
492 pgen=part->P();
493 pxgen=part->Px();
494 pygen=part->Py();
495 pzgen=part->Pz();
496 etagen=part->Eta();
497 phigen=part->Phi();
498 }
499
500 Int_t indexn=0;
501 xnt[indexn++]=pttrack;
502 xnt[indexn++]=pxtrack;
503 xnt[indexn++]=pytrack;
504 xnt[indexn++]=pztrack;
505 xnt[indexn++]=track->Eta();
506 xnt[indexn++]=track->Phi();
507 xnt[indexn++]=impactXY;
508 xnt[indexn++]=impactZ;
509 for(Int_t iLay=0; iLay<4; iLay++) xnt[indexn++]=dedxlay[iLay];
510 xnt[indexn++]=dedx;
511 xnt[indexn++]=chi2;
512 xnt[indexn++]=status;
513 xnt[indexn++]=isITSrefit;
514 xnt[indexn++]=isTPCin;
515 xnt[indexn++]=isTPCrefit;
516 xnt[indexn++]=isPureSA;
517 xnt[indexn++]=nITSclus;
518 xnt[indexn++]=nTPCclus;
519 xnt[indexn++]=clumap;
520 for(Int_t iLay=0; iLay<6; iLay++) xnt[indexn++]=(Float_t)statusLay[iLay];
521 xnt[indexn++]=trlabel;
522 xnt[indexn++]=ptgen;
523 xnt[indexn++]=pxgen;
524 xnt[indexn++]=pygen;
525 xnt[indexn++]=pzgen;
526 xnt[indexn++]=etagen;
527 xnt[indexn++]=phigen;
528 xnt[indexn++]=ntracks;
529
530 if(indexn>ntSize) printf("AliAnalysisTaskITSsaTracks: ERROR ntuple insexout of range\n");
531 if(fFillNtuple) fNtupleTracks->Fill(xnt);
532
533 if(!(status&AliESDtrack::kITSrefit)) continue;
519eed9b 534
535 Int_t iTrackType=-1;
536 if(status&AliESDtrack::kTPCin){
537 iTrackType=kTypeTPCITS;
538 }else{
539 if(status&AliESDtrack::kITSpureSA) iTrackType=kTypeITSpureSA;
540 else iTrackType=kTypeITSsa;
541 }
542 if(iTrackType<0 || iTrackType>=kNtrackTypes) continue;
045571a6 543
045571a6 544 Int_t nSPDPoints=0;
545 for(Int_t i=0; i<2; i++){
546 if(clumap&(1<<i)) ++nSPDPoints;
547 }
045571a6 548 Int_t nPointsForPid=0;
549 for(Int_t i=2; i<6; i++){
550 if(clumap&(1<<i)) ++nPointsForPid;
551 }
2afdea69 552
519eed9b 553 fHistEtaPhiAny[iTrackType]->Fill(track->Eta(),track->Phi());
554 if(nSPDPoints>=1) fHistEtaPhi1SPD[iTrackType]->Fill(track->Eta(),track->Phi());
555 if(nSPDPoints>=1 && nPointsForPid>=3) fHistEtaPhi4Clu[iTrackType]->Fill(track->Eta(),track->Phi());
556 if(nITSclus==6) fHistEtaPhi6Clu[iTrackType]->Fill(track->Eta(),track->Phi());
557
558
559 if(nITSclus<fMinITSpts)continue;
560 if((chi2/nITSclus) > fMaxITSChi2Clu) continue;
561 if(nSPDPoints<fMinSPDpts) continue;
562
2afdea69 563 fHistPt[iTrackType]->Fill(pttrack);
564 fHistEtaPhi[iTrackType]->Fill(track->Eta(),track->Phi());
519eed9b 565
2afdea69 566 fHistChi2[iTrackType]->Fill(chi2/nITSclus);
567 fHistNclu[iTrackType]->Fill(nITSclus);
519eed9b 568 if(nPointsForPid==2){
569 fHistdedxvsP2cls[iTrackType]->Fill(ptrack,dedx);
570 }
571 else if(nPointsForPid==3){
572 fHistdedxvsP3cls[iTrackType]->Fill(ptrack,dedx);
2afdea69 573 }
519eed9b 574 else if(nPointsForPid==4){
575 fHistdedxvsP4cls[iTrackType]->Fill(ptrack,dedx);
2afdea69 576 }
577 if(fReadMC){
578 if(trlabel<0){
579 fHistPtFake[iTrackType]->Fill(pttrack);
580 fHistEtaPhiFake[iTrackType]->Fill(track->Eta(),track->Phi());
581 fHistChi2Fake[iTrackType]->Fill(chi2/nITSclus);
582 fHistNcluFake[iTrackType]->Fill(nITSclus);
045571a6 583 }else{
2afdea69 584 fHistPtGood[iTrackType]->Fill(pttrack);
585 fHistEtaPhiGood[iTrackType]->Fill(track->Eta(),track->Phi());
586 fHistChi2Good[iTrackType]->Fill(chi2/nITSclus);
587 fHistNcluGood[iTrackType]->Fill(nITSclus);
045571a6 588 }
589 }
590
2afdea69 591 Int_t hadronSpecie=-1;
045571a6 592 if(fReadMC && fUseMCId){
045571a6 593 if(trlabel>=0){
2afdea69 594 TParticle* part = stack->Particle(trlabel);
595 Int_t pdg=TMath::Abs(part->GetPdgCode());
596 if(pdg==211) hadronSpecie=kPion;
597 else if(pdg==321) hadronSpecie=kKaon;
598 else if(pdg==2212) hadronSpecie=kProton;
045571a6 599 }
045571a6 600 }else{
2afdea69 601 if(nPointsForPid<fMinPtsforPid)continue;
045571a6 602 AliITSPIDResponse pidits(fReadMC);
603 Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
604 if(nSigmaPion>-2. && nSigmaPion<2.){
605 hadronSpecie=kPion;
606 }else{
607 Float_t nSigmaKaon=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kKaon,nPointsForPid,isSA);
608 if(nSigmaKaon>-2. && nSigmaKaon<2.){
609 hadronSpecie=kKaon;
610 }else{
611 Float_t nSigmaProton=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kProton,nPointsForPid,isSA);
612 if(nSigmaProton>-2. && nSigmaProton<2.){
613 hadronSpecie=kProton;
614 }
615 }
616 }
617 }
618 if(hadronSpecie<0) continue;
2afdea69 619 if(iTrackType==kTypeTPCITS){ // TPC+ITS tracks
045571a6 620 fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
621 fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
622 fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
1686659d 623 fHistCluInLayTPCITS[hadronSpecie]->Fill(-1.,-1.);
045571a6 624 for(Int_t iBit=0; iBit<6; iBit++){
625 if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
626 }
2afdea69 627 }else if(iTrackType==kTypeITSpureSA){ // TPC+ITS tracks
628 fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
629 fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
630 fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
631 fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
632 fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
1686659d 633 fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.,-1.);
2afdea69 634 Int_t outerLay=-1;
635 for(Int_t iBit=0; iBit<6; iBit++){
636 if(clumap&(1<<iBit)){
637 fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
638 if(iBit>outerLay) outerLay=iBit;
045571a6 639 }
045571a6 640 }
2afdea69 641 fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);
642
045571a6 643
2afdea69 644 if(fReadMC){
645 if(trlabel>=0){
045571a6 646 TParticle* part = stack->Particle(trlabel);
2afdea69 647 ptgen=part->Pt();
045571a6 648 Float_t invpttrack=track->OneOverPt();
649 Float_t invptgen=0.;
650 if(ptgen>0.) invptgen=1./ptgen;
651 fHistMCPtResid[hadronSpecie]->Fill(pttrack,pttrack-ptgen);
652 fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
653 fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
654 fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);
2afdea69 655 Float_t deltaphi=track->Phi()-part->Phi();
656 if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
657 if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
658 fHistMCPhiResid->Fill(pttrack,deltaphi);
045571a6 659 }
2afdea69 660 }
661 }else if(iTrackType==kTypeITSsa){ // TPC+ITS tracks
662 fHistPtITSsa[hadronSpecie]->Fill(pttrack);
663 fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
664 fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
1686659d 665 fHistCluInLayITSsa[hadronSpecie]->Fill(-1.,-1.);
2afdea69 666 for(Int_t iBit=0; iBit<6; iBit++){
667 if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
045571a6 668 }
669 }
2afdea69 670
671 // match TPCITS with ITSpureSA
672 if(nITSclus<fMinITSpts || nTPCclus<fMinTPCpts) continue;
045571a6 673 Int_t idxMI[12],idxSA[12];
674 for(Int_t icl=0; icl<12; icl++){
675 idxMI[icl]=-1;
676 idxSA[icl]=-1;
677 }
678 Int_t ncls=track->GetClusters(0,idxMI);
679 if(fMinITSpts<6){
680 Bool_t accept=kTRUE;
681 for(Int_t ilay=0; ilay<6; ilay++){
682 if(fRequirePoint[ilay]){
683 Int_t mask = 1<<ilay;
684 if(!(clumap & mask)){
685 accept=kFALSE;
686 break;
687 }
688 }
689 }
690 if(!accept) continue;
691 }
692 // Sort
693 for(Int_t i=0;i<12;i++){
694 for(Int_t j=i+1;j<12;j++){
695 if(idxMI[j]>idxMI[i]){
696 Int_t tmp=idxMI[j];
697 idxMI[j]=idxMI[i];
698 idxMI[i]=tmp;
699 }
700 }
701 }
702 // for(Int_t i=0; i<12; i++) printf("%d ",idxMI[i]);
703 // printf("\n");
704 if(idxMI[0]<0 && idxMI[0]==idxMI[1]) continue;
705 Bool_t matched=kFALSE;
706 for (Int_t iTrack2 = 0; iTrack2 < ntracks; iTrack2++) {
707 if(matched) break;
708 if(iTrack2==iTrack) continue;
709 AliESDtrack* track2 = esd->GetTrack(iTrack2);
710 Int_t status2=track2->GetStatus();
711 if(!(status2&AliESDtrack::kITSrefit)) continue;
712 if(!(status2&AliESDtrack::kITSpureSA)) continue;
713 Int_t clumap2=track2->GetITSClusterMap();
714 Int_t nITSclus2=track2->GetNcls(0);
715 Int_t nTPCclus2=track2->GetNcls(1);
716 if(nITSclus2<fMinITSpts || nTPCclus2>0) continue;
717 Int_t ncls2=track2->GetClusters(0,idxSA);
718 if(ncls2!=ncls) continue;
719 if(fMinITSpts<6){
720 Bool_t accept=kTRUE;
721 for(Int_t ilay=0; ilay<6; ilay++){
722 if(fRequirePoint[ilay]){
723 Int_t mask = 1<<ilay;
724 if(!(clumap2 & mask)){
725 accept=kFALSE;
726 break;
727 }
728 }
729 }
730 if(!accept) continue;
731 }
732 // Sort
733 for(Int_t i=0;i<12;i++){
734 for(Int_t j=i+1;j<12;j++){
735 if(idxSA[j]>idxSA[i]){
736 Int_t tmp=idxSA[j];
737 idxSA[j]=idxSA[i];
738 idxSA[i]=tmp;
739 }
740 }
741 }
742 Int_t match=0;
743 for(Int_t icl=0; icl<ncls; icl++){
744 if(idxSA[icl]!=idxMI[icl]){
745 match=0;
746 break;
747 }
748 else match++;
749 }
750 if(match==ncls && match>0){
751 matched=kTRUE;
752 Float_t pt1=track->Pt();
753 Float_t pt2=track2->Pt();
754 Float_t ptm1=track->OneOverPt();
755 Float_t ptm2=track2->OneOverPt();
756 fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
757 fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
758 fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
2afdea69 759 fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
760 Float_t deltaphi=track->Phi()-track2->Phi();
761 if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
762 if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
763 fHistPhiResid->Fill(pt1,deltaphi);
045571a6 764 }
765 }
766 }
767
768 PostData(1,fOutput);
769
770}
771//______________________________________________________________________________
772void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
773{
774 // Terminate analysis
775 fOutput = dynamic_cast<TList*> (GetOutputData(1));
776 if (!fOutput) {
777 printf("ERROR: fOutput not available\n");
778 return;
779 }
2afdea69 780 fNtupleTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("ntupleTracks"));
045571a6 781 fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
045571a6 782
2afdea69 783 TString spname[3]={"Pion","Kaon","Proton"};
045571a6 784
785 for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
786 fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
787 fHistPtITSsa[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSsa%s",spname[iSpec].Data())));
788 fHistPtITSpureSA[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSpureSA%s",spname[iSpec].Data())));
789
790 fHistEtaPhiTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiTPCITS%s",spname[iSpec].Data())));
791 fHistEtaPhiITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSsa%s",spname[iSpec].Data())));
792 fHistEtaPhiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSpureSA%s",spname[iSpec].Data())));
793
794 fHistNcluTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluTPCITS%s",spname[iSpec].Data())));
795 fHistNcluITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSsa%s",spname[iSpec].Data())));
796 fHistNcluITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSpureSA%s",spname[iSpec].Data())));
797
798 fHistd0rphiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0rphiITSpureSA%s",spname[iSpec].Data())));
799 fHistd0zITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0zITSpureSA%s",spname[iSpec].Data())));
800
801 fHistCluInLayTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayTPCITS%s",spname[iSpec].Data())));
802 fHistCluInLayITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSsa%s",spname[iSpec].Data())));
803 fHistCluInLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSpureSA%s",spname[iSpec].Data())));
804
805 fHistOuterLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hOuterLayITSpureSA%s",spname[iSpec].Data())));
806
807 fHistPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtResid%s",spname[iSpec].Data())));
808 fHistPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtRelResid%s",spname[iSpec].Data())));
809 fHistInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtResid%s",spname[iSpec].Data())));
810 fHistInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtRelResid%s",spname[iSpec].Data())));
811 fHistMCPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtResid%s",spname[iSpec].Data())));
812 fHistMCPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtRelResid%s",spname[iSpec].Data())));
813 fHistMCInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtResid%s",spname[iSpec].Data())));
814 fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
815
816 }
2afdea69 817 fHistMCPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hMCPhiResid"));
818 fHistPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hPhiResid"));
045571a6 819 return;
820}
821
822
823
824
825