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