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