git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding new track QA task running on ESD and AOD, bug fixes to QAMC (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),
67 fTrackCuts(0),
68 fTrackType(0),
69 fFilterMask(0),
70 fPtMax(100.),
71 fIsPbPb(0),
72 fCentClass(10),
73 fNVariables(13),
74 fVariables(0x0),
75 fAvgTrials(1),
76 fNEventAll(0),
77 fNEventSel(0),
78 fNEventReject(0),
79 fh1Centrality(0x0),
80 fh1Xsec(0),
81 fh1Trials(0),
82 fh1PtHard(0),
83 fh1PtHardTrials(0),
84 fh1NTracksAll(0x0),
85 fh1NTracksReject(0x0),
86 fh1NTracksSel(0x0),
87 fPtAll(0),
88 fPtSel(0),
89 fPtPhi(0x0),
90 fPtEta(0x0),
91 fPtDCA2D(0x0),
92 fPtDCAZ(0x0),
93 fPtNClustersTPC(0x0),
94 fPtNPointITS(0x0),
95 fPtChi2C(0x0),
96 fPtNSigmaToVertex(0x0),
97 fPtRelUncertainty1Pt(0x0),
98 fPtChi2PerClusterTPC(0x0),
99 fPtNCrossedRows(0x0),
100 fPtNCrossedRowsNClusF(0x0),
101 fPtNCrRNCrRNClusF(0x0),
102 fHistList(0)
103{
104 SetNVariables(13);
105}
106//________________________________________________________________________
107AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name):
108 AliAnalysisTaskSE(name),
109 fDataType(kESD),
110 fEvent(0x0),
111 fESD(0x0),
112 fVtx(0x0),
113 fTrackCuts(),
114 fTrackType(0),
115 fFilterMask(0),
116 fPtMax(100.),
117 fIsPbPb(0),
118 fCentClass(10),
119 fNVariables(13),
120 fVariables(0x0),
121 fAvgTrials(1),
122 fNEventAll(0),
123 fNEventSel(0),
124 fNEventReject(0),
125 fh1Centrality(0x0),
126 fh1Xsec(0),
127 fh1Trials(0),
128 fh1PtHard(0),
129 fh1PtHardTrials(0),
130 fh1NTracksAll(0x0),
131 fh1NTracksReject(0x0),
132 fh1NTracksSel(0x0),
133 fPtAll(0),
134 fPtSel(0),
135 fPtPhi(0x0),
136 fPtEta(0x0),
137 fPtDCA2D(0x0),
138 fPtDCAZ(0x0),
139 fPtNClustersTPC(0x0),
140 fPtNPointITS(0x0),
141 fPtChi2C(0x0),
142 fPtNSigmaToVertex(0x0),
143 fPtRelUncertainty1Pt(0x0),
144 fPtChi2PerClusterTPC(0x0),
145 fPtNCrossedRows(0x0),
146 fPtNCrossedRowsNClusF(0x0),
147 fPtNCrRNCrRNClusF(0x0),
148 fHistList(0)
149{
150 //
151 // Constructor. Initialization of Inputs and Outputs
152 //
153 AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
154
155 SetNVariables(13);
156
157 // Input slot #0 works with a TChain ESD
158 DefineInput(0, TChain::Class());
159 // Output slot #1 write into a TList
160 DefineOutput(1, TList::Class());
161}
162
163//________________________________________________________________________
164void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
165 //Create output objects
166 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
167
168 Bool_t oldStatus = TH1::AddDirectoryStatus();
169 TH1::AddDirectory(kFALSE);
170
171 OpenFile(1);
172 fHistList = new TList();
173 fHistList->SetOwner(kTRUE);
174
175 Float_t fgkPtMin = 0.;
176 Float_t fgkPtMax = fPtMax;
177
178 Float_t ptBinEdges[2][2];
179 ptBinEdges[0][0] = 10.;
180 ptBinEdges[0][1] = 1.;
181 ptBinEdges[1][0] = 20.;
182 ptBinEdges[1][1] = 2.;
183 Float_t binWidth3 = 5.;
184 if(fPtMax>100.) {
185 ptBinEdges[0][0] = 100.;
186 ptBinEdges[0][1] = 5.;
187 ptBinEdges[1][0] = 300.;
188 ptBinEdges[1][1] = 10.;
189 binWidth3 = 20.;
190 }
191
192 const Float_t ptmin1 = fgkPtMin;
193 const Float_t ptmax1 = ptBinEdges[0][0];
194 const Float_t ptmin2 = ptmax1 ;
195 const Float_t ptmax2 = ptBinEdges[1][0];
196 const Float_t ptmin3 = ptmax2 ;
197 const Float_t ptmax3 = fgkPtMax;
198 const Int_t nbin11 = (int)((ptmax1-ptmin1)/ptBinEdges[0][1]);
199 const Int_t nbin12 = (int)((ptmax2-ptmin2)/ptBinEdges[1][1])+nbin11;
200 const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
201 Int_t fgkNPtBins=nbin13;
202 //Create array with low edges of each bin
203 Double_t *binsPt=new Double_t[fgkNPtBins+1];
204 for(Int_t i=0; i<=fgkNPtBins; i++) {
205 if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
206 if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
207 if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
208 }
209
210 Int_t fgkNPhiBins = 18*6;
211 Float_t kMinPhi = 0.;
212 Float_t kMaxPhi = 2.*TMath::Pi();
213 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
214 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
215
216 Int_t fgkNEtaBins=20;
217 Float_t fgkEtaMin = -1.;
218 Float_t fgkEtaMax = 1.;
219 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
220 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
221
222 Int_t fgkNNClustersTPCBins=80;
223 Float_t fgkNClustersTPCMin = 0.5;
224 Float_t fgkNClustersTPCMax = 160.5;
225 Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
226 for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
227
228 Int_t fgkNDCA2DBins=80;
229 Float_t fgkDCA2DMin = -0.2;
230 Float_t fgkDCA2DMax = 0.2;
231 Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
232 for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
233
234 Int_t fgkNDCAZBins=80;
235 Float_t fgkDCAZMin = -2.;
236 Float_t fgkDCAZMax = 2.;
237 Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
238 for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
239
240 Int_t fgkNNPointITSBins=9;
241 Float_t fgkNPointITSMin = -0.5;
242 Float_t fgkNPointITSMax = 8.5;
243 Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
244 for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
245
246 Int_t fgkNNSigmaToVertexBins=40;
247 Float_t fgkNSigmaToVertexMin = 0.;
248 Float_t fgkNSigmaToVertexMax = 8.;
249 Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
250 for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
251
252 Int_t fgkNChi2CBins=20;
253 Float_t fgkChi2CMin = 0.;
254 Float_t fgkChi2CMax = 10.;
255 Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
256 for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
257
258 Int_t fgkNRel1PtUncertaintyBins=30;
259 Float_t fgkRel1PtUncertaintyMin = 0.;
260 Float_t fgkRel1PtUncertaintyMax = 0.3;
261 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
262 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
263
264 Float_t fgkChi2PerClusMin = 0.;
265 Float_t fgkChi2PerClusMax = 4.;
266 Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
267 Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
268 for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
269
270 Int_t fgkNCrossedRowsNClusFBins = 50;
271 Float_t fgkNCrossedRowsNClusFMin = 0.;
272 Float_t fgkNCrossedRowsNClusFMax = 1.;
273 Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
274 for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
275
276 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
277 fHistList->Add(fNEventAll);
278 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
279 fHistList->Add(fNEventSel);
280 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
281 //Set labels
282 fNEventReject->Fill("noESD",0);
283 fNEventReject->Fill("Trigger",0);
284 fNEventReject->Fill("NTracks<2",0);
285 fNEventReject->Fill("noVTX",0);
286 fNEventReject->Fill("VtxStatus",0);
287 fNEventReject->Fill("NCont<2",0);
288 fNEventReject->Fill("ZVTX>10",0);
289 fNEventReject->Fill("cent",0);
290 fNEventReject->Fill("cent>90",0);
291 fHistList->Add(fNEventReject);
292
293 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
294 fHistList->Add(fh1Centrality);
295
296 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
297 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
298 fHistList->Add(fh1Xsec);
299
300 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
301 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
302 fHistList->Add(fh1Trials);
303
304 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
305 fHistList->Add(fh1PtHard);
306 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
307 fHistList->Add(fh1PtHardTrials);
308
309 fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
310 fHistList->Add(fh1NTracksAll);
311
312 fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
313 fh1NTracksReject->Fill("noESDtrack",0);
314 fh1NTracksReject->Fill("noTPCInner",0);
315 fh1NTracksReject->Fill("FillTPC",0);
316 fh1NTracksReject->Fill("noTPConly",0);
317 fh1NTracksReject->Fill("relate",0);
318 fh1NTracksReject->Fill("trackCuts",0);
319 fh1NTracksReject->Fill("laser",0);
320 fHistList->Add(fh1NTracksReject);
321
322 fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
323 fHistList->Add(fh1NTracksSel);
324
325
326 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
327 fHistList->Add(fPtAll);
328 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
329 fHistList->Add(fPtSel);
330
331 fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
332 fHistList->Add(fPtPhi);
333
334 fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
335 fHistList->Add(fPtEta);
336
337 fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
338 fHistList->Add(fPtDCA2D);
339
340 fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
341 fHistList->Add(fPtDCAZ);
342
343 fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
344 fHistList->Add(fPtNClustersTPC);
345
346 fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
347 fHistList->Add(fPtNPointITS);
348
349 fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
350 fHistList->Add(fPtChi2C);
351
352 fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
353 fHistList->Add(fPtNSigmaToVertex);
354
355 fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
356 fHistList->Add(fPtRelUncertainty1Pt);
357
358 fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
359 fHistList->Add(fPtChi2PerClusterTPC);
360
361 fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
362 fHistList->Add(fPtNCrossedRows);
363
364 fPtNCrossedRowsNClusF = new TH2F("fPtNCrossedRowsNClusF","fPtNCrossedRowsNClusF",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
365 fHistList->Add(fPtNCrossedRowsNClusF);
366
367 fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
368 fHistList->Add(fPtNCrRNCrRNClusF);
369
370 TH1::AddDirectory(oldStatus);
371
372 PostData(1, fHistList);
373
374 if(binsPhi) delete [] binsPhi;
375 if(binsPt) delete [] binsPt;
376 if(binsNClustersTPC) delete [] binsNClustersTPC;
377 if(binsDCA2D) delete [] binsDCA2D;
378 if(binsDCAZ) delete [] binsDCAZ;
379 if(binsNPointITS) delete [] binsNPointITS;
380 if(binsNSigmaToVertex) delete [] binsNSigmaToVertex;
381 if(binsChi2C) delete [] binsChi2C;
382 if(binsEta) delete [] binsEta;
383 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
384 if(binsChi2PerClus) delete [] binsChi2PerClus;
385 if(binsChi2PerClus) delete [] binsNCrossedRowsNClusF;
386}
387
388//________________________________________________________________________
389Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
390 //
391 // Decide if event should be selected for analysis
392 //
393
394 // Checks following requirements:
395 // - fEvent available
396 // - trigger info from AliPhysicsSelection
397 // - MCevent available
398 // - number of reconstructed tracks > 1
399 // - primary vertex reconstructed
400 // - z-vertex < 10 cm
401 // - centrality in case of PbPb
402
403 Bool_t selectEvent = kTRUE;
404
405 //fEvent object available?
406 if (!fEvent) {
407 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
408 fNEventReject->Fill("noAliVEvent",1);
409 selectEvent = kFALSE;
410 return selectEvent;
411 }
412
413 //Check if number of reconstructed tracks is larger than 1
414 if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2) {
415 fNEventReject->Fill("NTracks<2",1);
416 selectEvent = kFALSE;
417 return selectEvent;
418 }
419
420 //Check if vertex is reconstructed
421 if(fDataType==kESD) {
422 fVtx = dynamic_cast<AliESDEvent*>(fEvent)->GetPrimaryVertexSPD();
423
424 if(!fVtx) {
425 fNEventReject->Fill("noVTX",1);
426 selectEvent = kFALSE;
427 return selectEvent;
428 }
429
430 if(!fVtx->GetStatus()) {
431 fNEventReject->Fill("VtxStatus",1);
432 selectEvent = kFALSE;
433 return selectEvent;
434 }
435
436 // Need vertex cut
437 if(fVtx->GetNContributors()<2) {
438 fNEventReject->Fill("NCont<2",1);
439 selectEvent = kFALSE;
440 return selectEvent;
441 }
442
443 //Check if z-vertex < 10 cm
444 double primVtx[3];
445 fVtx->GetXYZ(primVtx);
446 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
447 fNEventReject->Fill("ZVTX>10",1);
448 selectEvent = kFALSE;
449 return selectEvent;
450 }
451 }
452 else if(fDataType==kAOD) {
453 const AliAODVertex *vtx = dynamic_cast<AliAODEvent*>(fEvent)->GetPrimaryVertexSPD();
454 if(!vtx) {
455 fNEventReject->Fill("noVTX",1);
456 selectEvent = kFALSE;
457 return selectEvent;
458 }
459
460 // Need vertex cut
461 if(vtx->GetNContributors()<2) {
462 fNEventReject->Fill("NCont<2",1);
463 selectEvent = kFALSE;
464 return selectEvent;
465 }
466
467 //Check if z-vertex < 10 cm
468 double primVtx[3];
469 vtx->GetXYZ(primVtx);
470 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
471 fNEventReject->Fill("ZVTX>10",1);
472 selectEvent = kFALSE;
473 return selectEvent;
474 }
475
476 }
477
478 //Centrality selection should only be done in case of PbPb
479 if(IsPbPb()) {
480 Float_t cent = 0.;
481 if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
482 fNEventReject->Fill("cent",1);
483 selectEvent = kFALSE;
484 return selectEvent;
485 }
486 else {
487 if(fDataType==kESD) {
488 if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
489 cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
490 }
491 }
492 else if(fDataType==kAOD) {
493 if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
494 cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
495 }
496 if(cent>90.) {
497 fNEventReject->Fill("cent>90",1);
498 selectEvent = kFALSE;
499 return selectEvent;
500 }
501 fh1Centrality->Fill(cent);
502 }
503 }
504
505 return selectEvent;
506
507}
508
509//________________________________________________________________________
510Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
511 if(fDataType==kESD)
512 return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
513 else if(fDataType==kAOD)
514 return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
515 else
516 return 5;
517}
518
519//________________________________________________________________________
520Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
521
522 Float_t cent = 999;
523 if(esd->GetCentrality()){
524 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
525 if(fDebug>3) cout << "centrality: " << cent << endl;
526 }
527
528 if(cent>80)return 4;
529 if(cent>50)return 3;
530 if(cent>30)return 2;
531 if(cent>10)return 1;
532 return 0;
533
534}
535
536//________________________________________________________________________
537Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
538
539 Float_t cent = aod->GetHeader()->GetCentrality();
540 cout << "centrality: " << cent << endl;
541 if(cent>80)return 4;
542 if(cent>50)return 3;
543 if(cent>30)return 2;
544 if(cent>10)return 1;
545 return 0;
546
547}
548
549//________________________________________________________________________
550void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
551 // Main loop
552 // Called for each event
553 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
554
555 fEvent = InputEvent();
556 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
557
558 // All events without selection
559 fNEventAll->Fill(0.);
560
561 if(!SelectEvent()) {
562 // Post output data
563 PostData(1, fHistList);
564 return;
565 }
566
567
568 //Need to keep track of selected events
569 fNEventSel->Fill(0.);
570
571 fVariables = new TArrayF(fNVariables);
572
573 if(fDataType==kESD) DoAnalysisESD();
574 if(fDataType==kAOD) DoAnalysisAOD();
575
576 //Delete old fVariables
577 if(fVariables) delete fVariables;
578
579 // Post output data
580 PostData(1, fHistList);
581
582}
583
584//________________________________________________________________________
585void AliPWG4HighPtTrackQA::DoAnalysisESD() {
586
587 if(!fESD) {
588 PostData(1, fHistList);
589 return;
590 }
591
592 // ---- Get MC Header information (for MC productions in pThard bins) ----
593 Double_t ptHard = 0.;
594 Double_t nTrials = 1; // trials for MC trigger weight for real data
595
596 AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
597 if (eventHandlerMC) {
598
599 if(eventHandlerMC->MCEvent()){
600 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
601 if(pythiaGenHeader){
602 nTrials = pythiaGenHeader->Trials();
603 ptHard = pythiaGenHeader->GetPtHard();
604
605 fh1PtHard->Fill(ptHard);
606 fh1PtHardTrials->Fill(ptHard,nTrials);
607
608 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
609 }
610 }
611 }
612
613 Int_t nTracks = fESD->GetNumberOfTracks();
614 AliDebug(2,Form("nTracks ESD%d", nTracks));
615
616 /*
617 Variables to be put in fVariables
618 0: pt
619 1: phi
620 2: eta
621 3: dca2D
622 4: dcaZ
623 5: nClustersTPC
624 6: nPointITS
625 7: chi2C
626 8: nSigmaToVertex
627 9: relUncertainty1Pt
628 10: chi2PerClusterTPC
629 11: #crossed rows
630 12: (#crossed rows)/(#findable clusters)
631 */
632
633 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
634 fh1NTracksAll->Fill(0.);
635
636 //Get track for analysis
637 AliESDtrack *track;
638 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
639 if(!esdtrack) {
640 fh1NTracksReject->Fill("noESDtrack",1);
641 continue;
642 }
643
644 if(fTrackType==1)
645 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
646 else if(fTrackType==2) {
647 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
648 if(!track) {
649 fh1NTracksReject->Fill("noTPConly",1);
650 delete track;
651 continue;
652 }
653 AliExternalTrackParam exParam;
654 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
655 if( !relate ) {
656 fh1NTracksReject->Fill("relate",1);
657 delete track;
658 continue;
659 }
660 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
661 }
662 else
663 track = esdtrack;
664
665 if(!track) {
666 if(fTrackType==1 || fTrackType==2) delete track;
667 continue;
668 }
669
670 fPtAll->Fill(track->Pt());
671
672 if (!(fTrackCuts->AcceptTrack(track))) {
673 fh1NTracksReject->Fill("trackCuts",1);
674 if(fTrackType==1 || fTrackType==2) delete track;
675 continue;
676 }
677 if(track->GetTPCsignal()<10) { //Cut on laser tracks
678 fh1NTracksReject->Fill("laser",1);
679 if(fTrackType==1 || fTrackType==2) delete track;
680 continue;
681 }
682
683 fh1NTracksSel->Fill(0.);
684
685 fVariables->Reset(0.);
686
687 fVariables->SetAt(track->Pt(),0);
688 fVariables->SetAt(track->Phi(),1);
689 fVariables->SetAt(track->Eta(),2);
690
691 Float_t dca2D = 0.;
692 Float_t dcaz = 0.;
693 if(fTrackType==0) { //Global
694 track->GetImpactParameters(dca2D,dcaz);
695 }
696 else if(fTrackType==1 || fTrackType==2) { //TPConly
697 track->GetImpactParametersTPC(dca2D,dcaz);
698 }
699 fVariables->SetAt(dca2D,3);
700 fVariables->SetAt(dcaz,5);
701
702 fVariables->SetAt((float)track->GetTPCNcls(),5);
703
704 Int_t nPointITS = 0;
705 UChar_t itsMap = track->GetITSClusterMap();
706 for (Int_t i=0; i < 6; i++) {
707 if (itsMap & (1 << i))
708 nPointITS ++;
709 }
710 fVariables->SetAt((float)nPointITS,6);
711 fVariables->SetAt(track->GetConstrainedChi2(),7);
712 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
713
714 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
715
716 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
717
718 //cout << "#crossed rows (1): " << track->GetTPCClusterInfo(1) << endl;
719 //cout << "#crossed rows (2): " << track->GetTPCClusterInfo(2) << endl;
720 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
721 Float_t crossedRowsTPCNClsF = track->GetTPCClusterInfo(2,0);
722 //if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
723 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
724
725 FillHistograms();
726
727 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
728
729 if(fTrackType==1 || fTrackType==2) delete track;
730
731 }//track loop
732
733}
734
735//________________________________________________________________________
736void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
737
738 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
739 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
740
741 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
742 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
743
744 fVariables->Reset(0.);
745
746 fVariables->SetAt(aodtrack->Pt(),0);
747 fVariables->SetAt(aodtrack->Phi(),1);
748 fVariables->SetAt(aodtrack->Eta(),2);
749
750 Double_t dca[2] = {1e6,1e6};
751 Double_t covar[3] = {1e6,1e6,1e6};
752 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
753 fVariables->SetAt(dca[0],3);
754 fVariables->SetAt(dca[1],4);
755 }
756
757 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
758 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
759 fVariables->SetAt(0.,7);
760 fVariables->SetAt(0.,8);
761 fVariables->SetAt(0.,9);
762 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
763 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
764 Float_t crossedRowsTPCNClsF = 0.;
765 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
766 fVariables->SetAt(crossedRowsTPCNClsF,12);
767
768 fPtAll->Fill(fVariables->At(0));
769
770 FillHistograms();
771
772 }
773
774}
775
776//________________________________________________________________________
777void AliPWG4HighPtTrackQA::FillHistograms() {
778
779 fPtSel->Fill(fVariables->At(0));
780 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
781 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
782 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
783 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
784 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
785 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
786 if(fDataType==kESD) {
787 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
788 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
789 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
790 }
791 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
792 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
793 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
794 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
795}
796
797//________________________________________________________________________
798Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
799 //
800 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
801 // This is to called in Notify and should provide the path to the AOD/ESD file
802 // Copied from AliAnalysisTaskJetSpectrum2
803 //
804
805 TString file(currFile);
806 fXsec = 0;
807 fTrials = 1;
808
809 if(file.Contains("root_archive.zip#")){
810 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
811 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
812 file.Replace(pos+1,20,"");
813 }
814 else {
815 // not an archive take the basename....
816 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
817 }
818 // Printf("%s",file.Data());
819
820
821 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
822 if(!fxsec){
823 // next trial fetch the histgram file
824 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
825 if(!fxsec){
826 // not a severe condition but inciate that we have no information
827 return kFALSE;
828 }
829 else{
830 // find the tlist we want to be independtent of the name so use the Tkey
831 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
832 if(!key){
833 fxsec->Close();
834 return kFALSE;
835 }
836 TList *list = dynamic_cast<TList*>(key->ReadObj());
837 if(!list){
838 fxsec->Close();
839 return kFALSE;
840 }
841 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
842 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
843 fxsec->Close();
844 }
845 } // no tree pyxsec.root
846 else {
847 TTree *xtree = (TTree*)fxsec->Get("Xsection");
848 if(!xtree){
849 fxsec->Close();
850 return kFALSE;
851 }
852 UInt_t ntrials = 0;
853 Double_t xsection = 0;
854 xtree->SetBranchAddress("xsection",&xsection);
855 xtree->SetBranchAddress("ntrials",&ntrials);
856 xtree->GetEntry(0);
857 fTrials = ntrials;
858 fXsec = xsection;
859 fxsec->Close();
860 }
861 return kTRUE;
862}
863//________________________________________________________________________
864Bool_t AliPWG4HighPtTrackQA::Notify()
865{
866 //
867 // Implemented Notify() to read the cross sections
868 // and number of trials from pyxsec.root
869 // Copied from AliAnalysisTaskJetSpectrum2
870 //
871
872 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
873 Float_t xsection = 0;
874 Float_t ftrials = 1;
875
876 fAvgTrials = 1;
877 if(tree){
878 TFile *curfile = tree->GetCurrentFile();
879 if (!curfile) {
880 Error("Notify","No current file");
881 return kFALSE;
882 }
883 if(!fh1Xsec||!fh1Trials){
884 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
885 return kFALSE;
886 }
887 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
888 fh1Xsec->Fill("<#sigma>",xsection);
889 // construct a poor man average trials
890 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
891 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
892 }
893 return kTRUE;
894}
895
896//________________________________________________________________________
897AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
898
899 if(!mcEvent)return 0;
900 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
901 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
902 if(!pythiaGenHeader){
903 // cocktail ??
904 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
905
906 if (!genCocktailHeader) {
907 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
908 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
909 return 0;
910 }
911 TList* headerList = genCocktailHeader->GetHeaders();
912 for (Int_t i=0; i<headerList->GetEntries(); i++) {
913 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
914 if (pythiaGenHeader)
915 break;
916 }
917 if(!pythiaGenHeader){
918 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
919 return 0;
920 }
921 }
922 return pythiaGenHeader;
923
924}
925
926//_______________________________________________________________________
927Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
928{
929 //MV: copied from AliESDtrack since method is not available in AliAODTrack
930
931 //
932 // TPC cluster information
933 // type 0: get fraction of found/findable clusters with neighbourhood definition
934 // 1: findable clusters with neighbourhood definition
935 // 2: found clusters
936 //
937 // definition of findable clusters:
938 // a cluster is defined as findable if there is another cluster
939 // within +- nNeighbours pad rows. The idea is to overcome threshold
940 // effects with a very simple algorithm.
941 //
942
943 TBits fTPCClusterMap = tr->GetTPCClusterMap();
944 if (type==2) return fTPCClusterMap.CountBits();
945
946 Int_t found=0;
947 Int_t findable=0;
948 Int_t last=-nNeighbours;
949
950 for (Int_t i=row0; i<row1; ++i){
951 //look to current row
952 if (fTPCClusterMap[i]) {
953 last=i;
954 ++found;
955 ++findable;
956 continue;
957 }
958 //look to nNeighbours before
959 if ((i-last)<=nNeighbours) {
960 ++findable;
961 continue;
962 }
963 //look to nNeighbours after
964 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
965 if (fTPCClusterMap[j]){
966 ++findable;
967 break;
968 }
969 }
970 }
971 if (type==1) return findable;
972
973 if (type==0){
974 Float_t fraction=0;
975 if (findable>0)
976 fraction=(Float_t)found/(Float_t)findable;
977 else
978 fraction=0;
979 return fraction;
980 }
981 return 0; // undefined type - default value
982}
983
984
985//________________________________________________________________________
986void AliPWG4HighPtTrackQA::Terminate(Option_t *)
987{
988 // The Terminate() function is the last function to be called during
989 // a query. It always runs on the client, it can be used to present
990 // the results graphically or save the results to file.
991
992}
993
994#endif