]> 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
Updated QA task for constraint global tracks (M. Verweij)
[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);
884 delete track;
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);
903 delete track;
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
910 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
911
912 }
913 }
914 }
915 }
d756027f 916 else
917 track = esdtrack;
918
919 if(!track) {
afc75e88 920 // if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
d756027f 921 continue;
922 }
923
aa3ba8d2 924 if(fTrackType==2 || fTrackType==4) {
2b553e6f 925 //Cut on chi2 of constrained fit
aa3ba8d2 926 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
2b553e6f 927 fh1NTracksReject->Fill("chi2",1);
928 delete track;
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);
940 if(fTrackType==1 || fTrackType==2) delete track;
941 continue;
942 }
aa3ba8d2 943
d756027f 944 fh1NTracksSel->Fill(0.);
945
946 fVariables->Reset(0.);
947
2b553e6f 948 fVariables->SetAt(track->Pt(),0);
949 fVariables->SetAt(track->Phi(),1);
950 fVariables->SetAt(track->Eta(),2);
d756027f 951
952 Float_t dca2D = 0.;
953 Float_t dcaz = 0.;
954 if(fTrackType==0) { //Global
955 track->GetImpactParameters(dca2D,dcaz);
956 }
aa3ba8d2 957 else if(fTrackType==1 || fTrackType==2 || fTrackType==4) { //TPConly
d756027f 958 track->GetImpactParametersTPC(dca2D,dcaz);
959 }
960 fVariables->SetAt(dca2D,3);
aa3ba8d2 961 fVariables->SetAt(dcaz,5);
d756027f 962
963 fVariables->SetAt((float)track->GetTPCNcls(),5);
964
965 Int_t nPointITS = 0;
966 UChar_t itsMap = track->GetITSClusterMap();
967 for (Int_t i=0; i < 6; i++) {
968 if (itsMap & (1 << i))
969 nPointITS ++;
970 }
971 fVariables->SetAt((float)nPointITS,6);
2b553e6f 972 Float_t chi2C = (float)track->GetConstrainedChi2();
aa3ba8d2 973 if(fTrackType==1 || fTrackType==2 || fTrackType==4)
2b553e6f 974 chi2C = (float)track->GetConstrainedChi2TPC();
975 fVariables->SetAt(chi2C,7);
d756027f 976 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
977
978 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
979
980 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
981
d756027f 982 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
44684f3b 983 Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
984 if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
d756027f 985 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
2b553e6f 986 fVariables->SetAt(track->GetSigmaY2(),13);
987 fVariables->SetAt(track->GetSigmaZ2(),14);
988 fVariables->SetAt(track->GetSigmaSnp2(),15);
989 fVariables->SetAt(track->GetSigmaTgl2(),16);
990 fVariables->SetAt(track->GetSigma1Pt2(),17);
d756027f 991
992 FillHistograms();
993
994 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
995
44684f3b 996 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5) delete track;
d756027f 997
998 }//track loop
999
1000}
1001
1002//________________________________________________________________________
1003void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1004
1005 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1ea145ef 1006 if(!aod)return;
d756027f 1007 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1008
1009 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1010 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1011
1012 fVariables->Reset(0.);
1013
1014 fVariables->SetAt(aodtrack->Pt(),0);
1015 fVariables->SetAt(aodtrack->Phi(),1);
1016 fVariables->SetAt(aodtrack->Eta(),2);
1017
1018 Double_t dca[2] = {1e6,1e6};
1019 Double_t covar[3] = {1e6,1e6,1e6};
1020 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1021 fVariables->SetAt(dca[0],3);
1022 fVariables->SetAt(dca[1],4);
1023 }
1024
1025 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1026 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
aa3ba8d2 1027 fVariables->SetAt(aodtrack->Chi2perNDF(),7);
d756027f 1028 fVariables->SetAt(0.,8);
1029 fVariables->SetAt(0.,9);
1030 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1031 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
1032 Float_t crossedRowsTPCNClsF = 0.;
1033 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1034 fVariables->SetAt(crossedRowsTPCNClsF,12);
1035
1036 fPtAll->Fill(fVariables->At(0));
1037
1038 FillHistograms();
1039
1040 }
1041
1042}
1043
1044//________________________________________________________________________
1045void AliPWG4HighPtTrackQA::FillHistograms() {
1046
1047 fPtSel->Fill(fVariables->At(0));
1048 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1049 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1050 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1051 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1052 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1053 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1054 if(fDataType==kESD) {
1055 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1056 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1057 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
2b553e6f 1058 fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)/fVariables->At(9));
1059 fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1060 fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1061 fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1062 fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1063 fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
aa3ba8d2 1064
1065 fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1066 fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1067 fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1068 fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1069 fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1070 fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1071 fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
d756027f 1072 }
1073 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1074 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1075 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1076 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1077}
1078
1079//________________________________________________________________________
1080Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1081 //
1082 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1083 // This is to called in Notify and should provide the path to the AOD/ESD file
1084 // Copied from AliAnalysisTaskJetSpectrum2
1085 //
1086
1087 TString file(currFile);
1088 fXsec = 0;
1089 fTrials = 1;
1090
1091 if(file.Contains("root_archive.zip#")){
1092 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1093 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1094 file.Replace(pos+1,20,"");
1095 }
1096 else {
1097 // not an archive take the basename....
1098 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1099 }
1100 // Printf("%s",file.Data());
1101
1102
1103 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
1104 if(!fxsec){
1105 // next trial fetch the histgram file
1106 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1107 if(!fxsec){
1108 // not a severe condition but inciate that we have no information
1109 return kFALSE;
1110 }
1111 else{
1112 // find the tlist we want to be independtent of the name so use the Tkey
1113 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1114 if(!key){
1115 fxsec->Close();
1116 return kFALSE;
1117 }
1118 TList *list = dynamic_cast<TList*>(key->ReadObj());
1119 if(!list){
1120 fxsec->Close();
1121 return kFALSE;
1122 }
1123 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1124 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1125 fxsec->Close();
1126 }
1127 } // no tree pyxsec.root
1128 else {
1129 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1130 if(!xtree){
1131 fxsec->Close();
1132 return kFALSE;
1133 }
1134 UInt_t ntrials = 0;
1135 Double_t xsection = 0;
1136 xtree->SetBranchAddress("xsection",&xsection);
1137 xtree->SetBranchAddress("ntrials",&ntrials);
1138 xtree->GetEntry(0);
1139 fTrials = ntrials;
1140 fXsec = xsection;
1141 fxsec->Close();
1142 }
1143 return kTRUE;
1144}
1145//________________________________________________________________________
1146Bool_t AliPWG4HighPtTrackQA::Notify()
1147{
1148 //
1149 // Implemented Notify() to read the cross sections
1150 // and number of trials from pyxsec.root
1151 // Copied from AliAnalysisTaskJetSpectrum2
1152 //
1153
1154 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1155 Float_t xsection = 0;
1156 Float_t ftrials = 1;
1157
1158 fAvgTrials = 1;
1159 if(tree){
1160 TFile *curfile = tree->GetCurrentFile();
1161 if (!curfile) {
1162 Error("Notify","No current file");
1163 return kFALSE;
1164 }
1165 if(!fh1Xsec||!fh1Trials){
1166 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1167 return kFALSE;
1168 }
1169 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1170 fh1Xsec->Fill("<#sigma>",xsection);
1171 // construct a poor man average trials
1172 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1173 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1174 }
1175 return kTRUE;
1176}
1177
1178//________________________________________________________________________
1179AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1180
1181 if(!mcEvent)return 0;
1182 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1183 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1184 if(!pythiaGenHeader){
1185 // cocktail ??
1186 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1187
1188 if (!genCocktailHeader) {
1189 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1190 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1191 return 0;
1192 }
1193 TList* headerList = genCocktailHeader->GetHeaders();
1194 for (Int_t i=0; i<headerList->GetEntries(); i++) {
1195 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1196 if (pythiaGenHeader)
1197 break;
1198 }
1199 if(!pythiaGenHeader){
1200 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1201 return 0;
1202 }
1203 }
1204 return pythiaGenHeader;
1205
1206}
1207
1208//_______________________________________________________________________
1209Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1210{
1211 //MV: copied from AliESDtrack since method is not available in AliAODTrack
1212
1213 //
1214 // TPC cluster information
1215 // type 0: get fraction of found/findable clusters with neighbourhood definition
1216 // 1: findable clusters with neighbourhood definition
1217 // 2: found clusters
1218 //
1219 // definition of findable clusters:
1220 // a cluster is defined as findable if there is another cluster
1221 // within +- nNeighbours pad rows. The idea is to overcome threshold
1222 // effects with a very simple algorithm.
1223 //
1224
1225 TBits fTPCClusterMap = tr->GetTPCClusterMap();
1226 if (type==2) return fTPCClusterMap.CountBits();
1227
1228 Int_t found=0;
1229 Int_t findable=0;
1230 Int_t last=-nNeighbours;
1231
1232 for (Int_t i=row0; i<row1; ++i){
1233 //look to current row
1234 if (fTPCClusterMap[i]) {
1235 last=i;
1236 ++found;
1237 ++findable;
1238 continue;
1239 }
1240 //look to nNeighbours before
1241 if ((i-last)<=nNeighbours) {
1242 ++findable;
1243 continue;
1244 }
1245 //look to nNeighbours after
1246 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1247 if (fTPCClusterMap[j]){
1248 ++findable;
1249 break;
1250 }
1251 }
1252 }
1253 if (type==1) return findable;
1254
1255 if (type==0){
1256 Float_t fraction=0;
1257 if (findable>0)
1258 fraction=(Float_t)found/(Float_t)findable;
1259 else
1260 fraction=0;
1261 return fraction;
1262 }
1263 return 0; // undefined type - default value
1264}
1265
aa3ba8d2 1266void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1267
1268 fSystTrackCuts->Fill("noCut",1);
1269 if(TMath::Abs(track->Eta())>0.9)
1270 fSystTrackCuts->Fill("eta",1);
1271 if(track->Pt()<0.15 || track->Pt()>1e10)
1272 fSystTrackCuts->Fill("0.15<pT<1e10",1);
1273 if(track->GetKinkIndex(0)>0)
1274 fSystTrackCuts->Fill("kink",1);
1275 if(track->GetTPCclusters(0)<70)
1276 fSystTrackCuts->Fill("NClusterTPC",1);
1277 if(track->GetTPCclusters(0)>0.) {
1278 if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1279 fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1280 }
1281
1282 Float_t dcaToVertexXY = 0.;
1283 Float_t dcaToVertexZ = 0.;
1284 track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1285 Float_t fCutMaxDCAToVertexXY = 2.4;
1286 Float_t fCutMaxDCAToVertexZ = 3.2;
1287 Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1288 if(dcaToVertex>1.)
1289 fSystTrackCuts->Fill("DCA2D",1);
1290
1291}
d756027f 1292
1293//________________________________________________________________________
1294void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1295{
1296 // The Terminate() function is the last function to be called during
1297 // a query. It always runs on the client, it can be used to present
1298 // the results graphically or save the results to file.
1299
1300}
1301
1302#endif