Add flag to select hybrids with ITS reffit (Leticia)
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtQAMC.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 compares the global reconstruction with the MC information
18 // Momentum resolution is stored as function of track cuts and pt.
19 // Output: Histograms for different set of cuts
20 //-----------------------------------------------------------------------
21 // Author : Marta Verweij - UU
22 //-----------------------------------------------------------------------
23
24 #ifndef ALIPWG4HighPtQAMC_CXX
25 #define ALIPWG4HighPtQAMC_CXX
26
27 #include "AliPWG4HighPtQAMC.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
40 #include "AliAnalysisTask.h"
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 "AliAnalysisHelperJetTasks.h"
53
54 using namespace std; //required for resolving the 'cout' symbol
55
56 ClassImp(AliPWG4HighPtQAMC)
57
58 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
59 : AliAnalysisTask("AliPWG4HighPtQAMC", ""), 
60   fESD(0), 
61   fMC(0),
62   fStack(0),
63   fVtx(0x0),
64   fTrackCuts(0), 
65   fTrackCutsReject(),
66   fTrackType(0),
67   fSigmaConstrainedMax(1e6),
68   fPtMax(100.),
69   fAvgTrials(1),
70   fNEventAll(0),
71   fNEventSel(0),
72   fNEventReject(0),
73   fh1Xsec(0),
74   fh1Trials(0),
75   fh1PtHard(0),
76   fh1PtHardTrials(0),
77   fPtAll(0),  
78   fPtSel(0),  
79   fPtSelFakes(0),
80   fNPointTPCFakes(0),
81   fPtSelLargeLabel(0),
82   fMultRec(0),
83   fNPointTPCMultRec(0),
84   fDeltaPtMultRec(0),
85   fPtAllvsPtMC(0),
86   fPtAllminPtMCvsPtMC(0),
87   fPtAllminPtMCvsPtAll(0),
88   fPtAllvsPtMCvsMult(0),
89   fPtAllminPtMCvsPtAllvsMult(0),
90   fPtAllminPtMCvsPtAllNPointTPC(0),
91   fPtAllminPtMCvsPtAllNPointTPCIter1(0),
92   fPtAllminPtMCvsPtAllChi2TPC(0),
93   fPtAllminPtMCvsPtAllChi2TPCIter1(0),
94   fPtAllminPtMCvsPtAllDCAR(0),
95   fPtAllminPtMCvsPtAllDCAZ(0),
96   fPtAllminPtMCvsPtAllPhi(0),
97   fPtAllminPtMCvsPtAllNPointITS(0),
98   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
99   fPtAllminPtMCvsPtAllChi2C(0),
100   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
101   fPtAllMC(0),
102   fPtSelMC(0),
103   fHistList(0)
104 {
105
106   fPtBinEdges[0][0] = 10.;
107   fPtBinEdges[0][1] = 1.;
108   fPtBinEdges[1][0] = 20.;
109   fPtBinEdges[1][1] = 2.;
110   fPtBinEdges[2][0] = 100.;
111   fPtBinEdges[2][1] = 10.;
112
113 }
114 //________________________________________________________________________
115 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name): 
116   AliAnalysisTask(name,""), 
117   fESD(0),
118   fMC(0),
119   fStack(0),
120   fVtx(0x0),
121   fTrackCuts(),
122   fTrackCutsReject(),
123   fTrackType(0),
124   fSigmaConstrainedMax(1e6),
125   fPtMax(100.),
126   fAvgTrials(1),
127   fNEventAll(0),
128   fNEventSel(0),
129   fNEventReject(0),
130   fh1Xsec(0),
131   fh1Trials(0),
132   fh1PtHard(0),
133   fh1PtHardTrials(0),
134   fPtAll(0),
135   fPtSel(0),
136   fPtSelFakes(0),
137   fNPointTPCFakes(0),
138   fPtSelLargeLabel(0),
139   fMultRec(0),
140   fNPointTPCMultRec(0),
141   fDeltaPtMultRec(0),
142   fPtAllvsPtMC(0),
143   fPtAllminPtMCvsPtMC(0),
144   fPtAllminPtMCvsPtAll(0),
145   fPtAllvsPtMCvsMult(0),
146   fPtAllminPtMCvsPtAllvsMult(0),
147   fPtAllminPtMCvsPtAllNPointTPC(0),
148   fPtAllminPtMCvsPtAllNPointTPCIter1(0),
149   fPtAllminPtMCvsPtAllChi2TPC(0),
150   fPtAllminPtMCvsPtAllChi2TPCIter1(0),
151   fPtAllminPtMCvsPtAllDCAR(0),
152   fPtAllminPtMCvsPtAllDCAZ(0),
153   fPtAllminPtMCvsPtAllPhi(0),
154   fPtAllminPtMCvsPtAllNPointITS(0),
155   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
156   fPtAllminPtMCvsPtAllChi2C(0),
157   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
158   fPtAllMC(0),
159   fPtSelMC(0),
160   fHistList(0)
161 {
162   //
163   // Constructor. Initialization of Inputs and Outputs
164   //
165   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
166
167   fPtBinEdges[0][0] = 10.;
168   fPtBinEdges[0][1] = 1.;
169   fPtBinEdges[1][0] = 20.;
170   fPtBinEdges[1][1] = 2.;
171   fPtBinEdges[2][0] = 100.;
172   fPtBinEdges[2][1] = 10.;
173
174   // Input slot #0 works with a TChain ESD
175   DefineInput(0, TChain::Class());
176   // Output slot #0, #1 write into a TList
177   DefineOutput(0, TList::Class());
178 }
179
180 //________________________________________________________________________
181 void AliPWG4HighPtQAMC::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
182
183   if(region<3) {
184     fPtBinEdges[region][0] = ptmax;
185     fPtBinEdges[region][1] = ptBinWidth;
186   }
187   else {
188     AliError("Only 3 regions alowed. Use region 0/1/2\n");
189     return;
190   }
191
192 }
193
194 //________________________________________________________________________
195 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *) 
196 {
197   // Connect ESD and MC event handler here
198   // Called once
199   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
200   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
201   if (!tree) {
202     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
203     return;
204   }
205
206   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
207
208   if (!esdH) {
209     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
210     return;
211   } else
212     fESD = esdH->GetEvent();
213   
214   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
215   if (!eventHandler) {
216     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
217   }
218   else
219     fMC = eventHandler->MCEvent();
220
221 }
222
223
224 //________________________________________________________________________
225 void AliPWG4HighPtQAMC::CreateOutputObjects() {
226   //Create output objects
227   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
228
229   Bool_t oldStatus = TH1::AddDirectoryStatus();
230   TH1::AddDirectory(kFALSE); 
231
232   OpenFile(0);
233   fHistList = new TList();
234   fHistList->SetOwner(kTRUE);
235
236   Float_t fgkPtMin = 0.;
237
238   //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
239   const Float_t ptmin1 =  fgkPtMin;
240   const Float_t ptmax1 =  fPtBinEdges[0][0];
241   const Float_t ptmin2 =  ptmax1 ;
242   const Float_t ptmax2 =  fPtBinEdges[1][0];
243   const Float_t ptmin3 =  ptmax2 ;
244   const Float_t ptmax3 =  fPtBinEdges[2][0];//fgkPtMax;
245   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
246   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
247   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
248   Int_t fgkNPtBins=nbin13;
249   //Create array with low edges of each bin
250   Double_t *binsPt=new Double_t[fgkNPtBins+1];
251   for(Int_t i=0; i<=fgkNPtBins; i++) {
252     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
253     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
254     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
255   }
256
257   Int_t fgkNPhiBins = 18*6;
258   Float_t kMinPhi   = 0.;
259   Float_t kMaxPhi   = 2.*TMath::Pi();
260   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
261   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
262
263   Int_t fgkResPtBins=80;
264   Float_t kMinResPt = -1.;
265   Float_t kMaxResPt = 1.;
266   Double_t *binsResPt = new Double_t[fgkResPtBins+1];
267   for(Int_t i=0; i<=fgkResPtBins; i++) binsResPt[i]=(Double_t)kMinResPt + (kMaxResPt-kMinResPt)/fgkResPtBins*(Double_t)i ;
268
269   Int_t fgkMultBins=50;
270   Float_t kMinMult = 0.;
271   Float_t kMaxMult = 3000.;
272   Double_t *binsMult = new Double_t[fgkMultBins+1];
273   for(Int_t i=0; i<=fgkMultBins; i++) binsMult[i]=(Double_t)kMinMult + (kMaxMult-kMinMult)/fgkMultBins*(Double_t)i ;
274
275   Int_t fgkNEtaBins=20;
276   Float_t fgkEtaMin = -1.;
277   Float_t fgkEtaMax = 1.;
278   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
279   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
280
281   Int_t fgkNNClustersTPCBins=80;
282   Float_t fgkNClustersTPCMin = 0.5;
283   Float_t fgkNClustersTPCMax = 160.5;
284   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
285   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
286
287   Int_t fgkNDCA2DBins=80;
288   Float_t fgkDCA2DMin = -0.2;
289   Float_t fgkDCA2DMax = 0.2;
290   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
291     fgkDCA2DMin = -2.;
292     fgkDCA2DMax = 2.;
293   }
294   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
295   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
296
297   Int_t fgkNDCAZBins=80;
298   Float_t fgkDCAZMin = -2.;
299   Float_t fgkDCAZMax = 2.;
300   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
301     fgkDCAZMin = -5.;
302     fgkDCAZMax = 5.;
303   }
304   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
305   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
306
307  Int_t fgkNNPointITSBins=9;
308   Float_t fgkNPointITSMin = -0.5;
309   Float_t fgkNPointITSMax = 8.5;
310   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
311   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
312
313   Int_t fgkNNSigmaToVertexBins=20;
314   Float_t fgkNSigmaToVertexMin = 0.;
315   Float_t fgkNSigmaToVertexMax = 8.;
316   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
317   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
318
319   Int_t fgkNChi2CBins=20;
320   Float_t fgkChi2CMin = 0.;
321   Float_t fgkChi2CMax = 100.;
322   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
323   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
324
325   Int_t fgkNRel1PtUncertaintyBins=50;
326   Float_t fgkRel1PtUncertaintyMin = 0.;
327   Float_t fgkRel1PtUncertaintyMax = 1.;
328   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
329   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
330
331   Float_t fgkChi2PerClusMin = 0.;
332   Float_t fgkChi2PerClusMax = 4.;
333   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
334   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
335   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
336
337
338   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
339   fHistList->Add(fNEventAll);
340   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
341   fHistList->Add(fNEventSel);
342   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
343   //Set labels
344   fNEventReject->Fill("noESD",0);
345   fNEventReject->Fill("Trigger",0);
346   fNEventReject->Fill("noMCEvent",0);
347   fNEventReject->Fill("noStack",0);
348   fNEventReject->Fill("NTracks<2",0);
349   fNEventReject->Fill("noVTX",0);
350   fNEventReject->Fill("VtxStatus",0);
351   fNEventReject->Fill("NCont<2",0);
352   fNEventReject->Fill("ZVTX>10",0);
353   fHistList->Add(fNEventReject);
354
355   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
356   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
357   fHistList->Add(fh1Xsec);
358
359   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
360   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
361   fHistList->Add(fh1Trials);
362
363   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
364   fHistList->Add(fh1PtHard);
365   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
366   fHistList->Add(fh1PtHardTrials);
367
368   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
369   fHistList->Add(fPtAll);
370   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
371   fHistList->Add(fPtSel);
372
373   fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, binsPt);
374   fHistList->Add(fPtSelFakes);
375
376   fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",fgkNNClustersTPCBins, binsNClustersTPC);
377   fHistList->Add(fNPointTPCFakes);
378
379   fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, binsPt);
380   fHistList->Add(fPtSelLargeLabel);
381
382   fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",fgkMultBins, binsMult);
383   fHistList->Add(fMultRec);
384
385   fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",fgkNNClustersTPCBins, binsNClustersTPC);
386   fHistList->Add(fNPointTPCMultRec);
387
388   fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",50,0.,50.,40,-20.,20.);
389   fHistList->Add(fDeltaPtMultRec);
390
391   fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, binsPt,fgkNPtBins, binsPt);
392   fHistList->Add(fPtAllvsPtMC);
393
394   fPtAllminPtMCvsPtMC = new TH2F("fPtAllminPtMCvsPtMC","PtAllminPtMCvsPtMC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
395   fPtAllminPtMCvsPtMC->SetXTitle("p_{t}^{MC}");
396   fPtAllminPtMCvsPtMC->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
397   fHistList->Add(fPtAllminPtMCvsPtMC);
398
399   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
400   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{rec}");
401   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
402   fHistList->Add(fPtAllminPtMCvsPtAll);
403
404   fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, binsPt,fgkNPtBins, binsPt,fgkMultBins,binsMult);
405   fHistList->Add(fPtAllvsPtMCvsMult);
406
407   fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkMultBins,binsMult);
408   fPtAllminPtMCvsPtAllvsMult->SetXTitle("p_{t}^{MC}");
409   fPtAllminPtMCvsPtAllvsMult->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
410   fPtAllminPtMCvsPtAllvsMult->SetZTitle("N_{tracks}");
411   fHistList->Add(fPtAllminPtMCvsPtAllvsMult);
412   
413   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
414   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
415   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
416   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{cls,TPC}");
417   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
418
419   fPtAllminPtMCvsPtAllNPointTPCIter1 = new TH3F("fPtAllminPtMCvsPtAllNPointTPCIter1","PtAllminPtMCvsPtAllNPointTPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
420   fPtAllminPtMCvsPtAllNPointTPCIter1->SetXTitle("p_{t}^{MC}");
421   fPtAllminPtMCvsPtAllNPointTPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
422   fPtAllminPtMCvsPtAllNPointTPCIter1->SetZTitle("N_{clsIter1,TPC}");
423   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPCIter1);
424
425   fPtAllminPtMCvsPtAllChi2TPC = new TH3F("fPtAllminPtMCvsPtAllChi2TPC","PtAllminPtMCvsPtAllChi2TPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
426   fPtAllminPtMCvsPtAllChi2TPC->SetXTitle("p_{t}^{MC}");
427   fPtAllminPtMCvsPtAllChi2TPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
428   fPtAllminPtMCvsPtAllChi2TPC->SetZTitle("#chi^{2} TPC");
429   fHistList->Add(fPtAllminPtMCvsPtAllChi2TPC);
430
431   fPtAllminPtMCvsPtAllChi2TPCIter1 = new TH3F("fPtAllminPtMCvsPtAllChi2TPCIter1","PtAllminPtMCvsPtAllChi2TPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
432   fPtAllminPtMCvsPtAllChi2TPCIter1->SetXTitle("p_{t}^{MC}");
433   fPtAllminPtMCvsPtAllChi2TPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
434   fPtAllminPtMCvsPtAllChi2TPCIter1->SetZTitle("#chi^{2} TPC Iter1");
435   fHistList->Add(fPtAllminPtMCvsPtAllChi2TPCIter1);
436
437   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCA2DBins,binsDCA2D);
438   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
439   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
440   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
441   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
442
443   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCAZBins,binsDCAZ);
444   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
445   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
446   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
447   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
448
449   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNPhiBins,binsPhi);
450   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
451   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
452   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
453   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
454
455   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNPointITSBins,binsNPointITS);
456   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
457   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
458   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
459   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
460   
461   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
462   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
463   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
464   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
465   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
466
467   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2CBins,binsChi2C);
468   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
469   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
470   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
471   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
472
473   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
474   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
475   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
476   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
477   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
478
479
480   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, binsPt);
481   fHistList->Add(fPtAllMC);
482   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, binsPt);
483   fHistList->Add(fPtSelMC);
484
485   TH1::AddDirectory(oldStatus); 
486
487   PostData(0, fHistList);
488
489   if(binsPt)                delete [] binsPt;
490   if(binsResPt)             delete [] binsResPt;
491   if(binsPhi)               delete [] binsPhi;
492   if(binsNClustersTPC)      delete [] binsNClustersTPC;
493   if(binsDCA2D)             delete [] binsDCA2D;
494   if(binsDCAZ)              delete [] binsDCAZ;
495   if(binsNPointITS)         delete [] binsNPointITS;
496   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
497   if(binsChi2C)             delete [] binsChi2C;
498   if(binsEta)               delete [] binsEta;
499   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
500   if(binsChi2PerClus)       delete [] binsChi2PerClus;
501   if(binsMult)              delete [] binsMult;
502   
503 }
504
505 //________________________________________________________________________
506 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
507   //
508   // Decide if event should be selected for analysis
509   //
510
511   // Checks following requirements:
512   // - fESD available
513   // - trigger info from AliPhysicsSelection
514   // - MCevent available
515   // - number of reconstructed tracks > 1
516   // - primary vertex reconstructed
517   // - z-vertex < 10 cm
518
519   Bool_t selectEvent = kTRUE;
520
521   //fESD object available?
522   if (!fESD) {
523     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
524     fNEventReject->Fill("noESD",1);
525     selectEvent = kFALSE;
526     return selectEvent;
527   }
528
529   //Trigger
530   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
531   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
532     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
533     fNEventReject->Fill("Trigger",1);
534     selectEvent = kFALSE;
535     return selectEvent;
536   }
537
538   //MCEvent available?
539   //if yes: get stack
540   if(fMC) {
541     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
542     fStack = fMC->Stack();                //Particles Stack
543     if(fStack) { 
544       AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
545     }
546     else {
547       AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
548       fNEventReject->Fill("noStack",1);
549       selectEvent = kFALSE;
550       return selectEvent;
551     }
552   } else {
553     AliDebug(2,Form("ERROR: Could not retrieve stack"));
554     fNEventReject->Fill("noMCEvent",1);
555     selectEvent = kFALSE;
556     return selectEvent;
557   }
558
559   //Check if number of reconstructed tracks is larger than 1
560   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
561     fNEventReject->Fill("NTracks<2",1);
562     selectEvent = kFALSE;
563     return selectEvent;
564   }
565
566   //Check if vertex is reconstructed
567   if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
568   else              fVtx = fESD->GetPrimaryVertexSPD();
569
570   if(!fVtx) {
571     fNEventReject->Fill("noVTX",1);
572     selectEvent = kFALSE;
573     return selectEvent;
574   }
575    
576   if(!fVtx->GetStatus()) {
577     fNEventReject->Fill("VtxStatus",1);
578     selectEvent = kFALSE;
579     return selectEvent;
580   }
581   
582   // Need vertex cut
583   //  TString vtxName(fVtx->GetName());
584   if(fVtx->GetNContributors()<2) {
585     fNEventReject->Fill("NCont<2",1);
586     selectEvent = kFALSE;
587     return selectEvent;
588   }
589   
590   //Check if z-vertex < 10 cm
591   double primVtx[3];
592   fVtx->GetXYZ(primVtx);
593   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
594     fNEventReject->Fill("ZVTX>10",1);
595     selectEvent = kFALSE;
596     return selectEvent;
597   }
598
599   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
600
601   return selectEvent;
602
603 }
604
605 //________________________________________________________________________
606 void AliPWG4HighPtQAMC::Exec(Option_t *) {  
607   // Main loop
608   // Called for each event
609   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));  
610   // All events without selection
611   fNEventAll->Fill(0.);
612
613   if(!SelectEvent()) {
614     // Post output data
615     PostData(0, fHistList);
616     return;
617   }
618
619   // ---- Get MC Header information (for MC productions in pThard bins) ----
620   Double_t ptHard = 0.;
621   Double_t nTrials = 1; // trials for MC trigger weight for real data
622   
623   if(fMC){
624     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
625      if(pythiaGenHeader){
626        nTrials = pythiaGenHeader->Trials();
627        ptHard  = pythiaGenHeader->GetPtHard();
628        
629        fh1PtHard->Fill(ptHard);
630        fh1PtHardTrials->Fill(ptHard,nTrials);
631        
632        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
633      }
634   }
635
636   //Need to keep track of selected events
637   fNEventSel->Fill(0.);
638
639   Int_t nTracks = fESD->GetNumberOfTracks();
640   AliDebug(2,Form("nTracks ESD%d", nTracks));
641
642   int nMCtracks = fStack->GetNtrack();
643
644   Float_t pt      = 0.;
645   Float_t ptMC    = 0.;
646   Float_t phi     = 0.;
647   Float_t dca2D   = 0.;
648   Float_t dcaZ    = 0.;
649   Int_t nPointITS = 0;
650   Float_t chi2C   = 0.;
651   Float_t nSigmaToVertex    = 0.;
652   Float_t relUncertainty1Pt = 0.;
653
654   int mult = fTrackCuts->CountAcceptedTracks(fESD);
655
656   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
657     
658     AliESDtrack *track = 0;
659     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
660     if(!esdtrack) continue;
661
662     if(fTrackType==1) {
663       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
664       if(!track) continue;
665     }
666     else if(fTrackType==2) {
667       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
668       if(!track) continue;
669       
670       AliExternalTrackParam exParam;
671       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
672       if( !relate ) {
673         delete track;
674         continue;
675       }
676       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
677     }
678     else if(fTrackType==7) {
679       //use global constrained track
680       track = new AliESDtrack(*esdtrack);
681       //      track = esdtrack;
682       //      track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
683     }
684     else
685       track = esdtrack;
686
687     if(!track) continue;
688     
689     if(fTrackType==2) {
690       //Cut on chi2 of constrained fit
691       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
692         delete track;
693         continue;
694       }
695     }
696
697     Int_t label = TMath::Abs(track->GetLabel());
698     if(label>=nMCtracks) {
699       if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
700       if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
701       continue;
702     }
703
704     TParticle *particle = fStack->Particle(label) ;
705     if(!particle) {
706       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
707         if(track) delete track;
708       }
709       continue;
710     }
711
712     ptMC = particle->Pt();
713
714     if(fTrackType==1 || fTrackType==2) 
715       track->GetImpactParametersTPC(dca2D,dcaZ);  //TPConly
716     else
717       track->GetImpactParameters(dca2D,dcaZ);     //Global
718     
719     UChar_t itsMap = track->GetITSClusterMap();
720     for (Int_t i=0; i < 6; i++) {
721       if (itsMap & (1 << i))
722         nPointITS ++;
723     }
724     
725     //    fPtAll->Fill(pt);
726     fPtAllMC->Fill(ptMC);
727
728     if (fTrackCuts->AcceptTrack(track)) {
729
730       if(fTrackType==7) {
731         if(fTrackCutsReject ) {
732           if(fTrackCutsReject->AcceptTrack(track) ) {
733             if(track) delete track;
734             continue;
735           }
736         }
737         
738         if(esdtrack->GetConstrainedParam()) 
739           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
740       }
741
742       pt  = track->Pt();
743       phi = track->Phi();
744
745       nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
746       chi2C = track->GetConstrainedChi2();
747       relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
748
749       fPtSel->Fill(pt);
750       if(track->GetLabel()<0) {
751         fPtSelFakes->Fill(pt);
752         fNPointTPCFakes->Fill(track->GetTPCNcls());
753       }
754       fPtSelMC->Fill(ptMC);
755       fPtAllvsPtMC->Fill(ptMC,pt);
756       fPtAllminPtMCvsPtMC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
757       fPtAllminPtMCvsPtAll->Fill(pt,(1./pt-1./ptMC)/(1./ptMC) );
758       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
759       fPtAllminPtMCvsPtAllNPointTPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNclsIter1());
760       if(track->GetTPCNcls()>0.) 
761         fPtAllminPtMCvsPtAllChi2TPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2()/track->GetTPCNcls());
762       if(track->GetTPCNclsIter1()>0.) 
763         fPtAllminPtMCvsPtAllChi2TPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2Iter1()/track->GetTPCNclsIter1());
764       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
765       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
766       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
767       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
768       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
769       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
770       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
771
772       fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
773       fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
774
775       //Check if track is reconstructed multiple times
776       /*
777       int multCounter = 1;
778       for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
779         //   AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
780         AliESDtrack *track2;
781         AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
782         if(!esdtrack2) continue;
783       
784         if(fTrackType==1)
785           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
786         else if(fTrackType==2) {
787           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
788           if(!track2) { 
789             continue;
790           }
791           AliExternalTrackParam exParam2;
792           Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
793           if( !relate ) {
794             delete track2;
795             continue;
796           }
797           track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
798         }
799         else
800           track2 = esdtrack2;
801         if(!track2) {
802           continue;
803         }
804       
805         if (fTrackCuts->AcceptTrack(track2)) {
806           Int_t label2 = TMath::Abs(track2->GetLabel());
807           if(label==label2) {
808             fNPointTPCMultRec->Fill(track->GetTPCNcls());
809             fNPointTPCMultRec->Fill(track2->GetTPCNcls());
810             fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
811             multCounter++;
812           }
813         }
814         if(fTrackType==1 || fTrackType==2) delete track2;
815       }//track2 loop
816       if(multCounter>1) fMultRec->Fill(multCounter);
817       */
818
819     }//fTrackCuts selection
820
821
822     if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
823       if(track) delete track;    
824     }
825
826   }//ESD track loop
827    
828   // Post output data
829   PostData(0, fHistList);
830   
831 }
832 //________________________________________________________________________
833 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
834   //
835   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
836   // This is to called in Notify and should provide the path to the AOD/ESD file
837   // Copied from AliAnalysisTaskJetSpectrum2
838   //
839
840   TString file(currFile);  
841   fXsec = 0;
842   fTrials = 1;
843
844   if(file.Contains("root_archive.zip#")){
845     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
846     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
847     file.Replace(pos+1,20,"");
848   }
849   else {
850     // not an archive take the basename....
851     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
852   }
853   //  Printf("%s",file.Data());
854    
855
856   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
857   if(!fxsec){
858     // next trial fetch the histgram file
859     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
860     if(!fxsec){
861         // not a severe condition but inciate that we have no information
862       return kFALSE;
863     }
864     else{
865       // find the tlist we want to be independtent of the name so use the Tkey
866       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
867       if(!key){
868         fxsec->Close();
869         return kFALSE;
870       }
871       TList *list = dynamic_cast<TList*>(key->ReadObj());
872       if(!list){
873         fxsec->Close();
874         return kFALSE;
875       }
876       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
877       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
878       fxsec->Close();
879     }
880   } // no tree pyxsec.root
881   else {
882     TTree *xtree = (TTree*)fxsec->Get("Xsection");
883     if(!xtree){
884       fxsec->Close();
885       return kFALSE;
886     }
887     UInt_t   ntrials  = 0;
888     Double_t  xsection  = 0;
889     xtree->SetBranchAddress("xsection",&xsection);
890     xtree->SetBranchAddress("ntrials",&ntrials);
891     xtree->GetEntry(0);
892     fTrials = ntrials;
893     fXsec = xsection;
894     fxsec->Close();
895   }
896   return kTRUE;
897 }
898 //________________________________________________________________________
899 Bool_t AliPWG4HighPtQAMC::Notify()
900 {
901   //
902   // Implemented Notify() to read the cross sections
903   // and number of trials from pyxsec.root
904   // Copied from AliAnalysisTaskJetSpectrum2
905   // 
906
907   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
908   Float_t xsection = 0;
909   Float_t ftrials  = 1;
910
911   fAvgTrials = 1;
912   if(tree){
913     TFile *curfile = tree->GetCurrentFile();
914     if (!curfile) {
915       Error("Notify","No current file");
916       return kFALSE;
917     }
918     if(!fh1Xsec||!fh1Trials){
919       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
920       return kFALSE;
921     }
922     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
923     fh1Xsec->Fill("<#sigma>",xsection);
924     // construct a poor man average trials 
925     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
926     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
927   }  
928   return kTRUE;
929 }
930
931 //________________________________________________________________________
932 AliGenPythiaEventHeader*  AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
933   
934   if(!mcEvent)return 0;
935   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
936   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
937   if(!pythiaGenHeader){
938     // cocktail ??
939     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
940     
941     if (!genCocktailHeader) {
942       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
943       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
944       return 0;
945     }
946     TList* headerList = genCocktailHeader->GetHeaders();
947     for (Int_t i=0; i<headerList->GetEntries(); i++) {
948       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
949       if (pythiaGenHeader)
950         break;
951     }
952     if(!pythiaGenHeader){
953       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
954       return 0;
955     }
956   }
957   return pythiaGenHeader;
958
959 }
960
961 //________________________________________________________________________
962 void AliPWG4HighPtQAMC::Terminate(Option_t *)
963 {
964   // The Terminate() function is the last function to be called during
965   // a query. It always runs on the client, it can be used to present
966   // the results graphically or save the results to file.
967
968 }
969
970 #endif