227ba57c63cd000b0f623ad4d72a0be724c6029f
[u/mrichter/AliRoot.git] / PWGJE / 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(25),
77   fVariables(0x0),
78   fITSClusterMap(0),
79   fAvgTrials(1),
80   fNEventAll(0),
81   fNEventSel(0),
82   fNEventReject(0),
83   fh1Centrality(0x0),
84   fh1Xsec(0),
85   fh1Trials(0),
86   fh1PtHard(0),
87   fh1PtHardTrials(0),
88   fh1NTracksAll(0x0),
89   fh1NTracksReject(0x0),
90   fh1NTracksSel(0x0),
91   fPtAll(0),  
92   fPtSel(0),    
93   fPtPhi(0x0),
94   fPtEta(0x0),
95   fPtDCA2D(0x0),
96   fPtDCAZ(0x0),
97   fPtNClustersTPC(0x0),
98   fPtNClustersTPCIter1(0x0),
99   fPtNClustersTPCIter1Phi(0x0),
100   fPtNClustersTPCShared(0x0),
101   fPtNClustersTPCSharedFrac(0x0),
102   fPtNPointITS(0x0),
103   fPtChi2C(0x0),
104   fPtNSigmaToVertex(0x0),
105   fPtRelUncertainty1Pt(0x0),
106   fPtRelUncertainty1PtNClus(0x0),
107   fPtRelUncertainty1PtNClusIter1(0x0),
108   fPtRelUncertainty1PtNPointITS(0x0),
109   fPtRelUncertainty1PtITSClusterMap(0x0),
110   fPtRelUncertainty1PtChi2(0x0),
111   fPtRelUncertainty1PtChi2Iter1(0x0),
112   fPtRelUncertainty1PtPhi(0x0),
113   fPtUncertainty1Pt(0x0),
114   fPtChi2PerClusterTPC(0x0),
115   fPtChi2PerClusterTPCIter1(0x0),
116   fPtNCrossedRows(0x0),
117   fPtNCrossedRowsPhi(0x0),
118   fPtNCrossedRowsNClusFPhi(0x0),
119   fPtNCrRNCrRNClusF(0x0),
120   fPtNCrossedRowsFit(0x0),
121   fPtNCrossedRowsFitPhi(0x0),
122   fPtNCrossedRowsNClusFFitPhi(0x0),
123   fNCrossedRowsNCrossedRowsFit(0x0),
124   fNClustersNCrossedRows(0x0),
125   fNClustersNCrossedRowsFit(0x0),
126   fPtRelUncertainty1PtNCrossedRows(0x0),
127   fPtRelUncertainty1PtNCrossedRowsFit(0x0),
128   fPtChi2Gold(0x0),
129   fPtChi2GGC(0x0),
130   fPtChi2GoldPhi(0x0),
131   fPtChi2GGCPhi(0x0),
132   fChi2GoldChi2GGC(0x0),
133   fPtSigmaY2(0x0),
134   fPtSigmaZ2(0x0),
135   fPtSigmaSnp2(0x0),
136   fPtSigmaTgl2(0x0),
137   fPtSigma1Pt2(0x0),
138   fProfPtSigmaY2(0x0),
139   fProfPtSigmaZ2(0x0),
140   fProfPtSigmaSnp2(0x0),
141   fProfPtSigmaTgl2(0x0),
142   fProfPtSigma1Pt2(0x0),
143   fProfPtSigma1Pt(0x0),
144   fProfPtPtSigma1Pt(0x0),
145   fHistList(0)
146 {
147   SetNVariables(25);
148
149   fPtBinEdges[0][0] = 10.;
150   fPtBinEdges[0][1] = 1.;
151   fPtBinEdges[1][0] = 20.;
152   fPtBinEdges[1][1] = 2.;
153   fPtBinEdges[2][0] = 100.;
154   fPtBinEdges[2][1] = 5.;
155
156 }
157 //________________________________________________________________________
158 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name): 
159   AliAnalysisTaskSE(name), 
160   fDataType(kESD),
161   fEvent(0x0),
162   fESD(0x0),
163   fVtx(0x0),
164   fTrackCuts(0x0),
165   fTrackCutsITSLoose(0x0), 
166   fTrackCutsTPConly(0x0), 
167   fTrackType(0),
168   fFilterMask(0),
169   fSigmaConstrainedMax(-1.),
170   fPtMax(100.),
171   fIsPbPb(0),
172   fCentClass(10),
173   fNVariables(25),
174   fVariables(0x0),
175   fITSClusterMap(0),
176   fAvgTrials(1),
177   fNEventAll(0),
178   fNEventSel(0),
179   fNEventReject(0),
180   fh1Centrality(0x0),
181   fh1Xsec(0),
182   fh1Trials(0),
183   fh1PtHard(0),
184   fh1PtHardTrials(0),
185   fh1NTracksAll(0x0),
186   fh1NTracksReject(0x0),
187   fh1NTracksSel(0x0),
188   fPtAll(0),
189   fPtSel(0),
190   fPtPhi(0x0),
191   fPtEta(0x0),
192   fPtDCA2D(0x0),
193   fPtDCAZ(0x0),
194   fPtNClustersTPC(0x0),
195   fPtNClustersTPCIter1(0x0),
196   fPtNClustersTPCIter1Phi(0x0),
197   fPtNClustersTPCShared(0x0),
198   fPtNClustersTPCSharedFrac(0x0),
199   fPtNPointITS(0x0),
200   fPtChi2C(0x0),
201   fPtNSigmaToVertex(0x0),
202   fPtRelUncertainty1Pt(0x0),
203   fPtRelUncertainty1PtNClus(0x0),
204   fPtRelUncertainty1PtNClusIter1(0x0),
205   fPtRelUncertainty1PtNPointITS(0x0),
206   fPtRelUncertainty1PtITSClusterMap(0x0),
207   fPtRelUncertainty1PtChi2(0x0),
208   fPtRelUncertainty1PtChi2Iter1(0x0),
209   fPtRelUncertainty1PtPhi(0x0),
210   fPtUncertainty1Pt(0x0),
211   fPtChi2PerClusterTPC(0x0),
212   fPtChi2PerClusterTPCIter1(0x0),
213   fPtNCrossedRows(0x0),
214   fPtNCrossedRowsPhi(0x0),
215   fPtNCrossedRowsNClusFPhi(0x0),
216   fPtNCrRNCrRNClusF(0x0),
217   fPtNCrossedRowsFit(0x0),
218   fPtNCrossedRowsFitPhi(0x0),
219   fPtNCrossedRowsNClusFFitPhi(0x0),
220   fNCrossedRowsNCrossedRowsFit(0x0),
221   fNClustersNCrossedRows(0x0),
222   fNClustersNCrossedRowsFit(0x0),
223   fPtRelUncertainty1PtNCrossedRows(0x0),
224   fPtRelUncertainty1PtNCrossedRowsFit(0x0),
225   fPtChi2Gold(0x0),
226   fPtChi2GGC(0x0),
227   fPtChi2GoldPhi(0x0),
228   fPtChi2GGCPhi(0x0),
229   fChi2GoldChi2GGC(0x0),
230   fPtSigmaY2(0x0),
231   fPtSigmaZ2(0x0),
232   fPtSigmaSnp2(0x0),
233   fPtSigmaTgl2(0x0),
234   fPtSigma1Pt2(0x0),
235   fProfPtSigmaY2(0x0),
236   fProfPtSigmaZ2(0x0),
237   fProfPtSigmaSnp2(0x0),
238   fProfPtSigmaTgl2(0x0),
239   fProfPtSigma1Pt2(0x0),
240   fProfPtSigma1Pt(0x0),
241   fProfPtPtSigma1Pt(0x0),
242   fHistList(0)
243 {
244   //
245   // Constructor. Initialization of Inputs and Outputs
246   //
247   AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
248
249   SetNVariables(25);
250
251   fPtBinEdges[0][0] = 10.;
252   fPtBinEdges[0][1] = 1.;
253   fPtBinEdges[1][0] = 20.;
254   fPtBinEdges[1][1] = 2.;
255   fPtBinEdges[2][0] = 100.;
256   fPtBinEdges[2][1] = 5.;
257
258   // Input slot #0 works with a TChain ESD
259   DefineInput(0, TChain::Class());
260   // Output slot #1 write into a TList
261   DefineOutput(1, TList::Class());
262 }
263
264 //________________________________________________________________________
265 void AliPWG4HighPtTrackQA::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
266
267   if(region<3) {
268     fPtBinEdges[region][0] = ptmax;
269     fPtBinEdges[region][1] = ptBinWidth;
270   }
271   else {
272     AliError("Only 3 regions alowed. Use region 0/1/2\n");
273     return;
274   }
275
276 }
277
278 //________________________________________________________________________
279 void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
280   //Create output objects
281   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
282
283   Bool_t oldStatus = TH1::AddDirectoryStatus();
284   TH1::AddDirectory(kFALSE); 
285
286   OpenFile(1);
287   fHistList = new TList();
288   fHistList->SetOwner(kTRUE);
289   
290   Float_t fgkPtMin = 0.;
291   //  Float_t fgkPtMax = fPtMax;
292
293   //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
294   const Float_t ptmin1 =  fgkPtMin;
295   const Float_t ptmax1 =  fPtBinEdges[0][0];
296   const Float_t ptmin2 =  ptmax1 ;
297   const Float_t ptmax2 =  fPtBinEdges[1][0];
298   const Float_t ptmin3 =  ptmax2 ;
299   const Float_t ptmax3 =  fPtBinEdges[2][0];//fgkPtMax;
300   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
301   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
302   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
303   Int_t fgkNPtBins=nbin13;
304   //Create array with low edges of each bin
305   Double_t *binsPt=new Double_t[fgkNPtBins+1];
306   for(Int_t i=0; i<=fgkNPtBins; i++) {
307     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
308     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
309     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
310   }
311
312   Int_t fgkNPhiBins = 18*6;
313   Float_t kMinPhi   = 0.;
314   Float_t kMaxPhi   = 2.*TMath::Pi();
315   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
316   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
317
318   Int_t fgkNEtaBins=20;
319   Float_t fgkEtaMin = -1.;
320   Float_t fgkEtaMax = 1.;
321   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
322   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
323
324   Int_t fgkNNClustersTPCBins=80;
325   Float_t fgkNClustersTPCMin = 0.5;
326   Float_t fgkNClustersTPCMax = 160.5;
327   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
328   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
329
330   Int_t fgkNDCA2DBins=80;
331   Float_t fgkDCA2DMin = -0.2;
332   Float_t fgkDCA2DMax = 0.2;
333   if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==7) {
334     fgkDCA2DMin = -2.;
335     fgkDCA2DMax = 2.;
336   }
337   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
338   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
339
340   Int_t fgkNDCAZBins=80;
341   Float_t fgkDCAZMin = -2.;
342   Float_t fgkDCAZMax = 2.;
343   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
344     fgkDCAZMin = -5.;
345     fgkDCAZMax = 5.;
346   }
347   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
348   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
349
350   Int_t fgkNNPointITSBins=9;
351   Float_t fgkNPointITSMin = -0.5;
352   Float_t fgkNPointITSMax = 8.5;
353   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
354   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
355
356   Int_t fgkNITSClusterMapBins=65;
357   Float_t fgkITSClusterMapMin = -0.5;
358   Float_t fgkITSClusterMapMax = 64.5;
359   Double_t *binsITSClusterMap=new Double_t[fgkNITSClusterMapBins+1];
360   for(Int_t i=0; i<=fgkNITSClusterMapBins; i++) binsITSClusterMap[i]=(Double_t)fgkITSClusterMapMin + (fgkITSClusterMapMax-fgkITSClusterMapMin)/fgkNITSClusterMapBins*(Double_t)i ;
361
362
363   Int_t fgkNNSigmaToVertexBins=9;
364   Float_t fgkNSigmaToVertexMin = 0.;
365   Float_t fgkNSigmaToVertexMax = 9.;
366   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
367   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
368
369   Int_t fgkNChi2CBins=10;
370   //  Float_t fgkChi2CMin = 0.;
371   //  Float_t fgkChi2CMax = 100.; //10 sigma
372   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
373   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i] = (Double_t)i * (Double_t)i;
374
375   Float_t fgkRel1PtUncertaintyMin = 0.;
376   Float_t fgkRel1PtUncertaintyMax = 1.;
377   Float_t binEdgeRel1PtUncertainty1= 0.3;
378   Int_t fgkNRel1PtUncertaintyBins1 = 45;
379   Float_t binWidthRel1PtUncertainty1 = (binEdgeRel1PtUncertainty1-fgkRel1PtUncertaintyMin)/((Float_t)fgkNRel1PtUncertaintyBins1);
380   Int_t fgkNRel1PtUncertaintyBins2 = 35;
381   Float_t binWidthRel1PtUncertainty2 = (fgkRel1PtUncertaintyMax-binEdgeRel1PtUncertainty1)/((Float_t)fgkNRel1PtUncertaintyBins2);
382   Int_t fgkNRel1PtUncertaintyBins = fgkNRel1PtUncertaintyBins1 + fgkNRel1PtUncertaintyBins2;
383   
384   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
385   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) {
386     if(i<=fgkNRel1PtUncertaintyBins1)
387       binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (Double_t)binWidthRel1PtUncertainty1*(Double_t)i ;
388     if(i<=fgkNRel1PtUncertaintyBins && i>fgkNRel1PtUncertaintyBins1)
389       binsRel1PtUncertainty[i]=(Double_t)binEdgeRel1PtUncertainty1 + (Double_t)binWidthRel1PtUncertainty2*(Double_t)(i-fgkNRel1PtUncertaintyBins1);
390   }
391
392   Int_t fgkNUncertainty1PtBins = 30;
393   Float_t fgkUncertainty1PtMin = 0.;
394   Float_t fgkUncertainty1PtMax = 0.1;
395   if(fTrackType==1 || fTrackType==2 || fTrackType==4) 
396     fgkUncertainty1PtMax = 0.2; 
397   Double_t *binsUncertainty1Pt=new Double_t[fgkNUncertainty1PtBins+1];
398   for(Int_t i=0; i<=fgkNUncertainty1PtBins; i++) binsUncertainty1Pt[i]=(Double_t)fgkUncertainty1PtMin + (fgkUncertainty1PtMax-fgkUncertainty1PtMin)/fgkNUncertainty1PtBins*(Double_t)i ;
399
400   Float_t fgkChi2PerClusMin = 0.;
401   Float_t fgkChi2PerClusMax = 4.;
402   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
403   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
404   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
405
406   Int_t fgkNCrossedRowsNClusFBins  = 45;
407   Float_t fgkNCrossedRowsNClusFMin = 0.;
408   Float_t fgkNCrossedRowsNClusFMax = 1.5;
409   Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
410   for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
411
412   Float_t fgk1PtMin = 0.;
413   Float_t fgk1PtMax = 6.;
414   Float_t binEdge1Pt1 = 1.;
415   Float_t binWidth1Pt1 = 0.05;
416   Int_t fgkN1PtBins1 = (int)((binEdge1Pt1-fgk1PtMin)/binWidth1Pt1);
417   Float_t binWidth1Pt2 = 0.1;
418   Int_t fgkN1PtBins2 = (int)((fgk1PtMax-binEdge1Pt1)/binWidth1Pt2);
419   Int_t fgkN1PtBins = fgkN1PtBins1+fgkN1PtBins2;
420   Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
421
422   for(Int_t i=0; i<=fgkN1PtBins; i++) {
423     if(i<=fgkN1PtBins1) 
424       bins1Pt[i]=(Double_t)fgk1PtMin + (Double_t)(binEdge1Pt1-fgk1PtMin)/(Double_t)fgkN1PtBins1*(Double_t)i;
425     if(i<=fgkN1PtBins && i>fgkN1PtBins1)
426       bins1Pt[i]=(Double_t)binEdge1Pt1 + (Double_t)(fgk1PtMax-binEdge1Pt1)/(Double_t)fgkN1PtBins2*(Double_t)(i-fgkN1PtBins1);
427   }
428
429   Int_t fgkNSigmaY2Bins = 50;
430   Float_t fgkSigmaY2Min = 0.;
431   Float_t fgkSigmaY2Max = 1.;
432   if(fTrackType==1) fgkSigmaY2Max = 4.;
433   if(fTrackType==2 || fTrackType==4) fgkSigmaY2Max = 0.1;
434   Double_t *binsSigmaY2=new Double_t[fgkNSigmaY2Bins+1];
435   for(Int_t i=0; i<=fgkNSigmaY2Bins; i++) binsSigmaY2[i]=(Double_t)fgkSigmaY2Min + (fgkSigmaY2Max-fgkSigmaY2Min)/fgkNSigmaY2Bins*(Double_t)i ;
436
437   Int_t fgkNSigmaZ2Bins = 50;
438   Float_t fgkSigmaZ2Min = 0.;
439   Float_t fgkSigmaZ2Max = 0.4;
440   Double_t *binsSigmaZ2=new Double_t[fgkNSigmaZ2Bins+1];
441   for(Int_t i=0; i<=fgkNSigmaZ2Bins; i++) binsSigmaZ2[i]=(Double_t)fgkSigmaZ2Min + (fgkSigmaZ2Max-fgkSigmaZ2Min)/fgkNSigmaZ2Bins*(Double_t)i ;
442
443   Int_t fgkNSigmaSnp2Bins = 50;
444   Float_t fgkSigmaSnp2Min = 0.;
445   Float_t fgkSigmaSnp2Max = 0.05;
446   if(fTrackType==1) fgkSigmaSnp2Max = 0.2;
447   if(fTrackType==2 || fTrackType==4) fgkSigmaSnp2Max = 0.1;
448   Double_t *binsSigmaSnp2=new Double_t[fgkNSigmaSnp2Bins+1];
449   for(Int_t i=0; i<=fgkNSigmaSnp2Bins; i++) binsSigmaSnp2[i]=(Double_t)fgkSigmaSnp2Min + (fgkSigmaSnp2Max-fgkSigmaSnp2Min)/fgkNSigmaSnp2Bins*(Double_t)i ;
450
451   Int_t fgkNSigmaTgl2Bins = 50;
452   Float_t fgkSigmaTgl2Min = 0.;
453   Float_t fgkSigmaTgl2Max = 0.1;
454   if(fTrackType==1) fgkSigmaTgl2Max = 0.2;
455   if(fTrackType==2 || fTrackType==4) fgkSigmaTgl2Max = 0.1;
456   Double_t *binsSigmaTgl2=new Double_t[fgkNSigmaTgl2Bins+1];
457   for(Int_t i=0; i<=fgkNSigmaTgl2Bins; i++) binsSigmaTgl2[i]=(Double_t)fgkSigmaTgl2Min + (fgkSigmaTgl2Max-fgkSigmaTgl2Min)/fgkNSigmaTgl2Bins*(Double_t)i ;
458
459   Int_t fgkNSigma1Pt2Bins = 50;
460   Float_t fgkSigma1Pt2Min = 0.;
461   Float_t fgkSigma1Pt2Max = 1.;
462   Double_t *binsSigma1Pt2=new Double_t[fgkNSigma1Pt2Bins+1];
463   for(Int_t i=0; i<=fgkNSigma1Pt2Bins; i++) binsSigma1Pt2[i]=(Double_t)fgkSigma1Pt2Min + (fgkSigma1Pt2Max-fgkSigma1Pt2Min)/fgkNSigma1Pt2Bins*(Double_t)i ;
464
465
466   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
467   fHistList->Add(fNEventAll);
468   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
469   fHistList->Add(fNEventSel);
470   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
471   //Set labels
472   fNEventReject->Fill("noESD",0);
473   fNEventReject->Fill("Trigger",0);
474   fNEventReject->Fill("NTracks<2",0);
475   fNEventReject->Fill("noVTX",0);
476   fNEventReject->Fill("VtxStatus",0);
477   fNEventReject->Fill("NCont<2",0);
478   fNEventReject->Fill("ZVTX>10",0);
479   fNEventReject->Fill("cent",0);
480   fNEventReject->Fill("cent>90",0);
481   fHistList->Add(fNEventReject);
482
483   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
484   fHistList->Add(fh1Centrality);
485
486   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
487   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
488   fHistList->Add(fh1Xsec);
489
490   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
491   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
492   fHistList->Add(fh1Trials);
493
494   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
495   fHistList->Add(fh1PtHard);
496   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
497   fHistList->Add(fh1PtHardTrials);
498
499   fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
500   fHistList->Add(fh1NTracksAll);
501
502   fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
503   fh1NTracksReject->Fill("noESDtrack",0);
504   fh1NTracksReject->Fill("noTPCInner",0);
505   fh1NTracksReject->Fill("FillTPC",0);
506   fh1NTracksReject->Fill("noTPConly",0);
507   fh1NTracksReject->Fill("relate",0);
508   fh1NTracksReject->Fill("trackCuts",0);
509   fh1NTracksReject->Fill("laser",0);
510   fh1NTracksReject->Fill("chi2",0);
511   fHistList->Add(fh1NTracksReject);
512
513   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
514   fHistList->Add(fh1NTracksSel);
515
516   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
517   fHistList->Add(fPtAll);
518   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
519   fHistList->Add(fPtSel);
520
521   fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
522   fHistList->Add(fPtPhi);
523  
524   fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
525   fHistList->Add(fPtEta);
526  
527   fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
528   fHistList->Add(fPtDCA2D);
529  
530   fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
531   fHistList->Add(fPtDCAZ);
532  
533   fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
534   fHistList->Add(fPtNClustersTPC);
535
536   fPtNClustersTPCIter1 = new TH2F("fPtNClustersTPCIter1","fPtNClustersTPCIter1",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
537   fHistList->Add(fPtNClustersTPCIter1);
538
539   fPtNClustersTPCIter1Phi = new TH3F("fPtNClustersTPCIter1Phi","fPtNClustersTPCIter1Phi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
540   fHistList->Add(fPtNClustersTPCIter1Phi);
541
542   fPtNClustersTPCShared = new TH2F("fPtNClustersTPCShared","fPtNClustersTPCShared",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
543   fHistList->Add(fPtNClustersTPCShared);
544
545   fPtNClustersTPCSharedFrac = new TH2F("fPtNClustersTPCSharedFrac","fPtNClustersTPCSharedFrac",fgkNPtBins,binsPt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
546   fHistList->Add(fPtNClustersTPCSharedFrac);
547  
548   fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
549   fHistList->Add(fPtNPointITS);
550  
551   fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
552   fHistList->Add(fPtChi2C);
553  
554   fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
555   fHistList->Add(fPtNSigmaToVertex);
556
557   fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
558   fHistList->Add(fPtRelUncertainty1Pt);
559
560   fPtRelUncertainty1PtNClus = new TH3F("fPtRelUncertainty1PtNClus","fPtRelUncertainty1PtNClus",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
561   fHistList->Add(fPtRelUncertainty1PtNClus);
562
563   fPtRelUncertainty1PtNClusIter1 = new TH3F("fPtRelUncertainty1PtNClusIter1","fPtRelUncertainty1PtNClusIter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
564   fHistList->Add(fPtRelUncertainty1PtNClusIter1);
565
566   fPtRelUncertainty1PtNPointITS = new TH3F("fPtRelUncertainty1PtNPointITS","fPtRelUncertainty1PtNPointITS",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNPointITSBins,binsNPointITS);
567   fHistList->Add(fPtRelUncertainty1PtNPointITS);
568
569   fPtRelUncertainty1PtITSClusterMap = new TH3F("fPtRelUncertainty1PtITSClusterMap","fPtRelUncertainty1PtITSClusterMap",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNITSClusterMapBins,binsITSClusterMap);
570   fHistList->Add(fPtRelUncertainty1PtITSClusterMap);
571
572   fPtRelUncertainty1PtChi2 = new TH3F("fPtRelUncertainty1PtChi2","fPtRelUncertainty1PtChi2",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
573   fHistList->Add(fPtRelUncertainty1PtChi2);
574
575   fPtRelUncertainty1PtChi2Iter1 = new TH3F("fPtRelUncertainty1PtChi2Iter1","fPtRelUncertainty1PtChi2Iter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
576   fHistList->Add(fPtRelUncertainty1PtChi2Iter1);
577
578   fPtRelUncertainty1PtPhi = new TH3F("fPtRelUncertainty1PtPhi","fPtRelUncertainty1PtPhi",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNPhiBins,binsPhi);
579   fHistList->Add(fPtRelUncertainty1PtPhi);
580
581   fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
582   fHistList->Add(fPtUncertainty1Pt);
583  
584   fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
585   fHistList->Add(fPtChi2PerClusterTPC);
586  
587   fPtChi2PerClusterTPCIter1 = new TH2F("fPtChi2PerClusterTPCIter1","fPtChi2PerClusterTPCIter1",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
588   fHistList->Add(fPtChi2PerClusterTPCIter1);
589
590   fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
591   fHistList->Add(fPtNCrossedRows);
592
593   fPtNCrossedRowsPhi = new TH3F("fPtNCrossedRowsPhi","fPtNCrossedRowsPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
594   fHistList->Add(fPtNCrossedRowsPhi);
595
596   fPtNCrossedRowsNClusFPhi = new TH3F("fPtNCrossedRowsNClusFPhi","fPtNCrossedRowsNClusFPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
597   fHistList->Add(fPtNCrossedRowsNClusFPhi);
598  
599   fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
600   fHistList->Add(fPtNCrRNCrRNClusF);
601
602   fPtNCrossedRowsFit = new TH2F("fPtNCrossedRowsFit","fPtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
603   fHistList->Add(fPtNCrossedRowsFit);
604
605   fPtNCrossedRowsFitPhi = new TH3F("fPtNCrossedRowsFitPhi","fPtNCrossedRowsFitPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
606   fHistList->Add(fPtNCrossedRowsFitPhi);
607
608   fPtNCrossedRowsNClusFFitPhi = new TH3F("fPtNCrossedRowsNClusFFitPhi","fPtNCrossedRowsNClusFFitPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
609   fHistList->Add(fPtNCrossedRowsNClusFFitPhi);
610
611   fNCrossedRowsNCrossedRowsFit = new TH2F("fNCrossedRowsNCrossedRowsFit","fNCrossedRowsNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
612   fHistList->Add(fNCrossedRowsNCrossedRowsFit);
613
614   fNClustersNCrossedRows = new TH2F("fNClustersNCrossedRows","fNClustersNCrossedRows",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
615   fHistList->Add(fNClustersNCrossedRows);
616
617   fNClustersNCrossedRowsFit = new TH2F("fNClustersNCrossedRowsFit","fNClustersNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
618   fHistList->Add(fNClustersNCrossedRowsFit);
619
620   fPtRelUncertainty1PtNCrossedRows = new TH3F("fPtRelUncertainty1PtNCrossedRows","fPtRelUncertainty1PtNCrossedRows",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
621   fHistList->Add(fPtRelUncertainty1PtNCrossedRows);
622
623   fPtRelUncertainty1PtNCrossedRowsFit = new TH3F("fPtRelUncertainty1PtNCrossedRowsFit","fPtRelUncertainty1PtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
624   fHistList->Add(fPtRelUncertainty1PtNCrossedRowsFit);
625
626   fPtChi2Gold = new TH2F("fPtChi2Gold","fPtChi2Gold",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
627   fHistList->Add(fPtChi2Gold);
628
629   fPtChi2GGC = new TH2F("fPtChi2GGC","fPtChi2GGC",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
630   fHistList->Add(fPtChi2GGC);
631
632   fPtChi2GoldPhi = new TH3F("fPtChi2GoldPhi","fPtChi2GoldPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
633   fHistList->Add(fPtChi2GoldPhi);
634
635   fPtChi2GGCPhi = new TH3F("fPtChi2GGCPhi","fPtChi2GGCPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
636   fHistList->Add(fPtChi2GGCPhi);
637
638   fChi2GoldChi2GGC = new TH2F("fChi2GoldChi2GGC","fChi2GoldChi2GGC;#chi^{2}_{gold};#chi^{2}_{ggc}",fgkNChi2CBins,binsChi2C,fgkNChi2CBins,binsChi2C);
639   fHistList->Add(fChi2GoldChi2GGC);
640
641
642   fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
643   fHistList->Add(fPtSigmaY2);
644
645   fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
646   fHistList->Add(fPtSigmaZ2);
647
648   fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
649   fHistList->Add(fPtSigmaSnp2);
650
651   fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
652   fHistList->Add(fPtSigmaTgl2);
653
654   fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
655   fHistList->Add(fPtSigma1Pt2);
656
657   fProfPtSigmaY2 = new TProfile("fProfPtSigmaY2","fProfPtSigmaY2",fgkN1PtBins,bins1Pt);
658   fHistList->Add(fProfPtSigmaY2);
659
660   fProfPtSigmaZ2 = new TProfile("fProfPtSigmaZ2","fProfPtSigmaZ2",fgkN1PtBins,bins1Pt);
661   fHistList->Add(fProfPtSigmaZ2);
662
663   fProfPtSigmaSnp2 = new TProfile("fProfPtSigmaSnp2","fProfPtSigmaSnp2",fgkN1PtBins,bins1Pt);
664   fHistList->Add(fProfPtSigmaSnp2);
665
666   fProfPtSigmaTgl2 = new TProfile("fProfPtSigmaTgl2","fProfPtSigmaTgl2",fgkN1PtBins,bins1Pt);
667   fHistList->Add(fProfPtSigmaTgl2);
668
669   fProfPtSigma1Pt2 = new TProfile("fProfPtSigma1Pt2","fProfPtSigma1Pt2",fgkN1PtBins,bins1Pt);
670   fHistList->Add(fProfPtSigma1Pt2);
671
672   fProfPtSigma1Pt = new TProfile("fProfPtSigma1Pt","fProfPtSigma1Pt;p_{T};#sigma(1/p_{T})",fgkNPtBins,binsPt);
673   fHistList->Add(fProfPtSigma1Pt);
674
675   fProfPtPtSigma1Pt = new TProfile("fProfPtPtSigma1Pt","fProfPtPtSigma1Pt;p_{T};p_{T}#sigma(1/p_{T})",fgkNPtBins,binsPt);
676   fHistList->Add(fProfPtPtSigma1Pt);
677
678   TH1::AddDirectory(oldStatus); 
679
680   PostData(1, fHistList);
681
682   if(binsPhi)               delete [] binsPhi;
683   if(binsPt)                delete [] binsPt;
684   if(binsNClustersTPC)      delete [] binsNClustersTPC;
685   if(binsDCA2D)             delete [] binsDCA2D;
686   if(binsDCAZ)              delete [] binsDCAZ;
687   if(binsNPointITS)         delete [] binsNPointITS;
688   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
689   if(binsChi2C)             delete [] binsChi2C;
690   if(binsEta)               delete [] binsEta;
691   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
692   if(binsUncertainty1Pt)    delete [] binsUncertainty1Pt;
693   if(binsChi2PerClus)       delete [] binsChi2PerClus;
694   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
695   if(bins1Pt)               delete [] bins1Pt;
696   if(binsSigmaY2)           delete [] binsSigmaY2;
697   if(binsSigmaZ2)           delete [] binsSigmaZ2;
698   if(binsSigmaSnp2)         delete [] binsSigmaSnp2;
699   if(binsSigmaTgl2)         delete [] binsSigmaTgl2;
700   if(binsSigma1Pt2)         delete [] binsSigma1Pt2;
701 }
702
703 //________________________________________________________________________
704 Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
705   //
706   // Decide if event should be selected for analysis
707   //
708
709   // Checks following requirements:
710   // - fEvent available
711   // - trigger info from AliPhysicsSelection
712   // - MCevent available
713   // - number of reconstructed tracks > 1
714   // - primary vertex reconstructed
715   // - z-vertex < 10 cm
716   // - centrality in case of PbPb
717
718   Bool_t selectEvent = kTRUE;
719
720   //fEvent object available?
721   if (!fEvent) {
722     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
723     fNEventReject->Fill("noAliVEvent",1);
724     selectEvent = kFALSE;
725     return selectEvent;
726   }
727
728   //Check if number of reconstructed tracks is larger than 1
729   if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2)  {
730     fNEventReject->Fill("NTracks<2",1);
731     selectEvent = kFALSE;
732     return selectEvent;
733   }
734
735   //Check if vertex is reconstructed
736   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
737     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexTracks();
738
739     if (!fVtx || !fVtx->GetStatus())
740       fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
741
742     if(!fVtx) {
743       fNEventReject->Fill("noVTX",1);
744       selectEvent = kFALSE;
745       return selectEvent;
746     }
747     
748     if(!fVtx->GetStatus()) {
749       fNEventReject->Fill("VtxStatus",1);
750       selectEvent = kFALSE;
751       return selectEvent;
752     }
753     
754     // Need vertex cut
755     if(fVtx->GetNContributors()<2) {
756       fNEventReject->Fill("NCont<2",1);
757       selectEvent = kFALSE;
758       return selectEvent;
759     }
760     
761     //Check if z-vertex < 10 cm
762     double primVtx[3];
763     fVtx->GetXYZ(primVtx);
764     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
765       fNEventReject->Fill("ZVTX>10",1);
766       selectEvent = kFALSE;
767       return selectEvent;
768     }
769   }
770   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
771     const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
772     if(!vtx) {
773       fNEventReject->Fill("noVTX",1);
774       selectEvent = kFALSE;
775       return selectEvent;
776     }
777     
778     // Need vertex cut
779     if(vtx->GetNContributors()<2) {
780       fNEventReject->Fill("NCont<2",1);
781       selectEvent = kFALSE;
782       return selectEvent;
783     }
784     
785     //Check if z-vertex < 10 cm
786     double primVtx[3];
787     vtx->GetXYZ(primVtx);
788     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
789       fNEventReject->Fill("ZVTX>10",1);
790       selectEvent = kFALSE;
791       return selectEvent;
792     }
793
794   }
795
796   //Centrality selection should only be done in case of PbPb
797   if(IsPbPb()) {
798     Float_t cent = 0.;
799     if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
800       fNEventReject->Fill("cent",1);
801       selectEvent = kFALSE;
802       return selectEvent;
803     }
804     else {
805       if(fDataType==kESD) {
806         if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
807           cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
808         }
809       }
810       else if(fDataType==kAOD) {
811         if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
812           cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
813        }
814       if(cent>90.) {
815         fNEventReject->Fill("cent>90",1);
816         selectEvent = kFALSE;
817         return selectEvent;     
818       }
819       fh1Centrality->Fill(cent);
820     }
821   }
822  
823   return selectEvent;
824
825 }
826
827 //________________________________________________________________________
828 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
829   //
830   // Get centrality from ESD or AOD
831   //
832
833   if(fDataType==kESD)
834     return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
835   else if(fDataType==kAOD)
836     return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
837   else
838     return 5;
839 }
840
841 //________________________________________________________________________
842 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
843   //
844   // Get centrality from ESD
845   //
846
847   Float_t cent = -1;
848
849   if(esd){
850     if(esd->GetCentrality()){
851       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
852       if(fDebug>3) printf("centrality: %f\n",cent);
853     }
854   }
855
856   return GetCentralityClass(cent);
857
858 }
859
860 //________________________________________________________________________
861 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
862   //
863   // Get centrality from AOD
864   //
865
866   if(!aod) return 5;
867   Float_t cent = aod->GetHeader()->GetCentrality();
868   if(fDebug>3) printf("centrality: %f\n",cent);
869
870   return GetCentralityClass(cent);
871
872 }
873
874 //________________________________________________________________________
875 Int_t AliPWG4HighPtTrackQA::GetCentralityClass(Float_t cent) {
876   //
877   // Get centrality class
878   //
879
880   if(cent<0)  return 5; // OB - cent sometimes negative
881   if(cent>80) return 4;
882   if(cent>50) return 3;
883   if(cent>30) return 2;
884   if(cent>10) return 1;
885   return 0;
886
887 }
888
889 //________________________________________________________________________
890 void AliPWG4HighPtTrackQA::UserExec(Option_t *) {  
891   // Main loop
892   // Called for each event
893   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
894   
895   fEvent = InputEvent();
896   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
897
898   // All events without selection
899   fNEventAll->Fill(0.);
900
901   if(!SelectEvent()) {
902     // Post output data
903     PostData(1, fHistList);
904     return;
905   }
906
907
908   //Need to keep track of selected events
909   fNEventSel->Fill(0.);
910
911   fVariables = new TArrayF(fNVariables);
912   
913   if(fDataType==kESD) DoAnalysisESD();
914   if(fDataType==kAOD) DoAnalysisAOD();
915
916   //Delete old fVariables
917   if(fVariables) delete fVariables;
918
919   // Post output data
920   PostData(1, fHistList);
921
922 }
923
924 //________________________________________________________________________
925 void AliPWG4HighPtTrackQA::DoAnalysisESD() {
926   //
927   // Run analysis on ESD
928   //
929
930   if(!fESD) {
931     PostData(1, fHistList);
932     return;
933   }
934
935   // ---- Get MC Header information (for MC productions in pThard bins) ----
936   Double_t ptHard = 0.;
937   Double_t nTrials = 1; // trials for MC trigger weight for real data
938   
939   AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
940   if (eventHandlerMC) {
941     
942     if(eventHandlerMC->MCEvent()){
943       AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
944       if(pythiaGenHeader){
945         nTrials = pythiaGenHeader->Trials();
946         ptHard  = pythiaGenHeader->GetPtHard();
947         
948         fh1PtHard->Fill(ptHard);
949         fh1PtHardTrials->Fill(ptHard,nTrials);
950         
951         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
952       }
953     }
954   }
955
956   Int_t nTracks = fESD->GetNumberOfTracks();
957   AliDebug(2,Form("nTracks ESD%d", nTracks));
958
959   /*
960     Variables to be put in fVariables
961     0: pt
962     1: phi
963     2: eta
964     3: dca2D
965     4: dcaZ 
966     5: nClustersTPC
967     6: nPointITS   
968     7: chi2C       
969     8: nSigmaToVertex
970     9: trackLengthTPC
971     10: chi2PerClusterTPC
972     11: #crossed rows
973     12: (#crossed rows)/(#findable clusters)
974     13: SigmaY2
975     14: SigmaZ2
976     15: SigmaSnp2
977     16: SigmaTgl2
978     17: Sigma1Pt2
979     18: NClustersTPCIter1
980     19: Chi2TPCIter1
981     20: nClustersTPCShared
982     21: Golden Chi2 - global vs TPC constrained
983     22: Chi2 between global and global constrained
984   */
985
986   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
987     fh1NTracksAll->Fill(0.);
988
989     //Get track for analysis
990     AliESDtrack *track = 0x0;
991     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
992     if(!esdtrack) {
993       fh1NTracksReject->Fill("noESDtrack",1);
994       continue;
995     }
996     AliESDtrack *origtrack = new AliESDtrack(*esdtrack);
997     if(!origtrack)
998       continue;
999
1000     if(fTrackType==4) {
1001       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
1002         fh1NTracksReject->Fill("trackCuts",1);
1003         if(origtrack) delete origtrack;
1004         continue;
1005       }
1006     }
1007
1008     if(fTrackType==1)
1009       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1010     else if(fTrackType==2 || fTrackType==4) {
1011       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1012       if(!track) {
1013         fh1NTracksReject->Fill("noTPConly",1);
1014         if(origtrack) delete origtrack;
1015         continue;
1016       }
1017       AliExternalTrackParam exParam;
1018       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1019       if( !relate ) {
1020         fh1NTracksReject->Fill("relate",1);
1021         if(track) delete track;
1022         if(origtrack) delete origtrack;
1023         continue;
1024       }
1025       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1026     }
1027     else if(fTrackType==5 || fTrackType==6) {
1028       if(fTrackCuts->AcceptTrack(esdtrack)) {
1029         if(origtrack) delete origtrack;
1030         continue;
1031       }
1032       else {
1033         if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
1034
1035           if(fTrackType==5) {
1036             //use TPConly constrained track
1037             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1038             if(!track) {
1039               fh1NTracksReject->Fill("noTPConly",1);
1040               if(origtrack) delete origtrack;
1041               continue;
1042             }
1043             AliExternalTrackParam exParam;
1044             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1045             if( !relate ) {
1046               fh1NTracksReject->Fill("relate",1);
1047               if(track) delete track;
1048               if(origtrack) delete origtrack;
1049               continue;
1050             }
1051             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1052           }
1053           else if(fTrackType==6) {
1054             //use global constrained track
1055             track = new AliESDtrack(*esdtrack);
1056             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1057
1058           }
1059         }
1060       }
1061     }
1062     else if(fTrackType==7) {
1063       //use global constrained track
1064       track = new AliESDtrack(*esdtrack);
1065     }
1066     else
1067       track = esdtrack;
1068     
1069     if(!track) {
1070       if(origtrack) delete origtrack;
1071       continue;
1072     }
1073
1074     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
1075       //Cut on chi2 of constrained fit
1076       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
1077         fh1NTracksReject->Fill("chi2",1);
1078         if(track) delete track;
1079         if(origtrack) delete origtrack;
1080         continue;
1081       }
1082     }
1083
1084     fPtAll->Fill(track->Pt());
1085
1086     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
1087       fh1NTracksReject->Fill("trackCuts",1);
1088       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
1089         if(track) delete track;
1090       }
1091       if(origtrack) delete origtrack;
1092       continue;
1093     }
1094
1095     if(fTrackType==7) {
1096       if(fTrackCutsITSLoose ) {
1097         if(fTrackCutsITSLoose->AcceptTrack(track) ) {
1098           if(track) delete track;
1099           if(origtrack) delete origtrack;
1100           continue;
1101         }
1102       }
1103       
1104       if(esdtrack->GetConstrainedParam()) 
1105         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1106     }
1107
1108     if(!track) {
1109       if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1110         if(track) delete track;
1111       }
1112       if(origtrack) delete origtrack;
1113       continue;
1114     }
1115     
1116     fh1NTracksSel->Fill(0.);
1117
1118     fVariables->Reset(0.);
1119       
1120     fVariables->SetAt(track->Pt(),0); 
1121     fVariables->SetAt(track->Phi(),1); 
1122     fVariables->SetAt(track->Eta(),2); 
1123
1124     Float_t dca2D = 0.;
1125     Float_t dcaz  = 0.;
1126
1127     if(fTrackType==1 || fTrackType==2 || fTrackType==4) {  
1128       track->GetImpactParametersTPC(dca2D,dcaz); //TPConly
1129     }
1130     else
1131       track->GetImpactParameters(dca2D,dcaz);    //Global
1132
1133     fVariables->SetAt(dca2D,3);
1134     fVariables->SetAt(dcaz,4);
1135
1136     fVariables->SetAt((float)track->GetTPCNcls(),5);
1137
1138     Int_t nPointITS = 0;
1139     fITSClusterMap = track->GetITSClusterMap();
1140     UChar_t itsMap = track->GetITSClusterMap();
1141     for (Int_t i=0; i < 6; i++) {
1142       if (itsMap & (1 << i))
1143         nPointITS ++;
1144     }
1145     fVariables->SetAt((float)nPointITS,6);
1146     Float_t chi2C = (float)track->GetConstrainedChi2();
1147     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
1148       chi2C = (float)track->GetConstrainedChi2TPC();
1149     fVariables->SetAt(chi2C,7);
1150     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
1151   
1152     fVariables->SetAt(GetTrackLengthTPC(track),9);
1153   
1154     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
1155     
1156     //fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
1157     fVariables->SetAt(track->GetTPCCrossedRows(),11); //#crossed rows
1158
1159     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
1160     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
1161     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
1162     fVariables->SetAt(track->GetSigmaY2(),13);
1163     fVariables->SetAt(track->GetSigmaZ2(),14);
1164     fVariables->SetAt(track->GetSigmaSnp2(),15);
1165     fVariables->SetAt(track->GetSigmaTgl2(),16);
1166     fVariables->SetAt(track->GetSigma1Pt2(),17);
1167
1168     fVariables->SetAt(track->GetTPCNclsIter1(),18);
1169     fVariables->SetAt(track->GetTPCchi2Iter1(),19);
1170
1171     fVariables->SetAt(track->GetTPCnclsS(),20);
1172
1173     Float_t chi2Gold = origtrack->GetChi2TPCConstrainedVsGlobal(fVtx);//GetGoldenChi2(origtrack);
1174     Float_t chi2GGC  = GetGGCChi2(origtrack);
1175
1176     fVariables->SetAt(chi2Gold,21);
1177     fVariables->SetAt(chi2GGC,22);
1178
1179     fVariables->SetAt(GetTPCClusterInfoFitMap(track,2,1),23);
1180     Float_t crossedRowsTPCNClsFFit = 1.;
1181     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/track->GetTPCNclsF();
1182     fVariables->SetAt(crossedRowsTPCNClsFFit,24);
1183     
1184     FillHistograms();
1185   
1186     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1187
1188     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1189       if(track) delete track;
1190     }
1191     if(origtrack) delete origtrack;
1192     
1193   }//track loop
1194
1195 }
1196
1197 //________________________________________________________________________
1198 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1199
1200   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1201   if(!aod)return;
1202   AliExternalTrackParam *exParam = new  AliExternalTrackParam();
1203   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1204
1205     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1206     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1207
1208     fVariables->Reset(0.);
1209
1210     fVariables->SetAt(aodtrack->Pt(),0);
1211     fVariables->SetAt(aodtrack->Phi(),1);
1212     fVariables->SetAt(aodtrack->Eta(),2);
1213
1214     Double_t dca[2] = {1e6,1e6};
1215     Double_t covar[3] = {1e6,1e6,1e6};
1216     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1217       fVariables->SetAt(dca[0],3);
1218       fVariables->SetAt(dca[1],4);
1219     }
1220
1221     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1222     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1223     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1224     fVariables->SetAt(0.,8);
1225     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1226     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1227     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1),11);
1228     Float_t crossedRowsTPCNClsF = 0.;
1229     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1230     fVariables->SetAt(crossedRowsTPCNClsF,12);
1231
1232     //get covariance matrix
1233     Double_t cov[21] = {0,};
1234     aodtrack->GetCovMatrix(cov);
1235     Double_t pxpypz[3] = {0,};
1236     aodtrack->PxPyPz(pxpypz);
1237     Double_t xyz[3] = {0,};
1238     aodtrack->GetXYZ(xyz);
1239     Short_t sign = aodtrack->Charge();
1240     exParam->Set(xyz,pxpypz,cov,sign);
1241
1242     fVariables->SetAt(exParam->GetSigmaY2(),13);
1243     fVariables->SetAt(exParam->GetSigmaZ2(),14);
1244     fVariables->SetAt(exParam->GetSigmaSnp2(),15);
1245     fVariables->SetAt(exParam->GetSigmaTgl2(),16);
1246     fVariables->SetAt(exParam->GetSigma1Pt2(),17);
1247
1248     fVariables->SetAt(0.,18);
1249     fVariables->SetAt(0.,19);
1250
1251     TBits sharedClusterMap = aodtrack->GetTPCSharedMap();
1252     fVariables->SetAt(sharedClusterMap.CountBits(),20);
1253     
1254     fVariables->SetAt(0.,21); //not available in AOD
1255     fVariables->SetAt(0.,22); //not available in AOD
1256
1257     fVariables->SetAt(0.,23); //not available in AOD
1258     fVariables->SetAt(0.,24); //not available in AOD
1259
1260     fPtAll->Fill(fVariables->At(0));
1261
1262     FillHistograms();
1263
1264   }
1265
1266 }
1267
1268 //________________________________________________________________________
1269 void AliPWG4HighPtTrackQA::FillHistograms() {
1270
1271   fPtSel->Fill(fVariables->At(0));
1272   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1273   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1274   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1275   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1276   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1277   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1278
1279   
1280   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1281   fPtNClustersTPCIter1Phi->Fill(fVariables->At(0),fVariables->At(18),fVariables->At(1));
1282   fPtNClustersTPCShared->Fill(fVariables->At(0),fVariables->At(20));
1283   if(fVariables->At(5)>0.)
1284     fPtNClustersTPCSharedFrac->Fill(fVariables->At(0),fVariables->At(20)/fVariables->At(5));
1285
1286   if(fVariables->At(18)>0.)
1287     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1288   
1289   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1290   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1291   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1292   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1293   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1294   fPtRelUncertainty1PtNPointITS->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(6));
1295
1296   fPtRelUncertainty1PtITSClusterMap->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),(int)fITSClusterMap);
1297
1298   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1299   if(fVariables->At(18)>0.)
1300     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1301   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1302   
1303   fPtUncertainty1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1304   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1305   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1306   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1307   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1308   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1309
1310   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1311   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1312   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1313   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1314   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1315   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1316   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1317
1318   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1319   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1320   fPtNCrossedRowsPhi->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(1));
1321   fPtNCrossedRowsNClusFPhi->Fill(fVariables->At(0),fVariables->At(12),fVariables->At(1));
1322   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1323
1324   fPtChi2Gold->Fill(fVariables->At(0),fVariables->At(21));
1325   fPtChi2GGC->Fill(fVariables->At(0),fVariables->At(22));
1326
1327   fPtChi2GoldPhi->Fill(fVariables->At(0),fVariables->At(21),fVariables->At(1));
1328   fPtChi2GGCPhi->Fill(fVariables->At(0),fVariables->At(22),fVariables->At(1));
1329
1330   fChi2GoldChi2GGC->Fill(fVariables->At(21),fVariables->At(22));
1331
1332   fPtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(23));
1333   fPtNCrossedRowsFitPhi->Fill(fVariables->At(0),fVariables->At(23),fVariables->At(1));
1334   fPtNCrossedRowsNClusFFitPhi->Fill(fVariables->At(0),fVariables->At(24),fVariables->At(1));
1335   fNCrossedRowsNCrossedRowsFit->Fill(fVariables->At(11),fVariables->At(23));
1336
1337   fNClustersNCrossedRows->Fill(fVariables->At(5),fVariables->At(11));
1338   fNClustersNCrossedRowsFit->Fill(fVariables->At(5),fVariables->At(23));
1339
1340   fPtRelUncertainty1PtNCrossedRows->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(11));
1341   fPtRelUncertainty1PtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(23));
1342
1343 }
1344
1345 //________________________________________________________________________
1346 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1347   //
1348   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1349   // This is to called in Notify and should provide the path to the AOD/ESD file
1350   // Copied from AliAnalysisTaskJetSpectrum2
1351   //
1352
1353   TString file(currFile);  
1354   fXsec = 0;
1355   fTrials = 1;
1356
1357   if(file.Contains("root_archive.zip#")){
1358     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1359     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1360     file.Replace(pos+1,20,"");
1361   }
1362   else {
1363     // not an archive take the basename....
1364     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1365   }
1366   //  Printf("%s",file.Data());
1367    
1368
1369   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
1370   if(!fxsec){
1371     // next trial fetch the histgram file
1372     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1373     if(!fxsec){
1374         // not a severe condition but inciate that we have no information
1375       return kFALSE;
1376     }
1377     else{
1378       // find the tlist we want to be independtent of the name so use the Tkey
1379       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1380       if(!key){
1381         fxsec->Close();
1382         return kFALSE;
1383       }
1384       TList *list = dynamic_cast<TList*>(key->ReadObj());
1385       if(!list){
1386         fxsec->Close();
1387         return kFALSE;
1388       }
1389       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1390       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1391       fxsec->Close();
1392     }
1393   } // no tree pyxsec.root
1394   else {
1395     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1396     if(!xtree){
1397       fxsec->Close();
1398       return kFALSE;
1399     }
1400     UInt_t   ntrials  = 0;
1401     Double_t  xsection  = 0;
1402     xtree->SetBranchAddress("xsection",&xsection);
1403     xtree->SetBranchAddress("ntrials",&ntrials);
1404     xtree->GetEntry(0);
1405     fTrials = ntrials;
1406     fXsec = xsection;
1407     fxsec->Close();
1408   }
1409   return kTRUE;
1410 }
1411 //________________________________________________________________________
1412 Bool_t AliPWG4HighPtTrackQA::Notify()
1413 {
1414   //
1415   // Implemented Notify() to read the cross sections
1416   // and number of trials from pyxsec.root
1417   // Copied from AliAnalysisTaskJetSpectrum2
1418   // 
1419
1420   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1421   Float_t xsection = 0;
1422   Float_t ftrials  = 1;
1423
1424   fAvgTrials = 1;
1425   if(tree){
1426     TFile *curfile = tree->GetCurrentFile();
1427     if (!curfile) {
1428       Error("Notify","No current file");
1429       return kFALSE;
1430     }
1431     if(!fh1Xsec||!fh1Trials){
1432       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1433       return kFALSE;
1434     }
1435     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1436     fh1Xsec->Fill("<#sigma>",xsection);
1437     // construct a poor man average trials 
1438     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1439     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1440   }  
1441   return kTRUE;
1442 }
1443
1444 //________________________________________________________________________
1445 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1446   
1447   if(!mcEvent)return 0;
1448   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1449   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1450   if(!pythiaGenHeader){
1451     // cocktail ??
1452     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1453     
1454     if (!genCocktailHeader) {
1455       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1456       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1457       return 0;
1458     }
1459     TList* headerList = genCocktailHeader->GetHeaders();
1460     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1461       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1462       if (pythiaGenHeader)
1463         break;
1464     }
1465     if(!pythiaGenHeader){
1466       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1467       return 0;
1468     }
1469   }
1470   return pythiaGenHeader;
1471
1472 }
1473
1474 //_______________________________________________________________________
1475 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1476 {
1477   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1478
1479   //
1480   // TPC cluster information
1481   // type 0: get fraction of found/findable clusters with neighbourhood definition
1482   //      1: findable clusters with neighbourhood definition
1483   //      2: found clusters
1484   //
1485   // definition of findable clusters:
1486   //            a cluster is defined as findable if there is another cluster
1487   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1488   //           effects with a very simple algorithm.
1489   //
1490
1491   TBits fTPCClusterMap = tr->GetTPCClusterMap(); 
1492   if (type==2) return fTPCClusterMap.CountBits();
1493
1494   Int_t found=0;
1495   Int_t findable=0;
1496   Int_t last=-nNeighbours;
1497
1498   for (Int_t i=row0; i<row1; ++i){
1499     //look to current row
1500     if (fTPCClusterMap[i]) {
1501       last=i;
1502       ++found;
1503       ++findable;
1504       continue;
1505     }
1506     //look to nNeighbours before
1507     if ((i-last)<=nNeighbours) {
1508       ++findable;
1509       continue;
1510     }
1511     //look to nNeighbours after
1512     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1513       if (fTPCClusterMap[j]){
1514         ++findable;
1515         break;
1516       }
1517     }
1518   }
1519   if (type==1) return findable;
1520
1521   if (type==0){
1522     Float_t fraction=0;
1523     if (findable>0)
1524       fraction=(Float_t)found/(Float_t)findable;
1525     else
1526       fraction=0;
1527     return fraction;
1528   }
1529   return 0;  // undefined type - default value
1530 }
1531
1532 //_______________________________________________________________________
1533 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfoFitMap(AliESDtrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1534 {
1535   //
1536   // TPC cluster information from fit map
1537   // type 0: get fraction of found/findable clusters with neighbourhood definition
1538   //      1: findable clusters with neighbourhood definition
1539   //      2: found clusters
1540   //
1541   // definition of findable clusters:
1542   //            a cluster is defined as findable if there is another cluster
1543   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1544   //           effects with a very simple algorithm.
1545   //
1546
1547   TBits fTPCFitMap = tr->GetTPCFitMap();
1548   if (type==2) return fTPCFitMap.CountBits();
1549
1550   Int_t found=0;
1551   Int_t findable=0;
1552   Int_t last=-nNeighbours;
1553
1554   for (Int_t i=row0; i<row1; ++i){
1555     //look to current row
1556     if (fTPCFitMap[i]) {
1557       last=i;
1558       ++found;
1559       ++findable;
1560       continue;
1561     }
1562     //look to nNeighbours before
1563     if ((i-last)<=nNeighbours) {
1564       ++findable;
1565       continue;
1566     }
1567     //look to nNeighbours after
1568     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1569       if (fTPCFitMap[j]){
1570         ++findable;
1571         break;
1572       }
1573     }
1574   }
1575   if (type==1) return findable;
1576
1577   if (type==0){
1578     Float_t fraction=0;
1579     if (findable>0)
1580       fraction=(Float_t)found/(Float_t)findable;
1581     else
1582       fraction=0;
1583     return fraction;
1584   }
1585   return 0;  // undefined type - default value
1586 }
1587
1588 //_______________________________________________________________________
1589 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(AliESDtrack *track) {
1590   //
1591   // returns distance between 1st and last hit in TPC
1592   // distance given in number of padrows
1593   //
1594
1595   TBits fTPCClusterMap = track->GetTPCClusterMap(); 
1596   int firstHit = 0;
1597   int lastHit = 0;
1598
1599   for(int i=0; i<=159; i++) {
1600     if(fTPCClusterMap[i]>0) firstHit = i;
1601   }
1602   for(int i=159; i>=0; i--) {
1603     if(fTPCClusterMap[i]>0) lastHit = i;
1604   }
1605
1606   Int_t trackLength = lastHit - firstHit;
1607
1608   return trackLength;
1609 }
1610
1611 //_______________________________________________________________________
1612 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(AliAODTrack *track) {
1613   //
1614   // returns distance between 1st and last hit in TPC
1615   // distance given in number of padrows
1616   //
1617
1618   TBits fTPCClusterMap = track->GetTPCClusterMap(); 
1619   int firstHit = 0;
1620   int lastHit = 0;
1621
1622   for(int i=0; i<=159; i++) {
1623     if(fTPCClusterMap[i]>0) firstHit = i;
1624   }
1625   for(int i=159; i>=0; i--) {
1626     if(fTPCClusterMap[i]>0) lastHit = i;
1627   }
1628
1629   Int_t trackLength = lastHit - firstHit;
1630
1631   return trackLength;
1632 }
1633
1634 //_______________________________________________________________________
1635 Float_t AliPWG4HighPtTrackQA::GetGoldenChi2(AliESDtrack *origtrack) {
1636   //
1637   // Return chi2 between global and TPC constrained track
1638   // track should be the global unconstrained track
1639   //
1640
1641   Float_t chi2Gold = 0.;
1642
1643   AliESDtrack *tpcTrack = 0x0;
1644   tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,origtrack->GetID());
1645   if(tpcTrack) {
1646     AliExternalTrackParam exParam;
1647     Bool_t relate = tpcTrack->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1648     if( relate ) {
1649       tpcTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1650       //          Double_t pTPC[2],covTPC[3];     tpcTrack->PropagateToDCA(fVtx, fESD->GetMagneticField(), 10000,  pTPC, covTPC);
1651     }
1652   
1653     tpcTrack->Propagate(origtrack->GetAlpha(), origtrack->GetX(), fESD->GetMagneticField());
1654     chi2Gold = (Float_t)origtrack->GetPredictedChi2(tpcTrack);
1655   }
1656
1657   if(tpcTrack) delete tpcTrack;
1658
1659   return chi2Gold;
1660
1661 }
1662
1663 //_______________________________________________________________________
1664 Float_t AliPWG4HighPtTrackQA::GetGGCChi2(AliESDtrack *origtrack) {
1665   //
1666   // Return chi2 between global and global constrained track
1667   // track should be the global unconstrained track
1668   //
1669
1670   Float_t chi2GGC = 0.;
1671
1672   AliESDtrack *esdtrackC = new AliESDtrack(*origtrack);
1673   if(esdtrackC) {
1674     if(origtrack->GetConstrainedParam()) {
1675       esdtrackC->Set(origtrack->GetConstrainedParam()->GetX(),origtrack->GetConstrainedParam()->GetAlpha(),origtrack->GetConstrainedParam()->GetParameter(),origtrack->GetConstrainedParam()->GetCovariance());
1676       chi2GGC = (Float_t)origtrack->GetPredictedChi2(esdtrackC);
1677     }
1678     delete esdtrackC;
1679   }
1680   
1681   return chi2GGC;
1682
1683 }
1684
1685 //________________________________________________________________________
1686 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1687 {
1688   // The Terminate() function is the last function to be called during
1689   // a query. It always runs on the client, it can be used to present
1690   // the results graphically or save the results to file.
1691
1692 }
1693
1694 #endif