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