1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------------
24 #ifndef ALIPWG4HighPtQAMC_CXX
25 #define ALIPWG4HighPtQAMC_CXX
27 #include "AliPWG4HighPtQAMC.h"
40 #include "AliAnalysisTask.h"
41 #include "AliAnalysisManager.h"
42 #include "AliESDInputHandler.h"
43 #include "AliMCEvent.h"
44 #include "AliMCEventHandler.h"
46 #include "AliESDtrack.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliExternalTrackParam.h"
50 #include "AliGenPythiaEventHeader.h"
51 #include "AliGenCocktailEventHeader.h"
52 //#include "AliAnalysisHelperJetTasks.h"
54 using namespace std; //required for resolving the 'cout' symbol
56 ClassImp(AliPWG4HighPtQAMC)
58 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
59 : AliAnalysisTask("AliPWG4HighPtQAMC", ""),
67 fSigmaConstrainedMax(1e6),
86 fPtAllminPtMCvsPtAll(0),
87 fPtAllvsPtMCvsMult(0),
88 fPtAllminPtMCvsPtAllvsMult(0),
89 fPtAllminPtMCvsPtAllNPointTPC(0),
90 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
91 fPtAllminPtMCvsPtAllChi2TPC(0),
92 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
93 fPtAllminPtMCvsPtAllDCAR(0),
94 fPtAllminPtMCvsPtAllDCAZ(0),
95 fPtAllminPtMCvsPtAllPhi(0),
96 fPtAllminPtMCvsPtAllNPointITS(0),
97 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
98 fPtAllminPtMCvsPtAllChi2C(0),
99 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
105 fPtBinEdges[0][0] = 10.;
106 fPtBinEdges[0][1] = 1.;
107 fPtBinEdges[1][0] = 20.;
108 fPtBinEdges[1][1] = 2.;
109 fPtBinEdges[2][0] = 100.;
110 fPtBinEdges[2][1] = 10.;
113 //________________________________________________________________________
114 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
115 AliAnalysisTask(name,""),
123 fSigmaConstrainedMax(1e6),
139 fNPointTPCMultRec(0),
142 fPtAllminPtMCvsPtAll(0),
143 fPtAllvsPtMCvsMult(0),
144 fPtAllminPtMCvsPtAllvsMult(0),
145 fPtAllminPtMCvsPtAllNPointTPC(0),
146 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
147 fPtAllminPtMCvsPtAllChi2TPC(0),
148 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
149 fPtAllminPtMCvsPtAllDCAR(0),
150 fPtAllminPtMCvsPtAllDCAZ(0),
151 fPtAllminPtMCvsPtAllPhi(0),
152 fPtAllminPtMCvsPtAllNPointITS(0),
153 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
154 fPtAllminPtMCvsPtAllChi2C(0),
155 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
161 // Constructor. Initialization of Inputs and Outputs
163 AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
165 fPtBinEdges[0][0] = 10.;
166 fPtBinEdges[0][1] = 1.;
167 fPtBinEdges[1][0] = 20.;
168 fPtBinEdges[1][1] = 2.;
169 fPtBinEdges[2][0] = 100.;
170 fPtBinEdges[2][1] = 10.;
172 // Input slot #0 works with a TChain ESD
173 DefineInput(0, TChain::Class());
174 // Output slot #0, #1 write into a TList
175 DefineOutput(0, TList::Class());
178 //________________________________________________________________________
179 void AliPWG4HighPtQAMC::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
182 fPtBinEdges[region][0] = ptmax;
183 fPtBinEdges[region][1] = ptBinWidth;
186 AliError("Only 3 regions alowed. Use region 0/1/2\n");
192 //________________________________________________________________________
193 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
195 // Connect ESD and MC event handler here
197 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
198 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
200 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
204 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
207 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
210 fESD = esdH->GetEvent();
212 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
214 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
217 fMC = eventHandler->MCEvent();
222 //________________________________________________________________________
223 void AliPWG4HighPtQAMC::CreateOutputObjects() {
224 //Create output objects
225 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
227 Bool_t oldStatus = TH1::AddDirectoryStatus();
228 TH1::AddDirectory(kFALSE);
231 fHistList = new TList();
232 fHistList->SetOwner(kTRUE);
234 Float_t fgkPtMin = 0.;
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) ;
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 ;
261 Int_t fgkResPtBins=80;
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 ;
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 ;
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 ;
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 ;
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) {
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 ;
295 Int_t fgkNDCAZBins=80;
296 Float_t fgkDCAZMin = -2.;
297 Float_t fgkDCAZMax = 2.;
298 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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);
340 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
342 fNEventReject->Fill("noESD",0);
343 fNEventReject->Fill("Trigger",0);
344 fNEventReject->Fill("noMCEvent",0);
345 fNEventReject->Fill("noStack",0);
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);
351 fHistList->Add(fNEventReject);
353 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
354 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
355 fHistList->Add(fh1Xsec);
357 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
358 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
359 fHistList->Add(fh1Trials);
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);
366 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
367 fHistList->Add(fPtAll);
368 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
369 fHistList->Add(fPtSel);
371 fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, binsPt);
372 fHistList->Add(fPtSelFakes);
374 fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",fgkNNClustersTPCBins, binsNClustersTPC);
375 fHistList->Add(fNPointTPCFakes);
377 fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, binsPt);
378 fHistList->Add(fPtSelLargeLabel);
380 fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",fgkMultBins, binsMult);
381 fHistList->Add(fMultRec);
383 fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",fgkNNClustersTPCBins, binsNClustersTPC);
384 fHistList->Add(fNPointTPCMultRec);
386 fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",50,0.,50.,40,-20.,20.);
387 fHistList->Add(fDeltaPtMultRec);
389 fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, binsPt,fgkNPtBins, binsPt);
390 fHistList->Add(fPtAllvsPtMC);
392 fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
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);
397 fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, binsPt,fgkNPtBins, binsPt,fgkMultBins,binsMult);
398 fHistList->Add(fPtAllvsPtMCvsMult);
400 fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkMultBins,binsMult);
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);
406 fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
407 fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
408 fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
409 fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{cls,TPC}");
410 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
412 fPtAllminPtMCvsPtAllNPointTPCIter1 = new TH3F("fPtAllminPtMCvsPtAllNPointTPCIter1","PtAllminPtMCvsPtAllNPointTPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
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);
418 fPtAllminPtMCvsPtAllChi2TPC = new TH3F("fPtAllminPtMCvsPtAllChi2TPC","PtAllminPtMCvsPtAllChi2TPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
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);
424 fPtAllminPtMCvsPtAllChi2TPCIter1 = new TH3F("fPtAllminPtMCvsPtAllChi2TPCIter1","PtAllminPtMCvsPtAllChi2TPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
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);
430 fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCA2DBins,binsDCA2D);
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);
436 fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCAZBins,binsDCAZ);
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);
442 fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNPhiBins,binsPhi);
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);
448 fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNPointITSBins,binsNPointITS);
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);
454 fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
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);
460 fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2CBins,binsChi2C);
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);
466 fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
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);
473 fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, binsPt);
474 fHistList->Add(fPtAllMC);
475 fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, binsPt);
476 fHistList->Add(fPtSelMC);
478 TH1::AddDirectory(oldStatus);
480 PostData(0, fHistList);
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;
498 //________________________________________________________________________
499 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
501 // Decide if event should be selected for analysis
504 // Checks following requirements:
506 // - trigger info from AliPhysicsSelection
507 // - MCevent available
508 // - number of reconstructed tracks > 1
509 // - primary vertex reconstructed
510 // - z-vertex < 10 cm
512 Bool_t selectEvent = kTRUE;
514 //fESD object available?
516 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
517 fNEventReject->Fill("noESD",1);
518 selectEvent = kFALSE;
523 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
524 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
525 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
526 fNEventReject->Fill("Trigger",1);
527 selectEvent = kFALSE;
534 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
535 fStack = fMC->Stack(); //Particles Stack
537 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
540 AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
541 fNEventReject->Fill("noStack",1);
542 selectEvent = kFALSE;
546 AliDebug(2,Form("ERROR: Could not retrieve stack"));
547 fNEventReject->Fill("noMCEvent",1);
548 selectEvent = kFALSE;
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;
559 //Check if vertex is reconstructed
560 if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
561 else fVtx = fESD->GetPrimaryVertexSPD();
564 fNEventReject->Fill("noVTX",1);
565 selectEvent = kFALSE;
569 if(!fVtx->GetStatus()) {
570 fNEventReject->Fill("VtxStatus",1);
571 selectEvent = kFALSE;
576 // TString vtxName(fVtx->GetName());
577 if(fVtx->GetNContributors()<2) {
578 fNEventReject->Fill("NCont<2",1);
579 selectEvent = kFALSE;
583 //Check if z-vertex < 10 cm
585 fVtx->GetXYZ(primVtx);
586 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
587 fNEventReject->Fill("ZVTX>10",1);
588 selectEvent = kFALSE;
592 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
598 //________________________________________________________________________
599 void AliPWG4HighPtQAMC::Exec(Option_t *) {
601 // Called for each event
602 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
603 // All events without selection
604 fNEventAll->Fill(0.);
608 PostData(0, fHistList);
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
617 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
619 nTrials = pythiaGenHeader->Trials();
620 ptHard = pythiaGenHeader->GetPtHard();
622 fh1PtHard->Fill(ptHard);
623 fh1PtHardTrials->Fill(ptHard,nTrials);
625 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
629 //Need to keep track of selected events
630 fNEventSel->Fill(0.);
632 Int_t nTracks = fESD->GetNumberOfTracks();
633 AliDebug(2,Form("nTracks ESD%d", nTracks));
635 int nMCtracks = fStack->GetNtrack();
644 Float_t nSigmaToVertex = 0.;
645 Float_t relUncertainty1Pt = 0.;
647 int mult = fTrackCuts->CountAcceptedTracks(fESD);
649 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
651 AliESDtrack *track = 0;
652 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
653 if(!esdtrack) continue;
656 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
659 else if(fTrackType==2) {
660 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
663 AliExternalTrackParam exParam;
664 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
669 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
671 else if(fTrackType==7) {
672 //use global constrained track
673 track = new AliESDtrack(*esdtrack);
675 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
683 //Cut on chi2 of constrained fit
684 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
690 Int_t label = TMath::Abs(track->GetLabel());
691 if(label>=nMCtracks) {
692 if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
693 if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
697 TParticle *particle = fStack->Particle(label) ;
699 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
700 if(track) delete track;
705 ptMC = particle->Pt();
707 if(fTrackType==1 || fTrackType==2)
708 track->GetImpactParametersTPC(dca2D,dcaZ); //TPConly
710 track->GetImpactParameters(dca2D,dcaZ); //Global
712 UChar_t itsMap = track->GetITSClusterMap();
713 for (Int_t i=0; i < 6; i++) {
714 if (itsMap & (1 << i))
719 fPtAllMC->Fill(ptMC);
721 if (fTrackCuts->AcceptTrack(track)) {
724 if(fTrackCutsReject ) {
725 if(fTrackCutsReject->AcceptTrack(track) ) {
726 if(track) delete track;
731 if(esdtrack->GetConstrainedParam())
732 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
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;
743 if(track->GetLabel()<0) {
744 fPtSelFakes->Fill(pt);
745 fNPointTPCFakes->Fill(track->GetTPCNcls());
747 fPtSelMC->Fill(ptMC);
748 fPtAllvsPtMC->Fill(ptMC,pt);
749 fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
750 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
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());
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);
764 fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
765 fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
767 //Check if track is reconstructed multiple times
770 for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
771 // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
773 AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
774 if(!esdtrack2) continue;
777 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
778 else if(fTrackType==2) {
779 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
783 AliExternalTrackParam exParam2;
784 Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
789 track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
797 if (fTrackCuts->AcceptTrack(track2)) {
798 Int_t label2 = TMath::Abs(track2->GetLabel());
800 fNPointTPCMultRec->Fill(track->GetTPCNcls());
801 fNPointTPCMultRec->Fill(track2->GetTPCNcls());
802 fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
806 if(fTrackType==1 || fTrackType==2) delete track2;
808 if(multCounter>1) fMultRec->Fill(multCounter);
811 }//fTrackCuts selection
814 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
815 if(track) delete track;
821 PostData(0, fHistList);
824 //________________________________________________________________________
825 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
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
832 TString file(currFile);
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,"");
842 // not an archive take the basename....
843 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
845 // Printf("%s",file.Data());
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
850 // next trial fetch the histgram file
851 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
853 // not a severe condition but inciate that we have no information
857 // find the tlist we want to be independtent of the name so use the Tkey
858 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
863 TList *list = dynamic_cast<TList*>(key->ReadObj());
868 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
869 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
872 } // no tree pyxsec.root
874 TTree *xtree = (TTree*)fxsec->Get("Xsection");
880 Double_t xsection = 0;
881 xtree->SetBranchAddress("xsection",&xsection);
882 xtree->SetBranchAddress("ntrials",&ntrials);
890 //________________________________________________________________________
891 Bool_t AliPWG4HighPtQAMC::Notify()
894 // Implemented Notify() to read the cross sections
895 // and number of trials from pyxsec.root
896 // Copied from AliAnalysisTaskJetSpectrum2
899 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
900 Float_t xsection = 0;
905 TFile *curfile = tree->GetCurrentFile();
907 Error("Notify","No current file");
910 if(!fh1Xsec||!fh1Trials){
911 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
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;
923 //________________________________________________________________________
924 AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
926 if(!mcEvent)return 0;
927 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
928 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
929 if(!pythiaGenHeader){
931 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
933 if (!genCocktailHeader) {
934 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
935 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
938 TList* headerList = genCocktailHeader->GetHeaders();
939 for (Int_t i=0; i<headerList->GetEntries(); i++) {
940 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
944 if(!pythiaGenHeader){
945 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
949 return pythiaGenHeader;
953 //________________________________________________________________________
954 void AliPWG4HighPtQAMC::Terminate(Option_t *)
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.