git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
trigger filling of AOD extensions
[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
1ea145ef 421 if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
422 fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 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 }
1ea145ef 452 else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
453 const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 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
1ea145ef 522
d756027f 523 Float_t cent = 999;
1ea145ef 524
525 if(esd){
526 if(esd->GetCentrality()){
527 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
528 if(fDebug>3) cout << "centrality: " << cent << endl;
529 }
d756027f 530 }
531
532 if(cent>80)return 4;
533 if(cent>50)return 3;
534 if(cent>30)return 2;
535 if(cent>10)return 1;
536 return 0;
537
538}
539
540//________________________________________________________________________
541Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
542
1ea145ef 543 if(!aod)return 4;
d756027f 544 Float_t cent = aod->GetHeader()->GetCentrality();
545 cout << "centrality: " << cent << endl;
546 if(cent>80)return 4;
547 if(cent>50)return 3;
548 if(cent>30)return 2;
549 if(cent>10)return 1;
550 return 0;
551
552}
553
554//________________________________________________________________________
555void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
556 // Main loop
557 // Called for each event
558 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
559
560 fEvent = InputEvent();
561 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
562
563 // All events without selection
564 fNEventAll->Fill(0.);
565
566 if(!SelectEvent()) {
567 // Post output data
568 PostData(1, fHistList);
569 return;
570 }
571
572
573 //Need to keep track of selected events
574 fNEventSel->Fill(0.);
575
576 fVariables = new TArrayF(fNVariables);
577
578 if(fDataType==kESD) DoAnalysisESD();
579 if(fDataType==kAOD) DoAnalysisAOD();
580
581 //Delete old fVariables
582 if(fVariables) delete fVariables;
583
584 // Post output data
585 PostData(1, fHistList);
586
587}
588
589//________________________________________________________________________
590void AliPWG4HighPtTrackQA::DoAnalysisESD() {
591
592 if(!fESD) {
593 PostData(1, fHistList);
594 return;
595 }
596
597 // ---- Get MC Header information (for MC productions in pThard bins) ----
598 Double_t ptHard = 0.;
599 Double_t nTrials = 1; // trials for MC trigger weight for real data
600
601 AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
602 if (eventHandlerMC) {
603
604 if(eventHandlerMC->MCEvent()){
605 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
606 if(pythiaGenHeader){
607 nTrials = pythiaGenHeader->Trials();
608 ptHard = pythiaGenHeader->GetPtHard();
609
610 fh1PtHard->Fill(ptHard);
611 fh1PtHardTrials->Fill(ptHard,nTrials);
612
613 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
614 }
615 }
616 }
617
618 Int_t nTracks = fESD->GetNumberOfTracks();
619 AliDebug(2,Form("nTracks ESD%d", nTracks));
620
621 /*
622 Variables to be put in fVariables
623 0: pt
624 1: phi
625 2: eta
626 3: dca2D
627 4: dcaZ
628 5: nClustersTPC
629 6: nPointITS
630 7: chi2C
631 8: nSigmaToVertex
632 9: relUncertainty1Pt
633 10: chi2PerClusterTPC
634 11: #crossed rows
635 12: (#crossed rows)/(#findable clusters)
636 */
637
638 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
639 fh1NTracksAll->Fill(0.);
640
641 //Get track for analysis
642 AliESDtrack *track;
643 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
644 if(!esdtrack) {
645 fh1NTracksReject->Fill("noESDtrack",1);
646 continue;
647 }
648
649 if(fTrackType==1)
650 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
651 else if(fTrackType==2) {
652 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
653 if(!track) {
654 fh1NTracksReject->Fill("noTPConly",1);
d756027f 655 continue;
656 }
657 AliExternalTrackParam exParam;
658 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
659 if( !relate ) {
660 fh1NTracksReject->Fill("relate",1);
661 delete track;
662 continue;
663 }
664 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
665 }
666 else
667 track = esdtrack;
668
669 if(!track) {
670 if(fTrackType==1 || fTrackType==2) delete track;
671 continue;
672 }
673
674 fPtAll->Fill(track->Pt());
675
676 if (!(fTrackCuts->AcceptTrack(track))) {
677 fh1NTracksReject->Fill("trackCuts",1);
678 if(fTrackType==1 || fTrackType==2) delete track;
679 continue;
680 }
681 if(track->GetTPCsignal()<10) { //Cut on laser tracks
682 fh1NTracksReject->Fill("laser",1);
683 if(fTrackType==1 || fTrackType==2) delete track;
684 continue;
685 }
686
687 fh1NTracksSel->Fill(0.);
688
689 fVariables->Reset(0.);
690
691 fVariables->SetAt(track->Pt(),0);
692 fVariables->SetAt(track->Phi(),1);
693 fVariables->SetAt(track->Eta(),2);
694
695 Float_t dca2D = 0.;
696 Float_t dcaz = 0.;
697 if(fTrackType==0) { //Global
698 track->GetImpactParameters(dca2D,dcaz);
699 }
700 else if(fTrackType==1 || fTrackType==2) { //TPConly
701 track->GetImpactParametersTPC(dca2D,dcaz);
702 }
703 fVariables->SetAt(dca2D,3);
704 fVariables->SetAt(dcaz,5);
705
706 fVariables->SetAt((float)track->GetTPCNcls(),5);
707
708 Int_t nPointITS = 0;
709 UChar_t itsMap = track->GetITSClusterMap();
710 for (Int_t i=0; i < 6; i++) {
711 if (itsMap & (1 << i))
712 nPointITS ++;
713 }
714 fVariables->SetAt((float)nPointITS,6);
715 fVariables->SetAt(track->GetConstrainedChi2(),7);
716 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
717
718 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
719
720 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
721
722 //cout << "#crossed rows (1): " << track->GetTPCClusterInfo(1) << endl;
723 //cout << "#crossed rows (2): " << track->GetTPCClusterInfo(2) << endl;
724 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
725 Float_t crossedRowsTPCNClsF = track->GetTPCClusterInfo(2,0);
726 //if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
727 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
728
729 FillHistograms();
730
731 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
732
733 if(fTrackType==1 || fTrackType==2) delete track;
734
735 }//track loop
736
737}
738
739//________________________________________________________________________
740void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
741
742 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1ea145ef 743 if(!aod)return;
d756027f 744 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
745
746 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
747 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
748
749 fVariables->Reset(0.);
750
751 fVariables->SetAt(aodtrack->Pt(),0);
752 fVariables->SetAt(aodtrack->Phi(),1);
753 fVariables->SetAt(aodtrack->Eta(),2);
754
755 Double_t dca[2] = {1e6,1e6};
756 Double_t covar[3] = {1e6,1e6,1e6};
757 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
758 fVariables->SetAt(dca[0],3);
759 fVariables->SetAt(dca[1],4);
760 }
761
762 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
763 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
764 fVariables->SetAt(0.,7);
765 fVariables->SetAt(0.,8);
766 fVariables->SetAt(0.,9);
767 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
768 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
769 Float_t crossedRowsTPCNClsF = 0.;
770 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
771 fVariables->SetAt(crossedRowsTPCNClsF,12);
772
773 fPtAll->Fill(fVariables->At(0));
774
775 FillHistograms();
776
777 }
778
779}
780
781//________________________________________________________________________
782void AliPWG4HighPtTrackQA::FillHistograms() {
783
784 fPtSel->Fill(fVariables->At(0));
785 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
786 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
787 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
788 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
789 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
790 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
791 if(fDataType==kESD) {
792 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
793 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
794 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
795 }
796 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
797 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
798 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
799 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
800}
801
802//________________________________________________________________________
803Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
804 //
805 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
806 // This is to called in Notify and should provide the path to the AOD/ESD file
807 // Copied from AliAnalysisTaskJetSpectrum2
808 //
809
810 TString file(currFile);
811 fXsec = 0;
812 fTrials = 1;
813
814 if(file.Contains("root_archive.zip#")){
815 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
816 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
817 file.Replace(pos+1,20,"");
818 }
819 else {
820 // not an archive take the basename....
821 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
822 }
823 // Printf("%s",file.Data());
824
825
826 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
827 if(!fxsec){
828 // next trial fetch the histgram file
829 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
830 if(!fxsec){
831 // not a severe condition but inciate that we have no information
832 return kFALSE;
833 }
834 else{
835 // find the tlist we want to be independtent of the name so use the Tkey
836 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
837 if(!key){
838 fxsec->Close();
839 return kFALSE;
840 }
841 TList *list = dynamic_cast<TList*>(key->ReadObj());
842 if(!list){
843 fxsec->Close();
844 return kFALSE;
845 }
846 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
847 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
848 fxsec->Close();
849 }
850 } // no tree pyxsec.root
851 else {
852 TTree *xtree = (TTree*)fxsec->Get("Xsection");
853 if(!xtree){
854 fxsec->Close();
855 return kFALSE;
856 }
857 UInt_t ntrials = 0;
858 Double_t xsection = 0;
859 xtree->SetBranchAddress("xsection",&xsection);
860 xtree->SetBranchAddress("ntrials",&ntrials);
861 xtree->GetEntry(0);
862 fTrials = ntrials;
863 fXsec = xsection;
864 fxsec->Close();
865 }
866 return kTRUE;
867}
868//________________________________________________________________________
869Bool_t AliPWG4HighPtTrackQA::Notify()
870{
871 //
872 // Implemented Notify() to read the cross sections
873 // and number of trials from pyxsec.root
874 // Copied from AliAnalysisTaskJetSpectrum2
875 //
876
877 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
878 Float_t xsection = 0;
879 Float_t ftrials = 1;
880
881 fAvgTrials = 1;
882 if(tree){
883 TFile *curfile = tree->GetCurrentFile();
884 if (!curfile) {
885 Error("Notify","No current file");
886 return kFALSE;
887 }
888 if(!fh1Xsec||!fh1Trials){
889 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
890 return kFALSE;
891 }
892 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
893 fh1Xsec->Fill("<#sigma>",xsection);
894 // construct a poor man average trials
895 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
896 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
897 }
898 return kTRUE;
899}
900
901//________________________________________________________________________
902AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
903
904 if(!mcEvent)return 0;
905 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
906 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
907 if(!pythiaGenHeader){
908 // cocktail ??
909 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
910
911 if (!genCocktailHeader) {
912 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
913 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
914 return 0;
915 }
916 TList* headerList = genCocktailHeader->GetHeaders();
917 for (Int_t i=0; i<headerList->GetEntries(); i++) {
918 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
919 if (pythiaGenHeader)
920 break;
921 }
922 if(!pythiaGenHeader){
923 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
924 return 0;
925 }
926 }
927 return pythiaGenHeader;
928
929}
930
931//_______________________________________________________________________
932Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
933{
934 //MV: copied from AliESDtrack since method is not available in AliAODTrack
935
936 //
937 // TPC cluster information
938 // type 0: get fraction of found/findable clusters with neighbourhood definition
939 // 1: findable clusters with neighbourhood definition
940 // 2: found clusters
941 //
942 // definition of findable clusters:
943 // a cluster is defined as findable if there is another cluster
944 // within +- nNeighbours pad rows. The idea is to overcome threshold
945 // effects with a very simple algorithm.
946 //
947
948 TBits fTPCClusterMap = tr->GetTPCClusterMap();
949 if (type==2) return fTPCClusterMap.CountBits();
950
951 Int_t found=0;
952 Int_t findable=0;
953 Int_t last=-nNeighbours;
954
955 for (Int_t i=row0; i<row1; ++i){
956 //look to current row
957 if (fTPCClusterMap[i]) {
958 last=i;
959 ++found;
960 ++findable;
961 continue;
962 }
963 //look to nNeighbours before
964 if ((i-last)<=nNeighbours) {
965 ++findable;
966 continue;
967 }
968 //look to nNeighbours after
969 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
970 if (fTPCClusterMap[j]){
971 ++findable;
972 break;
973 }
974 }
975 }
976 if (type==1) return findable;
977
978 if (type==0){
979 Float_t fraction=0;
980 if (findable>0)
981 fraction=(Float_t)found/(Float_t)findable;
982 else
983 fraction=0;
984 return fraction;
985 }
986 return 0; // undefined type - default value
987}
988
989
990//________________________________________________________________________
991void AliPWG4HighPtTrackQA::Terminate(Option_t *)
992{
993 // The Terminate() function is the last function to be called during
994 // a query. It always runs on the client, it can be used to present
995 // the results graphically or save the results to file.
996
997}
998
999#endif