]> 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
Coverity warnings corrected.
[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.;
d889ce29 109 fPtBinEdges[2][0] = 100.;
110 fPtBinEdges[2][1] = 10.;
b1041e3b 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.;
d889ce29 169 fPtBinEdges[2][0] = 100.;
170 fPtBinEdges[2][1] = 10.;
b1041e3b 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
8da84c2c 655 if(fTrackType==1) {
71e77a79 656 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
8da84c2c 657 if(!track) continue;
658 }
36c36a0c 659 else if(fTrackType==2) {
660 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
8da84c2c 661 if(!track) continue;
662
36c36a0c 663 AliExternalTrackParam exParam;
664 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
665 if( !relate ) {
666 delete track;
667 continue;
668 }
669 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
670 }
d3a3f33d 671 else if(fTrackType==7) {
672 //use global constrained track
42881dab 673 track = new AliESDtrack(*esdtrack);
674 // track = esdtrack;
031aac8e 675 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
d3a3f33d 676 }
71e77a79 677 else
678 track = esdtrack;
83ff93cd 679
680 if(!track) continue;
71e77a79 681
aa3ba8d2 682 if(fTrackType==2) {
683 //Cut on chi2 of constrained fit
684 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
685 delete track;
686 continue;
687 }
688 }
689
fdceab34 690 Int_t label = TMath::Abs(track->GetLabel());
36c36a0c 691 if(label>=nMCtracks) {
692 if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
42881dab 693 if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
36c36a0c 694 continue;
695 }
031aac8e 696
cce400da 697 TParticle *particle = fStack->Particle(label) ;
36c36a0c 698 if(!particle) {
42881dab 699 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
b1041e3b 700 if(track) delete track;
701 }
36c36a0c 702 continue;
703 }
fdceab34 704
71e77a79 705 ptMC = particle->Pt();
706
031aac8e 707 if(fTrackType==1 || fTrackType==2)
0f76d8ae 708 track->GetImpactParametersTPC(dca2D,dcaZ); //TPConly
031aac8e 709 else
710 track->GetImpactParameters(dca2D,dcaZ); //Global
71e77a79 711
fdceab34 712 UChar_t itsMap = track->GetITSClusterMap();
fdceab34 713 for (Int_t i=0; i < 6; i++) {
714 if (itsMap & (1 << i))
715 nPointITS ++;
716 }
71e77a79 717
d889ce29 718 // fPtAll->Fill(pt);
fdceab34 719 fPtAllMC->Fill(ptMC);
71e77a79 720
fdceab34 721 if (fTrackCuts->AcceptTrack(track)) {
722
031aac8e 723 if(fTrackType==7) {
327d12da 724 if(fTrackCutsReject ) {
42881dab 725 if(fTrackCutsReject->AcceptTrack(track) ) {
726 if(track) delete track;
327d12da 727 continue;
42881dab 728 }
327d12da 729 }
730
031aac8e 731 if(esdtrack->GetConstrainedParam())
732 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
031aac8e 733 }
734
735 pt = track->Pt();
736 phi = track->Phi();
737
d889ce29 738 nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
739 chi2C = track->GetConstrainedChi2();
740 relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
741
fdceab34 742 fPtSel->Fill(pt);
36c36a0c 743 if(track->GetLabel()<0) {
744 fPtSelFakes->Fill(pt);
745 fNPointTPCFakes->Fill(track->GetTPCNcls());
746 }
fdceab34 747 fPtSelMC->Fill(ptMC);
57bee263 748 fPtAllvsPtMC->Fill(ptMC,pt);
fdceab34 749 fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
750 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
5a0bd31f 751 fPtAllminPtMCvsPtAllNPointTPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNclsIter1());
752 if(track->GetTPCNcls()>0.)
753 fPtAllminPtMCvsPtAllChi2TPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2()/track->GetTPCNcls());
754 if(track->GetTPCNclsIter1()>0.)
755 fPtAllminPtMCvsPtAllChi2TPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2Iter1()/track->GetTPCNclsIter1());
fdceab34 756 fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
757 fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
758 fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
759 fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
760 fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
761 fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
762 fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
36c36a0c 763
57bee263 764 fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
765 fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
766
36c36a0c 767 //Check if track is reconstructed multiple times
0f76d8ae 768 /*
36c36a0c 769 int multCounter = 1;
770 for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
771 // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
772 AliESDtrack *track2;
773 AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
774 if(!esdtrack2) continue;
775
776 if(fTrackType==1)
777 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
778 else if(fTrackType==2) {
779 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
780 if(!track2) {
36c36a0c 781 continue;
782 }
783 AliExternalTrackParam exParam2;
784 Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
785 if( !relate ) {
786 delete track2;
787 continue;
788 }
789 track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
790 }
791 else
792 track2 = esdtrack2;
793 if(!track2) {
36c36a0c 794 continue;
795 }
796
797 if (fTrackCuts->AcceptTrack(track2)) {
798 Int_t label2 = TMath::Abs(track2->GetLabel());
799 if(label==label2) {
800 fNPointTPCMultRec->Fill(track->GetTPCNcls());
801 fNPointTPCMultRec->Fill(track2->GetTPCNcls());
802 fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
803 multCounter++;
804 }
805 }
806 if(fTrackType==1 || fTrackType==2) delete track2;
807 }//track2 loop
808 if(multCounter>1) fMultRec->Fill(multCounter);
0f76d8ae 809 */
36c36a0c 810
811 }//fTrackCuts selection
812
5a0bd31f 813
42881dab 814 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
5a0bd31f 815 if(track) delete track;
816 }
36c36a0c 817
fdceab34 818 }//ESD track loop
819
820 // Post output data
821 PostData(0, fHistList);
5a0bd31f 822
fdceab34 823}
b4691ee7 824//________________________________________________________________________
825Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
826 //
827 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
828 // This is to called in Notify and should provide the path to the AOD/ESD file
829 // Copied from AliAnalysisTaskJetSpectrum2
830 //
831
832 TString file(currFile);
833 fXsec = 0;
834 fTrials = 1;
835
836 if(file.Contains("root_archive.zip#")){
837 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
838 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
839 file.Replace(pos+1,20,"");
840 }
841 else {
842 // not an archive take the basename....
843 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
844 }
cb76764e 845 // Printf("%s",file.Data());
b4691ee7 846
847
848 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
849 if(!fxsec){
850 // next trial fetch the histgram file
851 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
852 if(!fxsec){
853 // not a severe condition but inciate that we have no information
854 return kFALSE;
855 }
856 else{
857 // find the tlist we want to be independtent of the name so use the Tkey
858 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
859 if(!key){
860 fxsec->Close();
861 return kFALSE;
862 }
863 TList *list = dynamic_cast<TList*>(key->ReadObj());
864 if(!list){
865 fxsec->Close();
866 return kFALSE;
867 }
868 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
869 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
870 fxsec->Close();
871 }
872 } // no tree pyxsec.root
873 else {
874 TTree *xtree = (TTree*)fxsec->Get("Xsection");
875 if(!xtree){
876 fxsec->Close();
877 return kFALSE;
878 }
879 UInt_t ntrials = 0;
880 Double_t xsection = 0;
881 xtree->SetBranchAddress("xsection",&xsection);
882 xtree->SetBranchAddress("ntrials",&ntrials);
883 xtree->GetEntry(0);
884 fTrials = ntrials;
885 fXsec = xsection;
886 fxsec->Close();
887 }
888 return kTRUE;
889}
890//________________________________________________________________________
891Bool_t AliPWG4HighPtQAMC::Notify()
892{
893 //
894 // Implemented Notify() to read the cross sections
895 // and number of trials from pyxsec.root
896 // Copied from AliAnalysisTaskJetSpectrum2
897 //
898
899 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
900 Float_t xsection = 0;
901 Float_t ftrials = 1;
902
903 fAvgTrials = 1;
904 if(tree){
905 TFile *curfile = tree->GetCurrentFile();
906 if (!curfile) {
907 Error("Notify","No current file");
908 return kFALSE;
909 }
910 if(!fh1Xsec||!fh1Trials){
cb76764e 911 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
b4691ee7 912 return kFALSE;
913 }
914 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
915 fh1Xsec->Fill("<#sigma>",xsection);
916 // construct a poor man average trials
917 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
918 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
919 }
920 return kTRUE;
921}
922
923//________________________________________________________________________
924AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
925
926 if(!mcEvent)return 0;
927 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
928 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
929 if(!pythiaGenHeader){
930 // cocktail ??
931 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
932
933 if (!genCocktailHeader) {
cb76764e 934 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
b4691ee7 935 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
936 return 0;
937 }
938 TList* headerList = genCocktailHeader->GetHeaders();
939 for (Int_t i=0; i<headerList->GetEntries(); i++) {
940 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
941 if (pythiaGenHeader)
942 break;
943 }
944 if(!pythiaGenHeader){
945 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
946 return 0;
947 }
948 }
949 return pythiaGenHeader;
950
951}
952
fdceab34 953//________________________________________________________________________
954void AliPWG4HighPtQAMC::Terminate(Option_t *)
955{
df943115 956 // The Terminate() function is the last function to be called during
957 // a query. It always runs on the client, it can be used to present
958 // the results graphically or save the results to file.
fdceab34 959
fdceab34 960}
df943115 961
962#endif