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