]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/AliAnalysisTaskChargedHadronSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fixes for BUFFER_SIZE
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisTaskChargedHadronSpectra.cxx
CommitLineData
635c8b48 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// //
18// //
19// Analysis for identified particle spectra measured with TPC dE/dx. //
20// //
21// //
22///////////////////////////////////////////////////////////////////////////
23
24#include "Riostream.h"
25#include "TChain.h"
26#include "TTree.h"
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TH3F.h"
30#include "TList.h"
31#include "TMath.h"
32#include "TCanvas.h"
33#include "TObjArray.h"
34#include "TF1.h"
35#include "TFile.h"
36
37#include "AliAnalysisTaskSE.h"
38#include "AliAnalysisManager.h"
39
40#include "AliHeader.h"
41#include "AliGenPythiaEventHeader.h"
42#include "AliGenCocktailEventHeader.h"
43
44#include "AliPID.h"
45#include "AliESDtrackCuts.h"
46#include "AliESDVertex.h"
47#include "AliESDEvent.h"
48#include "AliESDInputHandler.h"
49#include "AliESDtrack.h"
10d100d4 50#include "AliESDpid.h"
635c8b48 51
52#include "AliMCEventHandler.h"
53#include "AliMCEvent.h"
54#include "AliStack.h"
55
56#include "AliLog.h"
57
58#include "AliAnalysisTaskChargedHadronSpectra.h"
59
60
61ClassImp(AliAnalysisTaskChargedHadronSpectra)
62
63//________________________________________________________________________
64AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra()
3c038b30 65: AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
66 fMCtrue(0),
67 fTrackingMode(0),
68 fAlephParameters(),
69 fHistPtMCKaon(0),
70 fHistPtMCProton(0),
71 fHistPtMCPion(0),
72 fHistPtEtaKaon(0),
73 fHistPtEtaKaonNoKink(0),
74 fHistPtEtaProton(0),
75 fHistPtEtaProtonDCA(0),
76 fHistPtEtaPion(0),
77 fHistClassicKaon(0),
78 fHistClassicProton(0),
79 fHistClassicPion(0),
80 fDeDx(0),
81 fHistTrackPerEvent(0),
82 fHistTrackPerEventMC(0),
83 fSecProtons(0),
84 fVertexZ(0),
85 fHistEtaNcls(0),
86 fHistEtaPhi(0),
87 fHistEffProton(0),
88 fHistEffProtonDCA(0),
89 fHistEffPion(0),
90 fHistEffKaon(0),
91 fHistRealTracks(0),
92 fHistMCparticles(0)
635c8b48 93{
94 // default Constructor
95}
96
97
98//________________________________________________________________________
99AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra(const char *name)
10d100d4 100 : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
635c8b48 101 fMCtrue(0),
3c038b30 102 fTrackingMode(0),
635c8b48 103 fAlephParameters(),
104 fHistPtMCKaon(0),
105 fHistPtMCProton(0),
106 fHistPtMCPion(0),
635c8b48 107 fHistPtEtaKaon(0),
108 fHistPtEtaKaonNoKink(0),
109 fHistPtEtaProton(0),
110 fHistPtEtaProtonDCA(0),
111 fHistPtEtaPion(0),
635c8b48 112 fHistClassicKaon(0),
113 fHistClassicProton(0),
114 fHistClassicPion(0),
635c8b48 115 fDeDx(0),
116 fHistTrackPerEvent(0),
117 fHistTrackPerEventMC(0),
118 fSecProtons(0),
119 fVertexZ(0),
120 fHistEtaNcls(0),
121 fHistEtaPhi(0),
122 fHistEffProton(0),
123 fHistEffProtonDCA(0),
124 fHistEffPion(0),
125 fHistEffKaon(0),
3c038b30 126 fHistRealTracks(0),
127 fHistMCparticles(0)
635c8b48 128{
129 //
3c038b30 130 // standard constructur which should be used
635c8b48 131 //
635c8b48 132 Printf("*** CONSTRUCTOR CALLED ****");
3c038b30 133 //
134 fMCtrue = kFALSE; // change three things for processing real data: 1. set this to false! / 2. change ALEPH parameters / 3. change parameters in the Exec()
135 fTrackingMode = 0;
136 /* real */
137 fAlephParameters[0] = 0.0283086;
138 fAlephParameters[1] = 2.63394e+01;
139 fAlephParameters[2] = 5.04114e-11;
140 fAlephParameters[3] = 2.12543e+00;
141 fAlephParameters[4] = 4.88663e+00;
142 //
143 // initialize PID object
144 //
145 fESDpid = new AliESDpid();
146 //
147 // create track cuts
148 //
149 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
150 //
151 Initialize();
635c8b48 152 // Output slot #0 writes into a TList container
c77fc043 153 DefineOutput(1, TList::Class());
635c8b48 154
3c038b30 155}
156
157
158//________________________________________________________________________
159void AliAnalysisTaskChargedHadronSpectra::Initialize()
160{
161 //
162 // updating parameters in case of changes
163 //
164 // 1. pid parameters
165 //
166 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
167 //
168 // 2. track cuts
169 //
170 //fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
171 //fESDtrackCuts->SetMaxNsigmaToVertex(3);
172 //fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
173 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
174 fESDtrackCuts->SetMinNClustersTPC(80);
175 fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
176 //fESDtrackCuts->SetMaxDCAToVertexXY(3);
177 //fESDtrackCuts->SetMaxDCAToVertexZ(3);
178 //
179 if (fTrackingMode == 1) {//for GLOBAL running IN ADDITION
180 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
181 fESDtrackCuts->SetRequireITSRefit(kTRUE);
182 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
183 //
184 //fESDtrackCuts->SetMinNClustersITS(3);
185 //
186 //fESDtrackCuts->SetMaxDCAToVertexXY(0.5);
187 //fESDtrackCuts->SetMaxDCAToVertexZ(0.05); // this is now done in this special function...
188 }
189
635c8b48 190
191}
192
635c8b48 193
194//________________________________________________________________________
c77fc043 195void AliAnalysisTaskChargedHadronSpectra::UserCreateOutputObjects()
635c8b48 196{
197 // Create histograms
198 // Called once
199 fListHist = new TList();
c77fc043 200 //fListHist->SetOwner(); // Whoever knows how the data handling is ...?
635c8b48 201
3c038b30 202 const Int_t kPtBins = 2*30;
203 const Double_t kPtMax = 2.0;
635c8b48 204 const Int_t kEtaBins = 4;
3c038b30 205 const Double_t kEtaMax = 1;
206 const Int_t kDeDxBins = 100;
635c8b48 207 const Double_t kDeDxMax = 1;
208 const Int_t kMultBins = 10;
209 const Int_t kMultMax = 300;
210
211 // sort pT-bins ..
3c038b30 212 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0};
635c8b48 213 Int_t indexes[kPtBins+1];
214 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
215 Double_t binsPt[kPtBins+1];
216 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
217
218
219 // MC histograms
220 fHistPtMCKaon = new TH3F("HistPtMCKaon", "PtEtaKaon; mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
221 fHistPtMCProton = new TH3F("HistPtMCProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
222 fHistPtMCPion = new TH3F("HistPtMCPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
635c8b48 223 //
224 fHistPtMCKaon->GetZaxis()->Set(kPtBins, binsPt);
225 fHistPtMCProton->GetZaxis()->Set(kPtBins, binsPt);
226 fHistPtMCPion->GetZaxis()->Set(kPtBins, binsPt);
635c8b48 227 //
228 fListHist->Add(fHistPtMCKaon);
229 fListHist->Add(fHistPtMCProton);
230 fListHist->Add(fHistPtMCPion);
3c038b30 231 //
635c8b48 232 // reconstructed particle histograms
233 fHistPtEtaKaon = new TH3F("HistPtEtaKaon", "PtEtaKaon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
234 fHistPtEtaKaonNoKink = new TH3F("HistPtEtaKaonNoKink", "PtEtaKaonNoKink;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
235 fHistPtEtaProton = new TH3F("HistPtEtaProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
3c038b30 236 fHistPtEtaProtonDCA = new TH3F("HistPtEtaProtonDCA", "PtEtaProton;DCA (cm); #eta; p_{T} (GeV)",1000,0,50,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
635c8b48 237 fHistPtEtaPion = new TH3F("HistPtEtaPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
635c8b48 238 //
239 fHistPtEtaKaon->GetZaxis()->Set(kPtBins, binsPt);
240 fHistPtEtaKaonNoKink->GetZaxis()->Set(kPtBins, binsPt);
241 fHistPtEtaProton->GetZaxis()->Set(kPtBins, binsPt);
242 fHistPtEtaProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
243 fHistPtEtaPion->GetZaxis()->Set(kPtBins, binsPt);
635c8b48 244 //
245 fListHist->Add(fHistPtEtaKaon);
246 fListHist->Add(fHistPtEtaKaonNoKink);
247 fListHist->Add(fHistPtEtaProton);
248 fListHist->Add(fHistPtEtaProtonDCA);
249 fListHist->Add(fHistPtEtaPion);
3c038b30 250
635c8b48 251 // histograms for the classical analysis
252 fHistClassicKaon = new TH3F("HistClassicKaon", "PtEtaKaon;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
253 fHistClassicProton = new TH3F("HistClassicProton", "PtEtaProton;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
254 fHistClassicPion = new TH3F("HistClassicPion", "PtEtaPion;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
635c8b48 255 //
256 fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
257 fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
258 fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
635c8b48 259 //
260 fListHist->Add(fHistClassicKaon);
261 fListHist->Add(fHistClassicProton);
262 fListHist->Add(fHistClassicPion);
635c8b48 263 // histograms of general interest
3c038b30 264 fDeDx = new TH3F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.1,100.,1000,0.,1000,2,-2,2);
635c8b48 265 BinLogX(fDeDx);
266 fListHist->Add(fDeDx);
3c038b30 267 fHistTrackPerEvent = new TH2F("HistTrackPerEvent", "Tracks per event;Number of Tracks;Number of Events; code",101,-0.5,100,10,-0.5,9.5);
635c8b48 268 fListHist->Add(fHistTrackPerEvent);
3c038b30 269 fHistTrackPerEventMC = new TH3F("HistTrackPerEventMC", "Tracks per event;code;TrackPerEvent;IsSelected",10,-0.5,9,101,-0.5,100,2,-0.5,1.5);
635c8b48 270 fListHist->Add(fHistTrackPerEventMC);
271 fSecProtons = new TH2F("SecProtons", "xy;x;y",100,-10,10,100,-10,10);
272 fListHist->Add(fSecProtons);
3c038b30 273 fVertexZ = new TH3F("VertexZ", "vertex position;x (cm); y (cm); z (cm); counts",100,-5,5,100,-5,5,100,-25,25);
635c8b48 274 fListHist->Add(fVertexZ);
275 //
276 fHistEtaNcls = new TH2F("HistEtaNcls", "EtaNcls;#eta;Ncls",100,-1.5,1.5,80,0,160);
277 fListHist->Add(fHistEtaNcls);
278 fHistEtaPhi = new TH2F("HistEtaPhi", "EtaNcls;#eta;#phi",100,-1.5,1.5,100,0,2*TMath::Pi());
279 fListHist->Add(fHistEtaPhi);
280
281 // histograms for a refined efficiency study
282 fHistEffProton = new TH3F("HistEffProton", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
283 fListHist->Add(fHistEffProton);
284 fHistEffProtonDCA = new TH3F("HistEffProtonDCA", "Eff;p_{T} (GeV); code",10,-0.5,9.5,200,0,15,kPtBins,-kPtMax,kPtMax);
285 fListHist->Add(fHistEffProtonDCA);
286 fHistEffPion = new TH3F("HistEffPion", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
287 fListHist->Add(fHistEffPion);
288 fHistEffKaon = new TH3F("HistEffKaon", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
289 fListHist->Add(fHistEffKaon);
290 //
291 fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
292 fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
293 fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
294 fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
295 //
3c038b30 296 // create the histograms with all necessary information
635c8b48 297 //
3c038b30 298 // (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rapPion, (5.) rapKaon, (6.) rapProton, (7.) dca, (8.) nclIter1, (9.) pullPion, (10.) pullKaon, (11.) pullProton
299 // 0, 1, 2, 3, 4 , 5, 6, 7, 8, 9, 10, 11
300 Int_t binsHistReal[12] = { 301, 16,20, 2, 16, 16, 16,60, 10,100,100,100};
301 Double_t xminHistReal[12] = { -0.5,-0.8, 0,-2,-0.8,-0.8,-0.8, 0, 60,-10,-10,-10};
302 Double_t xmaxHistReal[12] = {300.5, 0.8, 2, 2, 0.8, 0.8, 0.8,15,160, 10, 10, 10};
303 fHistRealTracks = new THnSparseS("fHistRealTracks","real tracks; (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rapPion, (5.) rapKaon, (6.) rapProton, (7.) dca, (8.) nclIter1, (9.) pullPion, (10.) pullKaon, (11.) pullProton",12,binsHistReal,xminHistReal,xmaxHistReal);
304 fListHist->Add(fHistRealTracks);
305 //
306 // (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rap, (5.) dca, (6.) ID, (7.) code
307 // 0, 1, 2, 3, 4, 5, 6, 7
308 Int_t binsHistMC[8] = { 301, 16,20, 2, 16,60, 5, 10};
309 Double_t xminHistMC[8] = { -0.5,-0.8, 0,-2,-0.8, 0,-0.5,-0.5};
310 Double_t xmaxHistMC[8] = {300.5, 0.8, 2, 2, 0.8,15, 4.5, 9.5};
311 fHistMCparticles = new THnSparseS("fHistMCparticles"," (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rap, (5.) dca, (6.) ID, (7.) code",8,binsHistMC,xminHistMC,xmaxHistMC);
312 fListHist->Add(fHistMCparticles);
635c8b48 313 //
635c8b48 314 //
635c8b48 315 //
635c8b48 316
3c038b30 317
635c8b48 318}
319
320//________________________________________________________________________
c77fc043 321void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *)
635c8b48 322{
3c038b30 323 //
635c8b48 324 // main event loop
3c038b30 325 //
635c8b48 326 AliLog::SetGlobalLogLevel(AliLog::kError);
327 //
328 // Check Monte Carlo information and other access first:
329 //
330 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
331 if (!eventHandler) {
332 //Printf("ERROR: Could not retrieve MC event handler");
333 if (fMCtrue) return;
334 }
3c038b30 335 //
336 AliMCEvent* mcEvent = 0x0;
337 AliStack* stack = 0x0;
338 if (eventHandler) mcEvent = eventHandler->MCEvent();
635c8b48 339 if (!mcEvent) {
340 //Printf("ERROR: Could not retrieve MC event");
341 if (fMCtrue) return;
342 }
3c038b30 343 if (fMCtrue) {
344 stack = mcEvent->Stack();
345 if (!stack) return;
346 }
347 //
c77fc043 348 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
635c8b48 349 if (!fESD) {
350 //Printf("ERROR: fESD not available");
351 return;
352 }
353
354 if (!fESDtrackCuts) {
355 Printf("ERROR: fESDtrackCuts not available");
356 return;
357 }
358
359 Int_t trackCounter = 0;
3c038b30 360 //
361 // check if event is selected by physics selection class
362 //
363 fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks(),0);
364 Bool_t isSelected = kFALSE;
365 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
366 if (isSelected) fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks(),1);
367 //
368 // monitor vertex position
369 //
370 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
371 if(vertex->GetNContributors()<1) {
372 // SPD vertex
373 vertex = fESD->GetPrimaryVertexSPD();
374 if(vertex->GetNContributors()<1) vertex = 0x0;
375 }
376
377 // 1st track loop to determine multiplicities
378 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
379 AliESDtrack *track = 0x0;
380 if (fTrackingMode == 0) track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
381 if (fTrackingMode == 1) track = new AliESDtrack(*fESD->GetTrack(i)); // alternative: global track GLOBAL GLOBAL GLOBAL GLOBAL GLOBAL GLOBAL GLOBAL
382 if (!track) continue;
383 if (fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta())< 0.9) {
384 trackCounter++;
385 //
386 }
387 delete track;
388 }
389
390 fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
391 fHistTrackPerEvent->Fill(trackCounter,2);
392
393 if (fMCtrue) {
394 //
395 //
396 //
397 AliHeader * header = mcEvent->Header();
398 Int_t processtype = GetPythiaEventProcessType(header);
399 // non diffractive
400 if (processtype !=92 && processtype !=93 && processtype != 94) fHistTrackPerEventMC->Fill(1., trackCounter, isSelected);
401 // single diffractive
402 if ((processtype == 92 || processtype == 93)) fHistTrackPerEventMC->Fill(2., trackCounter, isSelected);
403 // double diffractive
404 if (processtype == 94) fHistTrackPerEventMC->Fill(3., trackCounter, isSelected);
405 //
406 Int_t mult = trackCounter;//stack->GetNprimary();
407 //
408 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
409 TParticle * trackMC = stack->Particle(i);
410 Int_t pdg = trackMC->GetPdgCode();
411 //
412 Double_t xv = trackMC->Vx();
413 Double_t yv = trackMC->Vy();
414 Double_t zv = trackMC->Vz();
415 Double_t dxy = 0;
416 dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
417 Double_t dz = 0;
418 dz = TMath::Abs(zv); // so stupid to avoid warnings
419 //
420 // vertex cut - selection of primaries
421 //
422 //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
423 //
424 if (!stack->IsPhysicalPrimary(i)) continue;
425 //
426 if (pdg == 321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select K+ only
427 if (pdg == 211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select Pi+ only
428 if (pdg == 2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select p+ only
429 //
430 if (pdg == -321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select K- only
431 if (pdg == -211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select Pi- only
432 if (pdg == -2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select p- only
433 //
434 // special Kaon checks - those which decay before entering the TPC fiducial volume
435 // <-> should be revisited by looping from first to last daugther and check if UniqueID is 4
436 //
437 if (TMath::Abs(pdg)==321) {
438 TParticle * trackDaughterMC = stack->Particle(TMath::Abs(trackMC->GetFirstDaughter()));
439 Int_t pdgDaughter = trackDaughterMC->GetPdgCode();
440 if (TMath::Abs(pdgDaughter) == 13 && pdg == 321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Y()), trackMC->Pt()); // Kaon control hist
441 if (TMath::Abs(pdgDaughter) == 13 && pdg == -321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // Kaon control hist
442 }
443 }
444 }
445 //
446 //end MC tracks loop <-> now check event selection criteria
447 //
448 if (!isSelected) {
449 PostData(1, fListHist);
450 return;
451 }
452 //
453 if (!vertex) {
454 fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
455 PostData(1, fListHist);
456 return;
457 } else {
458 fVertexZ->Fill(vertex->GetXv(),vertex->GetYv(),vertex->GetZv());
459 if (TMath::Abs(vertex->GetZv()) > 10) {
460 fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
461 PostData(1, fListHist);
462 return;
463 }
464 }
465 //
466 // tarck loop
467 //
468 const Float_t kNsigmaCut = 3;
469 const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
470
471 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
472
473 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
474 //
475 // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
476 //
477 AliESDtrack *track = 0x0;
478 if (fTrackingMode == 0) track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
479 if (fTrackingMode == 1) track = new AliESDtrack(*fESD->GetTrack(i)); // alternative: global track GLOBAL GLOBAL GLOBAL GLOBAL
480 if (!track) continue;
481 AliESDtrack *trackDummy = fESD->GetTrack(i);
482 if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
483 delete track;
484 continue;
485 }
486 //
487 if (track->GetTPCNclsIter1() < 80) {
488 delete track;
489 continue;
490 }
491 //
492 AliExternalTrackParam trackIn(*trackDummy->GetInnerParam());
493 Double_t ptot = trackIn.GetP(); // momentum for dEdx determination
494 Double_t pT = track->Pt();
495 Double_t eta = TMath::Abs(track->Eta());
496 //
497 fHistEtaNcls->Fill(track->Eta(),track->GetTPCNcls());
498 fHistEtaPhi->Fill(track->Eta(),track->Phi());
499 //
500 // cut for dead regions in the detector
501 // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
502 //
503 // 2.a) apply some standard track cuts according to general recommendations
504 //
505 if (!fESDtrackCuts->AcceptTrack(track)) {
506 //
507 if (fMCtrue) {
508 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
509 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
510 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(6,eta,track->GetSign()*pT); // Kaon control hist
511 }
512 //
513 delete track;
514 continue;
515 }
516
517 if (fTrackingMode == 1) {
518 /*if (!SelectOnImpPar(track)) {
519 delete track;
520 continue;
521 }*/
522 }
523
524 //
525 // calculate rapidities and kinematics
526 //
527 //
528 Double_t pvec[3];
529 track->GetPxPyPz(pvec);
530 Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
531 Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
532 Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
533 //
534 Double_t rapPion = TMath::Abs(0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2])));
535 Double_t rapKaon = TMath::Abs(0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2])));
536 Double_t rapProton = TMath::Abs(0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2])));
537 //
538 // 3. make the PID
539 //
540 Double_t sign = track->GetSign();
541 Double_t tpcSignal = track->GetTPCsignal();
542 //
543 if (eta < 0.8 && track->GetTPCsignalN() > 70) fDeDx->Fill(ptot, tpcSignal, sign);
544 //
545 // 3.a. calculate expected signals
546 //
547 Float_t sigKaon = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kKaon);
548 Float_t sigProton = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kProton);
549 Float_t sigPion = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kPion);
550 //
551 // 3.b. fill histograms
552 //
553 fHistClassicKaon->Fill(sign*pT,rapKaon,(tpcSignal-sigKaon)/sigKaon);
554 fHistClassicProton->Fill(sign*pT,rapProton,(tpcSignal-sigProton)/sigProton);
555 fHistClassicPion->Fill(sign*pT,rapPion,(tpcSignal-sigPion)/sigPion);
556 //
557 /* 2sigma PID with 2sigma eff correction! */
558 // PION
559 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kPion)) < kNsigmaCut) {
560 fHistPtEtaPion->Fill(trackCounter, rapPion, sign*pT, k2sigmaCorr);
561 //
562 if (trackCounter < 300 && fMCtrue) {
563 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
564 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
565 if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,rapPion,sign*pT,k2sigmaCorr);
566 if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
567 fHistEffPion->Fill(1,rapPion,sign*pT,k2sigmaCorr);
568 if (trackMC->GetFirstMother() > 0) {
569 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
570 Int_t code = TMath::Abs(trackMother->GetPdgCode());
571 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffPion->Fill(3,rapPion,sign*pT,k2sigmaCorr);
572 }
573 }
574 if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,rapPion,sign*pT,k2sigmaCorr);
575 if (TMath::Abs(trackMC->GetPdgCode()) == 13) fHistEffPion->Fill(4,rapPion,sign*pT,k2sigmaCorr);
576 }
577 }
578 // KAON
579 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < kNsigmaCut) {
580 fHistPtEtaKaon->Fill(trackCounter, rapKaon, sign*pT, k2sigmaCorr);
581 //
582 if (trackCounter < 300 && fMCtrue) {
583 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
584 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
585 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,rapKaon,sign*pT,k2sigmaCorr);
586 if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
587 fHistEffKaon->Fill(1,rapKaon,sign*pT,k2sigmaCorr);
588 if (trackMC->GetFirstMother() > 0) {
589 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
590 Int_t code = TMath::Abs(trackMother->GetPdgCode());
591 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffKaon->Fill(3,rapKaon,sign*pT,k2sigmaCorr);
592 }
593 }
594 if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,rapKaon,sign*pT,k2sigmaCorr);
595 }
596 }
597 // KAON NO KINK
598 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < kNsigmaCut && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, rapKaon, sign*pT, k2sigmaCorr);
599 // PROTON
600 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)) < kNsigmaCut) {
601 fHistPtEtaProton->Fill(trackCounter, rapProton, sign*pT, k2sigmaCorr);
602 //
603 track->GetImpactParameters(dca,cov);
604 Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
605 fHistPtEtaProtonDCA->Fill(primVtxDCA, rapProton, sign*pT, k2sigmaCorr);
606 //
607 if (trackCounter < 300 && fMCtrue) {
608 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
609 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
610 if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
611 fHistEffProton->Fill(0.,rapProton,sign*pT,k2sigmaCorr);
612 if (rapProton < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
613 }
614 if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
615 fHistEffProton->Fill(1,rapProton,sign*pT,k2sigmaCorr);
616 if (rapProton < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
617 if (trackMC->GetFirstMother() > 0) {
618 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
619 Int_t code = TMath::Abs(trackMother->GetPdgCode());
620 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
621 fHistEffProton->Fill(3,rapProton,sign*pT,k2sigmaCorr);
622 if (rapProton < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
623 }
624 if (code!=3122&&code!=3222&&code!=3212&&code!=3112&&code!=3322&&code!=3312&&code!=3332&&code!=130&&code!=310) fSecProtons->Fill(trackMC->Vx(),trackMC->Vy());
625 }
626 }
627 if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,rapProton,sign*pT,k2sigmaCorr);
628 }
629 }
630
631 delete track;
632
633 } // end of track loop
634
635 // Post output data
636 PostData(1, fListHist);
637
635c8b48 638}
639
640
641//________________________________________________________________________
642void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
643{
644 // Draw result to the screen
645 // Called once at the end of the query
3c038b30 646 Printf("*** CONSTRUCTOR CALLED ****");
635c8b48 647
648}
649
650
3c038b30 651//________________________________________________________________________
652Bool_t AliAnalysisTaskChargedHadronSpectra::SelectOnImpPar(AliESDtrack* t) {
653 //
654 // cut on transverse impact parameter
655 //
656 Float_t d0z0[2],covd0z0[3];
657 t->GetImpactParameters(d0z0,covd0z0);
658 Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
659 Float_t d0max = 7.*sigma;
660 //
661 Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
662 if (t->Pt() > 1) sigmaZ = 0.0216;
663 Float_t d0maxZ = 5.*sigmaZ;
664 //
665 if (TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
666 return kFALSE;
667}
635c8b48 668
669
670//________________________________________________________________________
671void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
672
673 // Method for the correct logarithmic binning of histograms
674
675 TAxis *axis = h->GetXaxis();
676 int bins = axis->GetNbins();
677
678 Double_t from = axis->GetXmin();
679 Double_t to = axis->GetXmax();
680 Double_t *newBins = new Double_t[bins + 1];
681
682 newBins[0] = from;
683 Double_t factor = pow(to/from, 1./bins);
684
685 for (int i = 1; i <= bins; i++) {
686 newBins[i] = factor * newBins[i-1];
687 }
688 axis->Set(bins, newBins);
4ce766eb 689 delete [] newBins;
635c8b48 690
691}
692
693
694//________________________________________________________________________
695Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
696 //
697 // get the process type of the event.
698 //
699
700 // can only read pythia headers, either directly or from cocktalil header
701 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
702
703 if (!pythiaGenHeader) {
704
705 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
706 if (!genCocktailHeader) {
707 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
708 return -1;
709 }
710
711 TList* headerList = genCocktailHeader->GetHeaders();
712 if (!headerList) {
713 return -1;
714 }
715
716 for (Int_t i=0; i<headerList->GetEntries(); i++) {
717 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
718 if (pythiaGenHeader)
719 break;
720 }
721
722 if (!pythiaGenHeader) {
723 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
724 return -1;
725 }
726 }
727
728 if (adebug) {
729 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
730 }
731
732 return pythiaGenHeader->ProcessType();
733}
734
735
736//________________________________________________________________________
3c038b30 737TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBinLow, Int_t EtaBinHigh,const Double_t * AlephParams) {
635c8b48 738 //
739 // The multiple Gauss fitting is happening here...
740 // input histogram is of the form (Pt, eta, difference in dEdx)
741 //
742
743 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
744 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
745 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
746 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
747 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
748 //
749 foProton->SetParameters(AlephParams);
750 foPion->SetParameters(AlephParams);
751 foElec->SetParameters(AlephParams);
752 foKaon->SetParameters(AlephParams);
753 foMuon->SetParameters(AlephParams);
754
755
756 const Double_t kSigma = 0.06; // expected resolution (roughly)
757
3c038b30 758 const Int_t kPtBins = 2*30/*input->GetXaxis()->GetNbins()*/;
635c8b48 759 const Double_t kPtMax = input->GetXaxis()->GetXmax();
760 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
761
762 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
763 // sort pT-bins ..
3c038b30 764 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0};
635c8b48 765 Int_t indexes[kPtBins+1];
766 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
767 Double_t binsPt[kPtBins+1];
768 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
769 histPt->GetXaxis()->Set(kPtBins, binsPt);
770 //
771
772 TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
773 canvMany->Print("canvManyProton.ps[");
774
775 for(Int_t x=1; x < kPtBins+1; x++) {
776 Double_t pT = input->GetXaxis()->GetBinCenter(x);
3c038b30 777 cout << "Now fitting: " << pT << endl;
635c8b48 778 for(Int_t y=1; y < kEtaBins+1; y++) {
3c038b30 779 Double_t rapidity = input->GetYaxis()->GetBinCenter(y);
780 Double_t eta = TMath::ASinH(TMath::Sqrt(1 + (0.938*0.938)/(pT*pT))*TMath::SinH(rapidity));
635c8b48 781 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
782 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
3c038b30 783 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y + EtaBinHigh - EtaBinLow);
635c8b48 784 histDeDx->SetTitle(Form("pt_%f",pT));
785 Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
786 //
787 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
3c038b30 788 Double_t parameters[12] = {maxYield/2.,0,kSigma,
789 maxYield/2.,(foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),
790 maxYield/2.,(foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),
791 maxYield/2.,(foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot))};
635c8b48 792 gausFit->SetParameters(parameters);
3c038b30 793 if (TMath::Abs(pT) < 0.7) {
794 gausFit = new TF1("gausFit", "gaus(0)",-1/*-7*kSigma*/,1/*7*kSigma*/);
795 gausFit->SetRange(-6*kSigma,6*kSigma);
796 gausFit->SetParameters(parameters);
797 gausFit->SetParLimits(0,0.,maxYield);
798 //for(Int_t ipar = 3; ipar < 12; ipar++) gausFit->FixParameter(ipar,0.);
799 gausFit->SetParLimits(1,-0.15,0.15);
800 gausFit->SetParLimits(2,0.02,0.18);
801 //
802 } else {
803 Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
804 gausFit->SetParLimits(3,0.,maxYield);
805 gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
806 gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
807 //
808 Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
809 gausFit->SetParLimits(6,0.,maxYield);
810 gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
811 gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
812 //
813 Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
814 gausFit->SetParLimits(9,0.,maxYield);
815 gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
816 gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
817 }
635c8b48 818 histDeDx->Fit("gausFit","QNR");
819 gausFit->GetParameters(parameters);
820 // visualization
3c038b30 821 if (y == EtaBinLow) {
635c8b48 822 canvMany->cd(x);
823 gPad->SetLogy();
824 histDeDx->SetMarkerStyle(21);
825 histDeDx->SetMarkerSize(0.7);
826 histDeDx->Draw("E");
827 gausFit->SetLineWidth(2);
828 gausFit->Draw("same");
829 //
830 TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
831 TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
832 TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
833 TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
834 gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
835 gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
836 gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
837 gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
838 //
839 gausFit0.SetLineColor(4);
840 gausFit1.SetLineColor(2);
841 gausFit2.SetLineColor(6);
842 gausFit3.SetLineColor(8);
843 //
844 gausFit0.SetLineWidth(1);
845 gausFit1.SetLineWidth(1);
846 gausFit2.SetLineWidth(1);
847 gausFit3.SetLineWidth(1);
848 //
849 gausFit0.Draw("same");
850 gausFit1.Draw("same");
851 gausFit2.Draw("same");
852 gausFit3.Draw("same");
3c038b30 853 if (TMath::Abs(pT) < 2) canvMany->Print("canvManyProton.ps");
635c8b48 854 }
855
856 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
857 delete gausFit;
858 //
3c038b30 859 if (y == EtaBinLow && yield > 0) histPt->Fill(pT,yield);
635c8b48 860 //delete gausYield;
861 }
862 }
863 canvMany->Print("canvManyProton.ps]");
864
865 TCanvas * canvPt = new TCanvas();
866 canvPt->cd();
867 histPt->Draw();
868
869 return histPt;
870}
871
872
873
874//________________________________________________________________________
3c038b30 875TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBinLow, Int_t EtaBinHigh,const Double_t * AlephParams) {
635c8b48 876 //
877 // The multiple Gauss fitting is happening here...
3c038b30 878 // input histogram is of the form (Pt, eta, difference in dEdx)
635c8b48 879 //
3c038b30 880
635c8b48 881 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
882 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
883 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
884 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
885 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
886 //
887 foProton->SetParameters(AlephParams);
888 foPion->SetParameters(AlephParams);
889 foElec->SetParameters(AlephParams);
890 foKaon->SetParameters(AlephParams);
891 foMuon->SetParameters(AlephParams);
892
3c038b30 893
635c8b48 894 const Double_t kSigma = 0.06; // expected resolution (roughly)
895
3c038b30 896 const Int_t kPtBins = 2*30/*input->GetXaxis()->GetNbins()*/;
635c8b48 897 const Double_t kPtMax = input->GetXaxis()->GetXmax();
898 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
3c038b30 899
635c8b48 900 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
901 // sort pT-bins ..
3c038b30 902 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0};
635c8b48 903 Int_t indexes[kPtBins+1];
904 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
905 Double_t binsPt[kPtBins+1];
906 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
907 histPt->GetXaxis()->Set(kPtBins, binsPt);
908 //
909
3c038b30 910 TCanvas * canvMany = new TCanvas("canvManyPion","canvManyPion");
911 canvMany->Print("canvManyPion.ps[");
635c8b48 912
913 for(Int_t x=1; x < kPtBins+1; x++) {
914 Double_t pT = input->GetXaxis()->GetBinCenter(x);
915 for(Int_t y=1; y < kEtaBins+1; y++) {
3c038b30 916 Double_t rapidity = input->GetYaxis()->GetBinCenter(y);
917 Double_t eta = TMath::ASinH(TMath::Sqrt(1 + (0.1395*0.1395)/(pT*pT))*TMath::SinH(rapidity));
635c8b48 918 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
919 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
3c038b30 920 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y + EtaBinHigh - EtaBinLow);
921 histDeDx->SetTitle(Form("pt_%f",pT));
922 Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
635c8b48 923 //
3c038b30 924 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
925 Double_t parameters[12] = {maxYield/2.,0,kSigma,
926 maxYield/2.,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foPion->Eval(ptot)/foPion->Eval(ptot)),
927 maxYield/2.,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)),
928 maxYield/2.,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot))};
635c8b48 929 gausFit->SetParameters(parameters);
3c038b30 930 gausFit->SetParLimits(0,0.,maxYield);
931 gausFit->SetParLimits(1,-0.05,0.05);
932 gausFit->SetParLimits(2,0.04,0.08);
933 //
934 Double_t protonPosition = (foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
935 gausFit->SetParLimits(3,0.,maxYield);
936 gausFit->SetParLimits(4, 0.95*protonPosition, 1.05*protonPosition);
937 gausFit->SetParLimits(5,0.8*kSigma*(foProton->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foProton->Eval(ptot)/foPion->Eval(ptot)));
938 //
939 Double_t elecPosition = (foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
940 gausFit->SetParLimits(6,0.,maxYield);
941 gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
942 gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)));
635c8b48 943 //
3c038b30 944 Double_t kaonPosition = (foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
945 gausFit->SetParLimits(9,0.,maxYield);
946 gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
947 gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot)));
948 histDeDx->Fit("gausFit","QNR");
949 gausFit->GetParameters(parameters);
950 // visualization
951 if (y == EtaBinLow) {
952 canvMany->cd(x);
953 gPad->SetLogy();
954 histDeDx->SetMarkerStyle(21);
955 histDeDx->SetMarkerSize(0.7);
956 histDeDx->Draw("E");
957 gausFit->SetLineWidth(2);
635c8b48 958 gausFit->Draw("same");
3c038b30 959 //
960 TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
961 TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
962 TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
963 TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
964 gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
965 gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
966 gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
967 gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
968 //
969 gausFit0.SetLineColor(4);
970 gausFit1.SetLineColor(2);
971 gausFit2.SetLineColor(6);
972 gausFit3.SetLineColor(8);
973 //
974 gausFit0.SetLineWidth(1);
975 gausFit1.SetLineWidth(1);
976 gausFit2.SetLineWidth(1);
977 gausFit3.SetLineWidth(1);
978 //
979 gausFit0.Draw("same");
980 gausFit1.Draw("same");
981 gausFit2.Draw("same");
982 gausFit3.Draw("same");
983 if (TMath::Abs(pT) < 2) canvMany->Print("canvManyPion.ps");
635c8b48 984 }
3c038b30 985
986 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
987 delete gausFit;
635c8b48 988 //
3c038b30 989 if (y == EtaBinLow && yield > 0) histPt->Fill(pT,yield);
990 //delete gausYield;
635c8b48 991 }
992 }
3c038b30 993 canvMany->Print("canvManyPion.ps]");
994
635c8b48 995 TCanvas * canvPt = new TCanvas();
996 canvPt->cd();
997 histPt->Draw();
635c8b48 998
3c038b30 999 return histPt;
635c8b48 1000}
1001
1002
1003//________________________________________________________________________
3c038b30 1004TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBinLow, Int_t EtaBinHigh,const Double_t * AlephParams) {
635c8b48 1005 //
1006 // The multiple Gauss fitting is happening here...
3c038b30 1007 // input histogram is of the form (Pt, eta, difference in dEdx)
635c8b48 1008 //
1009
1010 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
1011 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
1012 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
1013 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
1014 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
1015 //
1016 foProton->SetParameters(AlephParams);
1017 foPion->SetParameters(AlephParams);
1018 foElec->SetParameters(AlephParams);
1019 foKaon->SetParameters(AlephParams);
1020 foMuon->SetParameters(AlephParams);
1021
3c038b30 1022
1023 const Double_t kSigma = 0.055; // expected resolution (roughly)
635c8b48 1024
3c038b30 1025 const Int_t kPtBins = 2*30/*input->GetXaxis()->GetNbins()*/;
635c8b48 1026 const Double_t kPtMax = input->GetXaxis()->GetXmax();
3c038b30 1027 //const Int_t kEtaBins = input->GetYaxis()->GetNbins();
635c8b48 1028
1029 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
1030 // sort pT-bins ..
3c038b30 1031 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0};
635c8b48 1032 Int_t indexes[kPtBins+1];
1033 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
1034 Double_t binsPt[kPtBins+1];
1035 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
1036 histPt->GetXaxis()->Set(kPtBins, binsPt);
1037 //
1038
3c038b30 1039 TCanvas * canvMany = new TCanvas("canvManyKaon","canvManyKaon");
1040 canvMany->Print("canvManyKaon.ps[");
635c8b48 1041
1042 for(Int_t x=1; x < kPtBins+1; x++) {
1043 Double_t pT = input->GetXaxis()->GetBinCenter(x);
3c038b30 1044 for(Int_t y=1; y < 2/*kEtaBins+1*/; y++) {
1045 Double_t rapidity = input->GetYaxis()->GetBinCenter(y);
1046 Double_t eta = TMath::ASinH(TMath::Sqrt(1 + (0.4936*0.4936)/(pT*pT))*TMath::SinH(rapidity));
635c8b48 1047 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1048 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
3c038b30 1049 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y + EtaBinHigh - EtaBinLow);
1050 histDeDx->SetTitle(Form("pt_%f",pT));
1051 Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
635c8b48 1052 //
3c038b30 1053 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1 /*-5*kSigma*/, 1/*5*kSigma*/);
1054 Double_t parameters[12] = {maxYield*0.1,0,kSigma,
1055 maxYield*0.8,(foPion->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)),
1056 maxYield*0.01,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)),
1057 maxYield*0.1,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot))};
635c8b48 1058 gausFit->SetParameters(parameters);
3c038b30 1059 gausFit->SetParLimits(0,0.,maxYield);
1060 gausFit->SetParLimits(1,-0.15,0.15);
1061 gausFit->SetParLimits(2,0.04,0.08);
1062 //
1063 Double_t pionPosition = (foPion->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
1064 //Double_t elecPosition = (foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
1065 Double_t protonPosition = (foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
635c8b48 1066 //
3c038b30 1067 if (TMath::Abs(pT) < 0.4) {
1068 cout << " fixing the parameters " << endl;
1069 gausFit->SetParameter(3,0.); gausFit->SetParameter(6,0.); gausFit->SetParameter(9,0.);
1070 gausFit->FixParameter(3,0); gausFit->FixParameter(6,0); gausFit->FixParameter(9,0);
1071 //gausFit->SetParLimits(0,0.,maxYield);
1072 //gausFit->SetParLimits(1,-0.15,0.15);
1073 //gausFit->SetParLimits(2,0.03,0.18);
1074 } else {
1075 //
1076 gausFit->SetParLimits(0,0.,0.5*maxYield);
1077 gausFit->SetParLimits(1,-0.05,0.05);
1078 gausFit->SetParLimits(2,0.04,0.08);
1079 //
1080 gausFit->SetParLimits(3,0.,maxYield);
1081 gausFit->SetParLimits(4, 0.5*pionPosition, 1.5*pionPosition);
1082 gausFit->SetParLimits(5,0.6*kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)));
1083 //
1084 gausFit->FixParameter(6,0);
1085 gausFit->FixParameter(7,0);
1086 gausFit->FixParameter(8,0);
1087 /*
1088 gausFit->SetParLimits(6,0.,0.2*maxYield);
1089 gausFit->SetParLimits(7, 0.9*elecPosition, 1.1*elecPosition);
1090 gausFit->SetParLimits(8,0.6*kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)));
1091 */
1092 //
1093 gausFit->SetParLimits(9,0.,0.2*maxYield);
1094 gausFit->SetParLimits(10, 0.9*protonPosition, 1.1*protonPosition);
1095 gausFit->SetParLimits(11,0.6*kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot)));
1096 }
1097 histDeDx->Fit("gausFit","QNR");
1098 gausFit->GetParameters(parameters);
1099 // visualization
1100 if (y == EtaBinLow) {
635c8b48 1101 canvMany->cd(x);
3c038b30 1102 //gPad->SetLogy();
1103 histDeDx->SetMarkerStyle(21);
1104 histDeDx->SetMarkerSize(0.7);
1105 histDeDx->Draw("E");
1106 gausFit->SetLineWidth(2);
635c8b48 1107 gausFit->Draw("same");
3c038b30 1108 //
1109 TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
1110 TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
1111 TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
1112 TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
1113 gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
1114 gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
1115 gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
1116 gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
1117 //
1118 gausFit0.SetLineColor(4);
1119 gausFit1.SetLineColor(2);
1120 gausFit2.SetLineColor(6);
1121 gausFit3.SetLineColor(8);
1122 //
1123 gausFit0.SetLineWidth(1);
1124 gausFit1.SetLineWidth(1);
1125 gausFit2.SetLineWidth(1);
1126 gausFit3.SetLineWidth(1);
1127 //
1128 gausFit0.Draw("same");
1129 gausFit1.Draw("same");
1130 gausFit2.Draw("same");
1131 gausFit3.Draw("same");
1132 if (TMath::Abs(pT) < 2) canvMany->Print("canvManyKaon.ps");
1133 //
1134 // propaganda slots for the paper
1135 //
1136 if (pT > 0.374 && pT < 0.376) {
1137 TFile f1("slice1.root","RECREATE");
1138 canvMany->Write();
1139 f1.Close();
1140 }
1141 if (pT > 0.524 && pT < 0.526) {
1142 TFile f1("slice2.root","RECREATE");
1143 canvMany->Write();
1144 f1.Close();
1145 }
1146 if (pT > 0.674 && pT < 0.676) {
1147 TFile f1("slice3.root","RECREATE");
1148 canvMany->Write();
1149 f1.Close();
1150 }
1151 if (pT > 0.624 && pT < 0.626) {
1152 TFile f1("slice4.root","RECREATE");
1153 canvMany->Write();
1154 f1.Close();
1155 }
635c8b48 1156 }
3c038b30 1157
1158 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
1159 cout << pT << " res.: " << parameters[2] << " " << " yield: " << yield << endl;
1160 delete gausFit;
635c8b48 1161 //
3c038b30 1162 if (y == EtaBinLow && yield > 0) histPt->Fill(pT,yield);
1163 //delete gausYield;
635c8b48 1164 }
1165 }
3c038b30 1166 canvMany->Print("canvManyKaon.ps]");
635c8b48 1167
1168 TCanvas * canvPt = new TCanvas();
1169 canvPt->cd();
1170 histPt->Draw();
1171
1172 return histPt;
635c8b48 1173}
1174
1175
1176//________________________________________________________________________
1177void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const TList * ListOfHistogramsData, const Char_t *filename) {
1178
1179 // Define ranges
1180 const Int_t kMultLow = 1;
1181 const Int_t kMultHigh = 10;
1182
1183 const Int_t kEtaHigh = 4;
1184
1185 // Extract histograms for the MC efficiency study
1186 TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1187 TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1188 TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
635c8b48 1189
1190 TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
635c8b48 1191 TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1192 TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
635c8b48 1193
1194 TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1195 TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1196 TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1197
1198 // Extract data for the final spectra
1199 TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
635c8b48 1200 TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1201 TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
635c8b48 1202
1203 // "MC" of the real data for consistency checks
1204 TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1205 TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1206 TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
635c8b48 1207
3c038b30 1208
635c8b48 1209 TFile spectraFile(filename,"RECREATE");
1210
1211 // 1. Protons
1212 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1213 TH1D *protonRaw = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1214 protonRaw->Sumw2();
1215 TH1D *protonMC = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1216 protonMC->Sumw2();
1217 TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1218 protonRawPrim->Sumw2();
1219 TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1220 protonSecReac->Sumw2();
1221 TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1222 protonSecWeak->Sumw2();
1223 TH1D *protonMis = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1224 protonMis->Sumw2();
1225 protonSecReac->Add(protonSecWeak,-1);
1226 //
1227 // 1. total (physical) efficiency
1228 //
1229 TH1D *protonEffTot = new TH1D(*protonRaw);
1230 protonEffTot->Divide(protonRaw,protonMC);
1231 protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1232 protonEffTot->Write();
1233 //
1234 // 2. contributions from secondaries created by interactions
1235 //
1236 TH1D *protonEffSecReac = new TH1D(*protonRaw);
1237 TH1D *dummy = new TH1D(*protonRaw);
1238 dummy->Add(protonSecReac,-1);
1239 protonEffSecReac->Divide(protonRaw,dummy);
1240 protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1241 protonEffSecReac->Write();
1242 delete dummy;
1243 //
1244 // 3. contributions from secondaries from weak decays
1245 //
1246 TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1247 TH1D *dummy2 = new TH1D(*protonRaw);
1248 dummy2->Add(protonSecWeak,-1);
1249 protonEffSecWeak->Divide(protonRaw,dummy2);
1250 protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1251 protonEffSecWeak->Write();
1252 delete dummy2;
1253 //
1254 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1255 //
1256 TH1D *protonEffAbs = new TH1D(*protonRaw);
1257 protonEffAbs->Divide(protonRawPrim,protonMC);
1258 protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1259 protonEffAbs->Write();
1260 //
1261 // 5. contamination
1262 //
1263 TH1D *protonEffMis = new TH1D(*protonMis);
1264 protonEffMis->Divide(protonMis,protonRaw);
1265 protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1266 protonEffMis->Write();
1267
1268 //
1269 // PROCESS THE REAL DATA FROM HERE ON
1270 //
1271 TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1272 histPtProtonEta->Write();
1273 //
1274 // Apply corrections as you would like them to have...
1275 //
1276 TH1D *protonFinal = new TH1D(*histPtProtonEta);
1277 protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1278 protonFinal->Sumw2();
1279 protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1280 protonFinal->Write();
1281 //
1282 // if available save also the MC truth for consistency checks
1283 //
1284 TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1285 histPtProtonEtaMcData->Write();
1286 }
1287
1288 // 2. Kaons
1289 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1290 TH1D *kaonRaw = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1291 kaonRaw->Sumw2();
1292 TH1D *kaonMC = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1293 kaonMC->Sumw2();
1294 TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1295 kaonRawPrim->Sumw2();
1296 TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1297 kaonSecReac->Sumw2();
1298 TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1299 kaonSecWeak->Sumw2();
1300 TH1D *kaonMis = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1301 kaonMis->Sumw2();
1302 kaonSecReac->Add(kaonSecWeak,-1);
1303 //
1304 // 1. total (physical) efficiency
1305 //
1306 TH1D *kaonEffTot = new TH1D(*kaonRaw);
1307 kaonEffTot->Divide(kaonRaw,kaonMC);
1308 kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1309 kaonEffTot->Write();
1310 //
1311 // 2. contributions from secondaries created by interactions
1312 //
1313 TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1314 TH1D *dummy = new TH1D(*kaonRaw);
1315 dummy->Add(kaonSecReac,-1);
1316 kaonEffSecReac->Divide(kaonRaw,dummy);
1317 kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1318 kaonEffSecReac->Write();
1319 delete dummy;
1320 //
1321 // 3. contributions from secondaries from weak decays
1322 //
1323 TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1324 TH1D *dummy2 = new TH1D(*kaonRaw);
1325 dummy2->Add(kaonSecWeak,-1);
1326 kaonEffSecWeak->Divide(kaonRaw,dummy2);
1327 kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1328 kaonEffSecWeak->Write();
1329 delete dummy2;
1330 //
1331 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1332 //
1333 TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1334 kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1335 kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1336 kaonEffAbs->Write();
1337 //
1338 // 5. contamination
1339 //
1340 TH1D *kaonEffMis = new TH1D(*kaonMis);
1341 kaonEffMis->Divide(kaonMis,kaonRaw);
1342 kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1343 kaonEffMis->Write();
1344
1345 //
1346 // PROCESS THE REAL DATA FROM HERE ON
1347 //
1348 TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1349 histPtKaonEta->Write();
1350 //
1351 // Apply corrections as you would like them to have...
1352 //
1353 TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1354 kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1355 kaonFinal->Sumw2();
1356 kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1357 kaonFinal->Write();
1358 //
1359 // if available save also the MC truth for consistency checks
1360 //
1361 TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1362 histPtKaonEtaMcData->Write();
1363 }
1364
1365
1366 // 3. Pions
1367 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1368 TH1D *pionRaw = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1369 pionRaw->Sumw2();
1370 TH1D *pionMC = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1371 pionMC->Sumw2();
1372 TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1373 pionRawPrim->Sumw2();
1374 TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1375 pionSecReac->Sumw2();
1376 TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1377 pionSecWeak->Sumw2();
1378 TH1D *pionMis = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1379 pionMis->Sumw2();
1380 pionSecReac->Add(pionSecWeak,-1);
1381 //
1382 // 1. total (physical) efficiency
1383 //
1384 TH1D *pionEffTot = new TH1D(*pionRaw);
1385 pionEffTot->Divide(pionRaw,pionMC);
1386 pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1387 pionEffTot->Write();
1388 //
1389 // 2. contributions from secondaries created by interactions
1390 //
1391 TH1D *pionEffSecReac = new TH1D(*pionRaw);
1392 TH1D *dummy = new TH1D(*pionRaw);
1393 dummy->Add(pionSecReac,-1);
1394 pionEffSecReac->Divide(pionRaw,dummy);
1395 pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1396 pionEffSecReac->Write();
1397 delete dummy;
1398 //
1399 // 3. contributions from secondaries from weak decays
1400 //
1401 TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1402 TH1D *dummy2 = new TH1D(*pionRaw);
1403 dummy2->Add(pionSecWeak,-1);
1404 pionEffSecWeak->Divide(pionRaw,dummy2);
1405 pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1406 pionEffSecWeak->Write();
1407 delete dummy2;
1408 //
1409 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1410 //
1411 TH1D *pionEffAbs = new TH1D(*pionRaw);
1412 pionEffAbs->Divide(pionRawPrim,pionMC);
1413 pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1414 pionEffAbs->Write();
1415 //
1416 // 5. contamination
1417 //
1418 TH1D *pionEffMis = new TH1D(*pionMis);
1419 pionEffMis->Divide(pionMis,pionRaw);
1420 pionEffMis->SetName(Form("PionContam_%i",iEta));
1421 pionEffMis->Write();
1422
1423 //
1424 // PROCESS THE REAL DATA FROM HERE ON
1425 //
1426 TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1427 histPtPionEta->Write();
1428 //
1429 // Apply corrections as you would like them to have...
1430 //
1431 TH1D *pionFinal = new TH1D(*histPtPionEta);
1432 pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1433 pionFinal->Sumw2();
1434 pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1435 pionFinal->Write();
1436 //
1437 // if available save also the MC truth for consistency checks
1438 //
1439 TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1440 histPtPionEtaMcData->Write();
1441
1442 }
1443
1444 //
1445 // Close the file after everthing is done
1446 //
1447 spectraFile.Close();
1448
1449}
1450
1451
1452
1453