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