]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliPWG4HighPtTrackQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Change in JETAN from AliEMCALGeoUtils to AliEMCALGeometry (M. Bregant)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtTrackQA.cxx
CommitLineData
d756027f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// This class stores QA variables as function of pT for different type
18// of tracks and track selection criteria
19// Output: Histograms for different set of cuts
20//-----------------------------------------------------------------------
21// Author : Marta Verweij - UU
22//-----------------------------------------------------------------------
23
24#ifndef ALIPWG4HIGHPTTRACKQA_CXX
25#define ALIPWG4HIGHPTTRACKQA_CXX
26
27#include "AliPWG4HighPtTrackQA.h"
28
29#include "TH1.h"
30#include "TH2.h"
31#include "TH3.h"
32#include "TProfile.h"
33#include "TList.h"
34#include "TFile.h"
35#include "TChain.h"
36#include "TH3F.h"
37#include "TKey.h"
38#include "TSystem.h"
39#include "TBits.h"
40
41#include "AliAnalysisManager.h"
42#include "AliESDInputHandler.h"
43#include "AliMCEvent.h"
44#include "AliMCEventHandler.h"
45#include "AliStack.h"
46#include "AliESDtrack.h"
47#include "AliESDtrackCuts.h"
48#include "AliExternalTrackParam.h"
49#include "AliLog.h"
50#include "AliGenPythiaEventHeader.h"
51#include "AliGenCocktailEventHeader.h"
52#include "AliCentrality.h"
53#include "AliAODVertex.h"
54#include "AliAODEvent.h"
55//#include "AliAnalysisHelperJetTasks.h"
56
57using namespace std; //required for resolving the 'cout' symbol
58
59ClassImp(AliPWG4HighPtTrackQA)
60
61AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA()
62: AliAnalysisTaskSE(),
63 fDataType(kESD),
64 fEvent(0x0),
65 fESD(0x0),
66 fVtx(0x0),
44684f3b 67 fTrackCuts(0x0),
68 fTrackCutsITSLoose(0x0),
69 fTrackCutsTPConly(0x0),
d756027f 70 fTrackType(0),
71 fFilterMask(0),
aa3ba8d2 72 fSigmaConstrainedMax(-1.),
d756027f 73 fPtMax(100.),
74 fIsPbPb(0),
75 fCentClass(10),
2b553e6f 76 fNVariables(18),
d756027f 77 fVariables(0x0),
78 fAvgTrials(1),
79 fNEventAll(0),
80 fNEventSel(0),
81 fNEventReject(0),
82 fh1Centrality(0x0),
83 fh1Xsec(0),
84 fh1Trials(0),
85 fh1PtHard(0),
86 fh1PtHardTrials(0),
87 fh1NTracksAll(0x0),
88 fh1NTracksReject(0x0),
89 fh1NTracksSel(0x0),
90 fPtAll(0),
91 fPtSel(0),
92 fPtPhi(0x0),
93 fPtEta(0x0),
94 fPtDCA2D(0x0),
95 fPtDCAZ(0x0),
96 fPtNClustersTPC(0x0),
97 fPtNPointITS(0x0),
98 fPtChi2C(0x0),
99 fPtNSigmaToVertex(0x0),
100 fPtRelUncertainty1Pt(0x0),
2b553e6f 101 fPtUncertainty1Pt(0x0),
d756027f 102 fPtChi2PerClusterTPC(0x0),
103 fPtNCrossedRows(0x0),
104 fPtNCrossedRowsNClusF(0x0),
105 fPtNCrRNCrRNClusF(0x0),
2b553e6f 106 fPtSigmaY2(0x0),
107 fPtSigmaZ2(0x0),
108 fPtSigmaSnp2(0x0),
109 fPtSigmaTgl2(0x0),
110 fPtSigma1Pt2(0x0),
aa3ba8d2 111 fProfPtSigmaY2(0x0),
112 fProfPtSigmaZ2(0x0),
113 fProfPtSigmaSnp2(0x0),
114 fProfPtSigmaTgl2(0x0),
115 fProfPtSigma1Pt2(0x0),
116 fProfPtSigma1Pt(0x0),
117 fProfPtPtSigma1Pt(0x0),
118 fSystTrackCuts(0x0),
d756027f 119 fHistList(0)
120{
2b553e6f 121 SetNVariables(18);
aa3ba8d2 122
123 fPtBinEdges[0][0] = 10.;
124 fPtBinEdges[0][1] = 1.;
125 fPtBinEdges[1][0] = 20.;
126 fPtBinEdges[1][1] = 2.;
127 fPtBinEdges[2][0] = 50.;
128 fPtBinEdges[2][1] = 5.;
129
d756027f 130}
131//________________________________________________________________________
132AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name):
133 AliAnalysisTaskSE(name),
134 fDataType(kESD),
135 fEvent(0x0),
136 fESD(0x0),
137 fVtx(0x0),
44684f3b 138 fTrackCuts(0x0),
139 fTrackCutsITSLoose(0x0),
140 fTrackCutsTPConly(0x0),
d756027f 141 fTrackType(0),
142 fFilterMask(0),
aa3ba8d2 143 fSigmaConstrainedMax(-1.),
d756027f 144 fPtMax(100.),
145 fIsPbPb(0),
146 fCentClass(10),
2b553e6f 147 fNVariables(18),
d756027f 148 fVariables(0x0),
149 fAvgTrials(1),
150 fNEventAll(0),
151 fNEventSel(0),
152 fNEventReject(0),
153 fh1Centrality(0x0),
154 fh1Xsec(0),
155 fh1Trials(0),
156 fh1PtHard(0),
157 fh1PtHardTrials(0),
158 fh1NTracksAll(0x0),
159 fh1NTracksReject(0x0),
160 fh1NTracksSel(0x0),
161 fPtAll(0),
162 fPtSel(0),
163 fPtPhi(0x0),
164 fPtEta(0x0),
165 fPtDCA2D(0x0),
166 fPtDCAZ(0x0),
167 fPtNClustersTPC(0x0),
168 fPtNPointITS(0x0),
169 fPtChi2C(0x0),
170 fPtNSigmaToVertex(0x0),
171 fPtRelUncertainty1Pt(0x0),
2b553e6f 172 fPtUncertainty1Pt(0x0),
d756027f 173 fPtChi2PerClusterTPC(0x0),
174 fPtNCrossedRows(0x0),
175 fPtNCrossedRowsNClusF(0x0),
176 fPtNCrRNCrRNClusF(0x0),
2b553e6f 177 fPtSigmaY2(0x0),
178 fPtSigmaZ2(0x0),
179 fPtSigmaSnp2(0x0),
180 fPtSigmaTgl2(0x0),
181 fPtSigma1Pt2(0x0),
aa3ba8d2 182 fProfPtSigmaY2(0x0),
183 fProfPtSigmaZ2(0x0),
184 fProfPtSigmaSnp2(0x0),
185 fProfPtSigmaTgl2(0x0),
186 fProfPtSigma1Pt2(0x0),
187 fProfPtSigma1Pt(0x0),
188 fProfPtPtSigma1Pt(0x0),
189 fSystTrackCuts(0x0),
d756027f 190 fHistList(0)
191{
192 //
193 // Constructor. Initialization of Inputs and Outputs
194 //
195 AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
196
2b553e6f 197 SetNVariables(18);
d756027f 198
aa3ba8d2 199 fPtBinEdges[0][0] = 10.;
200 fPtBinEdges[0][1] = 1.;
201 fPtBinEdges[1][0] = 20.;
202 fPtBinEdges[1][1] = 2.;
203 fPtBinEdges[2][0] = 50.;
204 fPtBinEdges[2][1] = 5.;
205
d756027f 206 // Input slot #0 works with a TChain ESD
207 DefineInput(0, TChain::Class());
208 // Output slot #1 write into a TList
209 DefineOutput(1, TList::Class());
210}
211
aa3ba8d2 212//________________________________________________________________________
213void AliPWG4HighPtTrackQA::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
214
215 if(region<3) {
216 fPtBinEdges[region][0] = ptmax;
217 fPtBinEdges[region][1] = ptBinWidth;
218 }
219 else {
220 AliError("Only 3 regions alowed. Use region 0/1/2\n");
221 return;
222 }
223
224}
225
d756027f 226//________________________________________________________________________
227void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
228 //Create output objects
229 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
230
231 Bool_t oldStatus = TH1::AddDirectoryStatus();
232 TH1::AddDirectory(kFALSE);
233
234 OpenFile(1);
235 fHistList = new TList();
236 fHistList->SetOwner(kTRUE);
237
238 Float_t fgkPtMin = 0.;
239 Float_t fgkPtMax = fPtMax;
240
aa3ba8d2 241 // Float_t ptBinEdges[2][2];
242 // ptBinEdges[0][0] = 10.;
243 // ptBinEdges[0][1] = 1.;
244 // ptBinEdges[1][0] = 20.;
245 // ptBinEdges[1][1] = 2.;
246 // Float_t binWidth3 = 5.;
247 // if(fPtMax>100.) {
248 // ptBinEdges[0][0] = 100.;
249 // ptBinEdges[0][1] = 5.;
250 // ptBinEdges[1][0] = 300.;
251 // ptBinEdges[1][1] = 10.;
252 // binWidth3 = 20.;
253 // }
254
255 //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
d756027f 256 const Float_t ptmin1 = fgkPtMin;
aa3ba8d2 257 const Float_t ptmax1 = fPtBinEdges[0][0];
d756027f 258 const Float_t ptmin2 = ptmax1 ;
aa3ba8d2 259 const Float_t ptmax2 = fPtBinEdges[1][0];
d756027f 260 const Float_t ptmin3 = ptmax2 ;
261 const Float_t ptmax3 = fgkPtMax;
aa3ba8d2 262 const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
263 const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
2793fe00 264 const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
d756027f 265 Int_t fgkNPtBins=nbin13;
266 //Create array with low edges of each bin
267 Double_t *binsPt=new Double_t[fgkNPtBins+1];
268 for(Int_t i=0; i<=fgkNPtBins; i++) {
269 if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
270 if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
271 if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
272 }
273
274 Int_t fgkNPhiBins = 18*6;
275 Float_t kMinPhi = 0.;
276 Float_t kMaxPhi = 2.*TMath::Pi();
277 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
278 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
279
280 Int_t fgkNEtaBins=20;
281 Float_t fgkEtaMin = -1.;
282 Float_t fgkEtaMax = 1.;
283 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
284 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
285
286 Int_t fgkNNClustersTPCBins=80;
287 Float_t fgkNClustersTPCMin = 0.5;
288 Float_t fgkNClustersTPCMax = 160.5;
289 Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
290 for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
291
292 Int_t fgkNDCA2DBins=80;
293 Float_t fgkDCA2DMin = -0.2;
294 Float_t fgkDCA2DMax = 0.2;
aa3ba8d2 295 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
2b553e6f 296 fgkDCA2DMin = -2.;
297 fgkDCA2DMax = 2.;
298 }
d756027f 299 Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
300 for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
301
302 Int_t fgkNDCAZBins=80;
303 Float_t fgkDCAZMin = -2.;
304 Float_t fgkDCAZMax = 2.;
aa3ba8d2 305 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
2b553e6f 306 fgkDCAZMin = -5.;
307 fgkDCAZMax = 5.;
308 }
d756027f 309 Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
310 for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
311
312 Int_t fgkNNPointITSBins=9;
313 Float_t fgkNPointITSMin = -0.5;
314 Float_t fgkNPointITSMax = 8.5;
315 Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
316 for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
317
318 Int_t fgkNNSigmaToVertexBins=40;
319 Float_t fgkNSigmaToVertexMin = 0.;
320 Float_t fgkNSigmaToVertexMax = 8.;
321 Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
322 for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
323
324 Int_t fgkNChi2CBins=20;
325 Float_t fgkChi2CMin = 0.;
2b553e6f 326 Float_t fgkChi2CMax = 100.;
d756027f 327 Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
328 for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
329
330 Int_t fgkNRel1PtUncertaintyBins=30;
331 Float_t fgkRel1PtUncertaintyMin = 0.;
332 Float_t fgkRel1PtUncertaintyMax = 0.3;
aa3ba8d2 333 if(fTrackType==1 || fTrackType==2 || fTrackType==4)
334 fgkRel1PtUncertaintyMax = 0.5;
d756027f 335 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
336 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
337
2b553e6f 338 Int_t fgkNUncertainty1PtBins = 30;
339 Float_t fgkUncertainty1PtMin = 0.;
340 Float_t fgkUncertainty1PtMax = 0.1;
aa3ba8d2 341 if(fTrackType==1 || fTrackType==2 || fTrackType==4)
342 fgkUncertainty1PtMax = 0.2;
2b553e6f 343 Double_t *binsUncertainty1Pt=new Double_t[fgkNUncertainty1PtBins+1];
344 for(Int_t i=0; i<=fgkNUncertainty1PtBins; i++) binsUncertainty1Pt[i]=(Double_t)fgkUncertainty1PtMin + (fgkUncertainty1PtMax-fgkUncertainty1PtMin)/fgkNUncertainty1PtBins*(Double_t)i ;
345
d756027f 346 Float_t fgkChi2PerClusMin = 0.;
347 Float_t fgkChi2PerClusMax = 4.;
348 Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
349 Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
350 for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
351
352 Int_t fgkNCrossedRowsNClusFBins = 50;
353 Float_t fgkNCrossedRowsNClusFMin = 0.;
354 Float_t fgkNCrossedRowsNClusFMax = 1.;
355 Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
356 for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
357
2b553e6f 358 Float_t fgk1PtMin = 0.;
359 Float_t fgk1PtMax = 6.;
aa3ba8d2 360 Float_t binEdge1Pt1 = 1.;
361 Float_t binWidth1Pt1 = 0.05;
362 Int_t fgkN1PtBins1 = (int)((binEdge1Pt1-fgk1PtMin)/binWidth1Pt1);
363 Float_t binWidth1Pt2 = 0.1;
364 Int_t fgkN1PtBins2 = (int)((fgk1PtMax-binEdge1Pt1)/binWidth1Pt2);
365 Int_t fgkN1PtBins = fgkN1PtBins1+fgkN1PtBins2;
2b553e6f 366 Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
aa3ba8d2 367
368 for(Int_t i=0; i<=fgkN1PtBins; i++) {
369 if(i<=fgkN1PtBins1)
370 bins1Pt[i]=(Double_t)fgk1PtMin + (Double_t)(binEdge1Pt1-fgk1PtMin)/(Double_t)fgkN1PtBins1*(Double_t)i;
371 if(i<=fgkN1PtBins && i>fgkN1PtBins1)
372 bins1Pt[i]=(Double_t)binEdge1Pt1 + (Double_t)(fgk1PtMax-binEdge1Pt1)/(Double_t)fgkN1PtBins2*(Double_t)(i-fgkN1PtBins1);
373 }
2b553e6f 374
375 Int_t fgkNSigmaY2Bins = 50;
376 Float_t fgkSigmaY2Min = 0.;
aa3ba8d2 377 Float_t fgkSigmaY2Max = 1.;
378 if(fTrackType==1) fgkSigmaY2Max = 4.;
379 if(fTrackType==2 || fTrackType==4) fgkSigmaY2Max = 0.1;
2b553e6f 380 Double_t *binsSigmaY2=new Double_t[fgkNSigmaY2Bins+1];
381 for(Int_t i=0; i<=fgkNSigmaY2Bins; i++) binsSigmaY2[i]=(Double_t)fgkSigmaY2Min + (fgkSigmaY2Max-fgkSigmaY2Min)/fgkNSigmaY2Bins*(Double_t)i ;
382
383 Int_t fgkNSigmaZ2Bins = 50;
384 Float_t fgkSigmaZ2Min = 0.;
aa3ba8d2 385 Float_t fgkSigmaZ2Max = 0.4;
2b553e6f 386 Double_t *binsSigmaZ2=new Double_t[fgkNSigmaZ2Bins+1];
387 for(Int_t i=0; i<=fgkNSigmaZ2Bins; i++) binsSigmaZ2[i]=(Double_t)fgkSigmaZ2Min + (fgkSigmaZ2Max-fgkSigmaZ2Min)/fgkNSigmaZ2Bins*(Double_t)i ;
388
389 Int_t fgkNSigmaSnp2Bins = 50;
390 Float_t fgkSigmaSnp2Min = 0.;
aa3ba8d2 391 Float_t fgkSigmaSnp2Max = 0.05;
392 if(fTrackType==1) fgkSigmaSnp2Max = 0.2;
393 if(fTrackType==2 || fTrackType==4) fgkSigmaSnp2Max = 0.1;
2b553e6f 394 Double_t *binsSigmaSnp2=new Double_t[fgkNSigmaSnp2Bins+1];
395 for(Int_t i=0; i<=fgkNSigmaSnp2Bins; i++) binsSigmaSnp2[i]=(Double_t)fgkSigmaSnp2Min + (fgkSigmaSnp2Max-fgkSigmaSnp2Min)/fgkNSigmaSnp2Bins*(Double_t)i ;
396
397 Int_t fgkNSigmaTgl2Bins = 50;
398 Float_t fgkSigmaTgl2Min = 0.;
aa3ba8d2 399 Float_t fgkSigmaTgl2Max = 0.1;
400 if(fTrackType==1) fgkSigmaTgl2Max = 0.2;
401 if(fTrackType==2 || fTrackType==4) fgkSigmaTgl2Max = 0.1;
2b553e6f 402 Double_t *binsSigmaTgl2=new Double_t[fgkNSigmaTgl2Bins+1];
403 for(Int_t i=0; i<=fgkNSigmaTgl2Bins; i++) binsSigmaTgl2[i]=(Double_t)fgkSigmaTgl2Min + (fgkSigmaTgl2Max-fgkSigmaTgl2Min)/fgkNSigmaTgl2Bins*(Double_t)i ;
404
405 Int_t fgkNSigma1Pt2Bins = 50;
406 Float_t fgkSigma1Pt2Min = 0.;
aa3ba8d2 407 Float_t fgkSigma1Pt2Max = 1.;
2b553e6f 408 Double_t *binsSigma1Pt2=new Double_t[fgkNSigma1Pt2Bins+1];
409 for(Int_t i=0; i<=fgkNSigma1Pt2Bins; i++) binsSigma1Pt2[i]=(Double_t)fgkSigma1Pt2Min + (fgkSigma1Pt2Max-fgkSigma1Pt2Min)/fgkNSigma1Pt2Bins*(Double_t)i ;
410
411
d756027f 412 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
413 fHistList->Add(fNEventAll);
414 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
415 fHistList->Add(fNEventSel);
416 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
417 //Set labels
418 fNEventReject->Fill("noESD",0);
419 fNEventReject->Fill("Trigger",0);
420 fNEventReject->Fill("NTracks<2",0);
421 fNEventReject->Fill("noVTX",0);
422 fNEventReject->Fill("VtxStatus",0);
423 fNEventReject->Fill("NCont<2",0);
424 fNEventReject->Fill("ZVTX>10",0);
425 fNEventReject->Fill("cent",0);
426 fNEventReject->Fill("cent>90",0);
427 fHistList->Add(fNEventReject);
428
429 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
430 fHistList->Add(fh1Centrality);
431
432 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
433 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
434 fHistList->Add(fh1Xsec);
435
436 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
437 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
438 fHistList->Add(fh1Trials);
439
440 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
441 fHistList->Add(fh1PtHard);
442 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
443 fHistList->Add(fh1PtHardTrials);
444
445 fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
446 fHistList->Add(fh1NTracksAll);
447
448 fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
449 fh1NTracksReject->Fill("noESDtrack",0);
450 fh1NTracksReject->Fill("noTPCInner",0);
451 fh1NTracksReject->Fill("FillTPC",0);
452 fh1NTracksReject->Fill("noTPConly",0);
453 fh1NTracksReject->Fill("relate",0);
454 fh1NTracksReject->Fill("trackCuts",0);
455 fh1NTracksReject->Fill("laser",0);
2b553e6f 456 fh1NTracksReject->Fill("chi2",0);
d756027f 457 fHistList->Add(fh1NTracksReject);
458
459 fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
460 fHistList->Add(fh1NTracksSel);
461
aa3ba8d2 462 fSystTrackCuts = new TH1F("fSystTrackCuts","fSystTrackCuts",1,-0.5,0.5);
463 fSystTrackCuts->Fill("noCut",0);
464 fSystTrackCuts->Fill("eta",0);
465 fSystTrackCuts->Fill("0.15<pT<1e10",0);
466 fSystTrackCuts->Fill("kink",0);
467 fSystTrackCuts->Fill("NClusterTPC",0);
468 fSystTrackCuts->Fill("Chi2PerNClusTPC",0);
469 fSystTrackCuts->Fill("DCA2D",0);
470 fHistList->Add(fSystTrackCuts);
471
d756027f 472
473 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
474 fHistList->Add(fPtAll);
475 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
476 fHistList->Add(fPtSel);
477
478 fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
479 fHistList->Add(fPtPhi);
480
481 fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
482 fHistList->Add(fPtEta);
483
484 fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
485 fHistList->Add(fPtDCA2D);
486
487 fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
488 fHistList->Add(fPtDCAZ);
489
490 fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
491 fHistList->Add(fPtNClustersTPC);
492
493 fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
494 fHistList->Add(fPtNPointITS);
495
496 fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
497 fHistList->Add(fPtChi2C);
498
499 fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
500 fHistList->Add(fPtNSigmaToVertex);
501
502 fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
503 fHistList->Add(fPtRelUncertainty1Pt);
2b553e6f 504
505 fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
506 fHistList->Add(fPtUncertainty1Pt);
d756027f 507
508 fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
509 fHistList->Add(fPtChi2PerClusterTPC);
510
511 fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
512 fHistList->Add(fPtNCrossedRows);
513
514 fPtNCrossedRowsNClusF = new TH2F("fPtNCrossedRowsNClusF","fPtNCrossedRowsNClusF",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
515 fHistList->Add(fPtNCrossedRowsNClusF);
516
517 fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
518 fHistList->Add(fPtNCrRNCrRNClusF);
2b553e6f 519
aa3ba8d2 520 fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
2b553e6f 521 fHistList->Add(fPtSigmaY2);
522
523 fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
524 fHistList->Add(fPtSigmaZ2);
525
aa3ba8d2 526 fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
2b553e6f 527 fHistList->Add(fPtSigmaSnp2);
528
aa3ba8d2 529 fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
2b553e6f 530 fHistList->Add(fPtSigmaTgl2);
531
aa3ba8d2 532 fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
2b553e6f 533 fHistList->Add(fPtSigma1Pt2);
534
aa3ba8d2 535 fProfPtSigmaY2 = new TProfile("fProfPtSigmaY2","fProfPtSigmaY2",fgkN1PtBins,bins1Pt);
536 fHistList->Add(fProfPtSigmaY2);
537
538 fProfPtSigmaZ2 = new TProfile("fProfPtSigmaZ2","fProfPtSigmaZ2",fgkN1PtBins,bins1Pt);
539 fHistList->Add(fProfPtSigmaZ2);
540
541 fProfPtSigmaSnp2 = new TProfile("fProfPtSigmaSnp2","fProfPtSigmaSnp2",fgkN1PtBins,bins1Pt);
542 fHistList->Add(fProfPtSigmaSnp2);
543
544 fProfPtSigmaTgl2 = new TProfile("fProfPtSigmaTgl2","fProfPtSigmaTgl2",fgkN1PtBins,bins1Pt);
545 fHistList->Add(fProfPtSigmaTgl2);
546
547 fProfPtSigma1Pt2 = new TProfile("fProfPtSigma1Pt2","fProfPtSigma1Pt2",fgkN1PtBins,bins1Pt);
548 fHistList->Add(fProfPtSigma1Pt2);
549
550 fProfPtSigma1Pt = new TProfile("fProfPtSigma1Pt","fProfPtSigma1Pt;p_{T};#sigma(1/p_{T})",fgkNPtBins,binsPt);
551 fHistList->Add(fProfPtSigma1Pt);
552
553 fProfPtPtSigma1Pt = new TProfile("fProfPtPtSigma1Pt","fProfPtPtSigma1Pt;p_{T};p_{T}#sigma(1/p_{T})",fgkNPtBins,binsPt);
554 fHistList->Add(fProfPtPtSigma1Pt);
555
d756027f 556 TH1::AddDirectory(oldStatus);
557
558 PostData(1, fHistList);
559
560 if(binsPhi) delete [] binsPhi;
561 if(binsPt) delete [] binsPt;
562 if(binsNClustersTPC) delete [] binsNClustersTPC;
563 if(binsDCA2D) delete [] binsDCA2D;
564 if(binsDCAZ) delete [] binsDCAZ;
565 if(binsNPointITS) delete [] binsNPointITS;
566 if(binsNSigmaToVertex) delete [] binsNSigmaToVertex;
567 if(binsChi2C) delete [] binsChi2C;
568 if(binsEta) delete [] binsEta;
569 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
2b553e6f 570 if(binsUncertainty1Pt) delete [] binsUncertainty1Pt;
d756027f 571 if(binsChi2PerClus) delete [] binsChi2PerClus;
572 if(binsChi2PerClus) delete [] binsNCrossedRowsNClusF;
2b553e6f 573 if(bins1Pt) delete [] bins1Pt;
574 if(binsSigmaY2) delete [] binsSigmaY2;
575 if(binsSigmaZ2) delete [] binsSigmaZ2;
576 if(binsSigmaSnp2) delete [] binsSigmaSnp2;
577 if(binsSigmaTgl2) delete [] binsSigmaTgl2;
578 if(binsSigma1Pt2) delete [] binsSigma1Pt2;
d756027f 579}
580
581//________________________________________________________________________
582Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
583 //
584 // Decide if event should be selected for analysis
585 //
586
587 // Checks following requirements:
588 // - fEvent available
589 // - trigger info from AliPhysicsSelection
590 // - MCevent available
591 // - number of reconstructed tracks > 1
592 // - primary vertex reconstructed
593 // - z-vertex < 10 cm
594 // - centrality in case of PbPb
595
596 Bool_t selectEvent = kTRUE;
597
598 //fEvent object available?
599 if (!fEvent) {
600 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
601 fNEventReject->Fill("noAliVEvent",1);
602 selectEvent = kFALSE;
603 return selectEvent;
604 }
605
606 //Check if number of reconstructed tracks is larger than 1
607 if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2) {
608 fNEventReject->Fill("NTracks<2",1);
609 selectEvent = kFALSE;
610 return selectEvent;
611 }
612
613 //Check if vertex is reconstructed
1ea145ef 614 if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
615 fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 616
617 if(!fVtx) {
618 fNEventReject->Fill("noVTX",1);
619 selectEvent = kFALSE;
620 return selectEvent;
621 }
622
623 if(!fVtx->GetStatus()) {
624 fNEventReject->Fill("VtxStatus",1);
625 selectEvent = kFALSE;
626 return selectEvent;
627 }
628
629 // Need vertex cut
630 if(fVtx->GetNContributors()<2) {
631 fNEventReject->Fill("NCont<2",1);
632 selectEvent = kFALSE;
633 return selectEvent;
634 }
635
636 //Check if z-vertex < 10 cm
637 double primVtx[3];
638 fVtx->GetXYZ(primVtx);
639 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
640 fNEventReject->Fill("ZVTX>10",1);
641 selectEvent = kFALSE;
642 return selectEvent;
643 }
644 }
1ea145ef 645 else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
646 const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 647 if(!vtx) {
648 fNEventReject->Fill("noVTX",1);
649 selectEvent = kFALSE;
650 return selectEvent;
651 }
652
653 // Need vertex cut
654 if(vtx->GetNContributors()<2) {
655 fNEventReject->Fill("NCont<2",1);
656 selectEvent = kFALSE;
657 return selectEvent;
658 }
659
660 //Check if z-vertex < 10 cm
661 double primVtx[3];
662 vtx->GetXYZ(primVtx);
663 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
664 fNEventReject->Fill("ZVTX>10",1);
665 selectEvent = kFALSE;
666 return selectEvent;
667 }
668
669 }
670
671 //Centrality selection should only be done in case of PbPb
672 if(IsPbPb()) {
673 Float_t cent = 0.;
674 if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
675 fNEventReject->Fill("cent",1);
676 selectEvent = kFALSE;
677 return selectEvent;
678 }
679 else {
680 if(fDataType==kESD) {
681 if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
682 cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
683 }
684 }
685 else if(fDataType==kAOD) {
686 if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
687 cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
688 }
689 if(cent>90.) {
690 fNEventReject->Fill("cent>90",1);
691 selectEvent = kFALSE;
692 return selectEvent;
693 }
694 fh1Centrality->Fill(cent);
695 }
696 }
697
698 return selectEvent;
699
700}
701
702//________________________________________________________________________
703Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
669e2312 704 //
705 // Get centrality from ESD or AOD
706 //
707
d756027f 708 if(fDataType==kESD)
709 return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
710 else if(fDataType==kAOD)
711 return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
712 else
713 return 5;
714}
715
716//________________________________________________________________________
717Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
669e2312 718 //
719 // Get centrality from ESD
720 //
d756027f 721
669e2312 722 Float_t cent = -1;
1ea145ef 723
724 if(esd){
725 if(esd->GetCentrality()){
726 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
669e2312 727 if(fDebug>3) printf("centrality: %f\n",cent);
1ea145ef 728 }
d756027f 729 }
730
669e2312 731 return GetCentralityClass(cent);
d756027f 732
733}
734
735//________________________________________________________________________
736Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
669e2312 737 //
738 // Get centrality from AOD
739 //
d756027f 740
669e2312 741 if(!aod) return 5;
d756027f 742 Float_t cent = aod->GetHeader()->GetCentrality();
669e2312 743 if(fDebug>3) printf("centrality: %f\n",cent);
744
745 return GetCentralityClass(cent);
746
747}
748
749//________________________________________________________________________
750Int_t AliPWG4HighPtTrackQA::GetCentralityClass(Float_t cent) {
751 //
752 // Get centrality class
753 //
754
755 if(cent<0) return 5; // OB - cent sometimes negative
756 if(cent>80) return 4;
757 if(cent>50) return 3;
758 if(cent>30) return 2;
759 if(cent>10) return 1;
d756027f 760 return 0;
761
762}
763
764//________________________________________________________________________
765void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
766 // Main loop
767 // Called for each event
768 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
769
770 fEvent = InputEvent();
771 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
772
773 // All events without selection
774 fNEventAll->Fill(0.);
775
776 if(!SelectEvent()) {
777 // Post output data
778 PostData(1, fHistList);
779 return;
780 }
781
782
783 //Need to keep track of selected events
784 fNEventSel->Fill(0.);
785
786 fVariables = new TArrayF(fNVariables);
787
788 if(fDataType==kESD) DoAnalysisESD();
789 if(fDataType==kAOD) DoAnalysisAOD();
790
791 //Delete old fVariables
792 if(fVariables) delete fVariables;
793
794 // Post output data
795 PostData(1, fHistList);
796
797}
798
799//________________________________________________________________________
800void AliPWG4HighPtTrackQA::DoAnalysisESD() {
801
802 if(!fESD) {
803 PostData(1, fHistList);
804 return;
805 }
806
807 // ---- Get MC Header information (for MC productions in pThard bins) ----
808 Double_t ptHard = 0.;
809 Double_t nTrials = 1; // trials for MC trigger weight for real data
810
811 AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
812 if (eventHandlerMC) {
813
814 if(eventHandlerMC->MCEvent()){
815 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
816 if(pythiaGenHeader){
817 nTrials = pythiaGenHeader->Trials();
818 ptHard = pythiaGenHeader->GetPtHard();
819
820 fh1PtHard->Fill(ptHard);
821 fh1PtHardTrials->Fill(ptHard,nTrials);
822
823 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
824 }
825 }
826 }
827
828 Int_t nTracks = fESD->GetNumberOfTracks();
829 AliDebug(2,Form("nTracks ESD%d", nTracks));
830
831 /*
832 Variables to be put in fVariables
833 0: pt
834 1: phi
835 2: eta
836 3: dca2D
837 4: dcaZ
838 5: nClustersTPC
839 6: nPointITS
840 7: chi2C
841 8: nSigmaToVertex
842 9: relUncertainty1Pt
843 10: chi2PerClusterTPC
844 11: #crossed rows
845 12: (#crossed rows)/(#findable clusters)
2b553e6f 846 13: SigmaY2
847 14: SigmaZ2
848 15: SigmaSnp2
849 16: SigmaTgl2
850 17: Sigma1Pt2
d756027f 851 */
852
853 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
854 fh1NTracksAll->Fill(0.);
855
856 //Get track for analysis
2b553e6f 857 AliESDtrack *track = 0x0;
d756027f 858 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
859 if(!esdtrack) {
860 fh1NTracksReject->Fill("noESDtrack",1);
861 continue;
862 }
863
aa3ba8d2 864 if(fTrackType==4) {
865 FillSystematicCutHist(esdtrack);
866 if (!(fTrackCuts->AcceptTrack(esdtrack))) {
867 fh1NTracksReject->Fill("trackCuts",1);
868 continue;
869 }
870 }
871
d756027f 872 if(fTrackType==1)
873 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
aa3ba8d2 874 else if(fTrackType==2 || fTrackType==4) {
d756027f 875 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
876 if(!track) {
877 fh1NTracksReject->Fill("noTPConly",1);
d756027f 878 continue;
879 }
880 AliExternalTrackParam exParam;
881 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
882 if( !relate ) {
883 fh1NTracksReject->Fill("relate",1);
fb4a2fc7 884 if(track) delete track;
d756027f 885 continue;
886 }
887 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
888 }
44684f3b 889 else if(fTrackType==5 || fTrackType==6) {
890 if(fTrackCuts->AcceptTrack(esdtrack)) {
891 continue;
892 }
893 else {
894 if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
895
896 if(fTrackType==5) {
897 //use TPConly constrained track
898 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
899 AliExternalTrackParam exParam;
900 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
901 if( !relate ) {
902 fh1NTracksReject->Fill("relate",1);
fb4a2fc7 903 if(track) delete track;
44684f3b 904 continue;
905 }
906 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
907 }
908 else if(fTrackType==6) {
909 //use global constrained track
fb4a2fc7 910 track = esdtrack;
44684f3b 911 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
912
913 }
914 }
915 }
916 }
d756027f 917 else
918 track = esdtrack;
919
920 if(!track) {
d756027f 921 continue;
922 }
923
fb4a2fc7 924 if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
2b553e6f 925 //Cut on chi2 of constrained fit
aa3ba8d2 926 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
2b553e6f 927 fh1NTracksReject->Fill("chi2",1);
fb4a2fc7 928 if(track) delete track;
2b553e6f 929 continue;
930 }
931 }
932
aa3ba8d2 933 if(fTrackType!=4)
934 FillSystematicCutHist(track);
935
d756027f 936 fPtAll->Fill(track->Pt());
937
44684f3b 938 if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
d756027f 939 fh1NTracksReject->Fill("trackCuts",1);
fb4a2fc7 940 if(fTrackType==1 || fTrackType==2) {
941 if(track) delete track;
942 }
d756027f 943 continue;
944 }
aa3ba8d2 945
d756027f 946 fh1NTracksSel->Fill(0.);
947
948 fVariables->Reset(0.);
949
2b553e6f 950 fVariables->SetAt(track->Pt(),0);
951 fVariables->SetAt(track->Phi(),1);
952 fVariables->SetAt(track->Eta(),2);
d756027f 953
954 Float_t dca2D = 0.;
955 Float_t dcaz = 0.;
956 if(fTrackType==0) { //Global
957 track->GetImpactParameters(dca2D,dcaz);
958 }
aa3ba8d2 959 else if(fTrackType==1 || fTrackType==2 || fTrackType==4) { //TPConly
d756027f 960 track->GetImpactParametersTPC(dca2D,dcaz);
961 }
962 fVariables->SetAt(dca2D,3);
aa3ba8d2 963 fVariables->SetAt(dcaz,5);
d756027f 964
965 fVariables->SetAt((float)track->GetTPCNcls(),5);
966
967 Int_t nPointITS = 0;
968 UChar_t itsMap = track->GetITSClusterMap();
969 for (Int_t i=0; i < 6; i++) {
970 if (itsMap & (1 << i))
971 nPointITS ++;
972 }
973 fVariables->SetAt((float)nPointITS,6);
2b553e6f 974 Float_t chi2C = (float)track->GetConstrainedChi2();
aa3ba8d2 975 if(fTrackType==1 || fTrackType==2 || fTrackType==4)
2b553e6f 976 chi2C = (float)track->GetConstrainedChi2TPC();
977 fVariables->SetAt(chi2C,7);
d756027f 978 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
979
980 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
981
982 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
983
d756027f 984 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
44684f3b 985 Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
986 if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
d756027f 987 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
2b553e6f 988 fVariables->SetAt(track->GetSigmaY2(),13);
989 fVariables->SetAt(track->GetSigmaZ2(),14);
990 fVariables->SetAt(track->GetSigmaSnp2(),15);
991 fVariables->SetAt(track->GetSigmaTgl2(),16);
992 fVariables->SetAt(track->GetSigma1Pt2(),17);
d756027f 993
994 FillHistograms();
995
996 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
997
fb4a2fc7 998 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5) {
999 if(track) delete track;
1000 }
d756027f 1001
1002 }//track loop
1003
1004}
1005
1006//________________________________________________________________________
1007void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1008
1009 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1ea145ef 1010 if(!aod)return;
d756027f 1011 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1012
1013 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1014 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1015
1016 fVariables->Reset(0.);
1017
1018 fVariables->SetAt(aodtrack->Pt(),0);
1019 fVariables->SetAt(aodtrack->Phi(),1);
1020 fVariables->SetAt(aodtrack->Eta(),2);
1021
1022 Double_t dca[2] = {1e6,1e6};
1023 Double_t covar[3] = {1e6,1e6,1e6};
1024 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1025 fVariables->SetAt(dca[0],3);
1026 fVariables->SetAt(dca[1],4);
1027 }
1028
1029 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1030 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
aa3ba8d2 1031 fVariables->SetAt(aodtrack->Chi2perNDF(),7);
d756027f 1032 fVariables->SetAt(0.,8);
1033 fVariables->SetAt(0.,9);
1034 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1035 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
1036 Float_t crossedRowsTPCNClsF = 0.;
1037 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1038 fVariables->SetAt(crossedRowsTPCNClsF,12);
1039
1040 fPtAll->Fill(fVariables->At(0));
1041
1042 FillHistograms();
1043
1044 }
1045
1046}
1047
1048//________________________________________________________________________
1049void AliPWG4HighPtTrackQA::FillHistograms() {
1050
1051 fPtSel->Fill(fVariables->At(0));
1052 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1053 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1054 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1055 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1056 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1057 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1058 if(fDataType==kESD) {
1059 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1060 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1061 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
fb4a2fc7 1062 fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*fVariables->At(9));
2b553e6f 1063 fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1064 fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1065 fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1066 fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1067 fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
aa3ba8d2 1068
1069 fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1070 fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1071 fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1072 fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1073 fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1074 fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1075 fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
d756027f 1076 }
1077 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1078 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1079 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1080 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1081}
1082
1083//________________________________________________________________________
1084Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1085 //
1086 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1087 // This is to called in Notify and should provide the path to the AOD/ESD file
1088 // Copied from AliAnalysisTaskJetSpectrum2
1089 //
1090
1091 TString file(currFile);
1092 fXsec = 0;
1093 fTrials = 1;
1094
1095 if(file.Contains("root_archive.zip#")){
1096 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1097 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1098 file.Replace(pos+1,20,"");
1099 }
1100 else {
1101 // not an archive take the basename....
1102 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1103 }
1104 // Printf("%s",file.Data());
1105
1106
1107 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1108 if(!fxsec){
1109 // next trial fetch the histgram file
1110 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1111 if(!fxsec){
1112 // not a severe condition but inciate that we have no information
1113 return kFALSE;
1114 }
1115 else{
1116 // find the tlist we want to be independtent of the name so use the Tkey
1117 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1118 if(!key){
1119 fxsec->Close();
1120 return kFALSE;
1121 }
1122 TList *list = dynamic_cast<TList*>(key->ReadObj());
1123 if(!list){
1124 fxsec->Close();
1125 return kFALSE;
1126 }
1127 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1128 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1129 fxsec->Close();
1130 }
1131 } // no tree pyxsec.root
1132 else {
1133 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1134 if(!xtree){
1135 fxsec->Close();
1136 return kFALSE;
1137 }
1138 UInt_t ntrials = 0;
1139 Double_t xsection = 0;
1140 xtree->SetBranchAddress("xsection",&xsection);
1141 xtree->SetBranchAddress("ntrials",&ntrials);
1142 xtree->GetEntry(0);
1143 fTrials = ntrials;
1144 fXsec = xsection;
1145 fxsec->Close();
1146 }
1147 return kTRUE;
1148}
1149//________________________________________________________________________
1150Bool_t AliPWG4HighPtTrackQA::Notify()
1151{
1152 //
1153 // Implemented Notify() to read the cross sections
1154 // and number of trials from pyxsec.root
1155 // Copied from AliAnalysisTaskJetSpectrum2
1156 //
1157
1158 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1159 Float_t xsection = 0;
1160 Float_t ftrials = 1;
1161
1162 fAvgTrials = 1;
1163 if(tree){
1164 TFile *curfile = tree->GetCurrentFile();
1165 if (!curfile) {
1166 Error("Notify","No current file");
1167 return kFALSE;
1168 }
1169 if(!fh1Xsec||!fh1Trials){
1170 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1171 return kFALSE;
1172 }
1173 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1174 fh1Xsec->Fill("<#sigma>",xsection);
1175 // construct a poor man average trials
1176 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1177 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1178 }
1179 return kTRUE;
1180}
1181
1182//________________________________________________________________________
1183AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1184
1185 if(!mcEvent)return 0;
1186 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1187 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1188 if(!pythiaGenHeader){
1189 // cocktail ??
1190 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1191
1192 if (!genCocktailHeader) {
1193 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1194 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1195 return 0;
1196 }
1197 TList* headerList = genCocktailHeader->GetHeaders();
1198 for (Int_t i=0; i<headerList->GetEntries(); i++) {
1199 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1200 if (pythiaGenHeader)
1201 break;
1202 }
1203 if(!pythiaGenHeader){
1204 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1205 return 0;
1206 }
1207 }
1208 return pythiaGenHeader;
1209
1210}
1211
1212//_______________________________________________________________________
1213Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1214{
1215 //MV: copied from AliESDtrack since method is not available in AliAODTrack
1216
1217 //
1218 // TPC cluster information
1219 // type 0: get fraction of found/findable clusters with neighbourhood definition
1220 // 1: findable clusters with neighbourhood definition
1221 // 2: found clusters
1222 //
1223 // definition of findable clusters:
1224 // a cluster is defined as findable if there is another cluster
1225 // within +- nNeighbours pad rows. The idea is to overcome threshold
1226 // effects with a very simple algorithm.
1227 //
1228
1229 TBits fTPCClusterMap = tr->GetTPCClusterMap();
1230 if (type==2) return fTPCClusterMap.CountBits();
1231
1232 Int_t found=0;
1233 Int_t findable=0;
1234 Int_t last=-nNeighbours;
1235
1236 for (Int_t i=row0; i<row1; ++i){
1237 //look to current row
1238 if (fTPCClusterMap[i]) {
1239 last=i;
1240 ++found;
1241 ++findable;
1242 continue;
1243 }
1244 //look to nNeighbours before
1245 if ((i-last)<=nNeighbours) {
1246 ++findable;
1247 continue;
1248 }
1249 //look to nNeighbours after
1250 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1251 if (fTPCClusterMap[j]){
1252 ++findable;
1253 break;
1254 }
1255 }
1256 }
1257 if (type==1) return findable;
1258
1259 if (type==0){
1260 Float_t fraction=0;
1261 if (findable>0)
1262 fraction=(Float_t)found/(Float_t)findable;
1263 else
1264 fraction=0;
1265 return fraction;
1266 }
1267 return 0; // undefined type - default value
1268}
1269
aa3ba8d2 1270void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1271
1272 fSystTrackCuts->Fill("noCut",1);
1273 if(TMath::Abs(track->Eta())>0.9)
1274 fSystTrackCuts->Fill("eta",1);
1275 if(track->Pt()<0.15 || track->Pt()>1e10)
1276 fSystTrackCuts->Fill("0.15<pT<1e10",1);
1277 if(track->GetKinkIndex(0)>0)
1278 fSystTrackCuts->Fill("kink",1);
1279 if(track->GetTPCclusters(0)<70)
1280 fSystTrackCuts->Fill("NClusterTPC",1);
1281 if(track->GetTPCclusters(0)>0.) {
1282 if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1283 fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1284 }
1285
1286 Float_t dcaToVertexXY = 0.;
1287 Float_t dcaToVertexZ = 0.;
1288 track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1289 Float_t fCutMaxDCAToVertexXY = 2.4;
1290 Float_t fCutMaxDCAToVertexZ = 3.2;
1291 Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1292 if(dcaToVertex>1.)
1293 fSystTrackCuts->Fill("DCA2D",1);
1294
1295}
d756027f 1296
1297//________________________________________________________________________
1298void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1299{
1300 // The Terminate() function is the last function to be called during
1301 // a query. It always runs on the client, it can be used to present
1302 // the results graphically or save the results to file.
1303
1304}
1305
1306#endif