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