Bugfix: memory leak in the histogram class
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisTaskChargedHadronSpectra.cxx
CommitLineData
635c8b48 1
2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17///////////////////////////////////////////////////////////////////////////
18// //
19// //
20// Analysis for identified particle spectra measured with TPC dE/dx. //
21// //
22// //
23///////////////////////////////////////////////////////////////////////////
24
25#include "Riostream.h"
26#include "TChain.h"
27#include "TTree.h"
28#include "TH1F.h"
29#include "TH2F.h"
30#include "TH3F.h"
31#include "TList.h"
32#include "TMath.h"
33#include "TCanvas.h"
34#include "TObjArray.h"
35#include "TF1.h"
36#include "TFile.h"
37
38#include "AliAnalysisTaskSE.h"
39#include "AliAnalysisManager.h"
40
41#include "AliHeader.h"
42#include "AliGenPythiaEventHeader.h"
43#include "AliGenCocktailEventHeader.h"
44
45#include "AliPID.h"
46#include "AliESDtrackCuts.h"
47#include "AliESDVertex.h"
48#include "AliESDEvent.h"
49#include "AliESDInputHandler.h"
50#include "AliESDtrack.h"
10d100d4 51#include "AliESDpid.h"
635c8b48 52
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliStack.h"
56
57#include "AliLog.h"
58
59#include "AliAnalysisTaskChargedHadronSpectra.h"
60
61
62ClassImp(AliAnalysisTaskChargedHadronSpectra)
63
64//________________________________________________________________________
65AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra()
10d100d4 66 : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
635c8b48 67 fMCtrue(0),
68 fAlephParameters(),
69 fHistPtMCKaon(0),
70 fHistPtMCProton(0),
71 fHistPtMCPion(0),
72 fHistPtMCElectron(0),
73 fHistPtMCMuon(0),
74 fHistPtEtaKaon(0),
75 fHistPtEtaKaonNoKink(0),
76 fHistPtEtaProton(0),
77 fHistPtEtaProtonDCA(0),
78 fHistPtEtaPion(0),
79 fHistPtEtaElectron(0),
80 fHistClassicKaon(0),
81 fHistClassicProton(0),
82 fHistClassicPion(0),
83 fHistClassicElectron(0),
84 fDeDx(0),
85 fHistTrackPerEvent(0),
86 fHistTrackPerEventMC(0),
87 fSecProtons(0),
88 fVertexZ(0),
89 fHistEtaNcls(0),
90 fHistEtaPhi(0),
91 fHistEffProton(0),
92 fHistEffProtonDCA(0),
93 fHistEffPion(0),
94 fHistEffKaon(0),
95 fHighPtElectrons(0),
96 fHighPtHadrons(0)
97{
98 // default Constructor
99}
100
101
102//________________________________________________________________________
103AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra(const char *name)
10d100d4 104 : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
635c8b48 105 fMCtrue(0),
106 fAlephParameters(),
107 fHistPtMCKaon(0),
108 fHistPtMCProton(0),
109 fHistPtMCPion(0),
110 fHistPtMCElectron(0),
111 fHistPtMCMuon(0),
112 fHistPtEtaKaon(0),
113 fHistPtEtaKaonNoKink(0),
114 fHistPtEtaProton(0),
115 fHistPtEtaProtonDCA(0),
116 fHistPtEtaPion(0),
117 fHistPtEtaElectron(0),
118 fHistClassicKaon(0),
119 fHistClassicProton(0),
120 fHistClassicPion(0),
121 fHistClassicElectron(0),
122 fDeDx(0),
123 fHistTrackPerEvent(0),
124 fHistTrackPerEventMC(0),
125 fSecProtons(0),
126 fVertexZ(0),
127 fHistEtaNcls(0),
128 fHistEtaPhi(0),
129 fHistEffProton(0),
130 fHistEffProtonDCA(0),
131 fHistEffPion(0),
132 fHistEffKaon(0),
133 fHighPtElectrons(0),
134 fHighPtHadrons(0)
135{
136 //
137 // standard constructur which should be used - PID objects is initialized
138 //
139
140 fMCtrue = kTRUE; // change two things for processing real data: 1. set this to false! / 2. change ALEPH parameters
141 fAlephParameters[0] = 4.23232575531564326e+00;//50*0.76176e-1;
142 fAlephParameters[1] = 8.68482806165147636e+00;//10.632;
143 fAlephParameters[2] = 1.34000000000000005e-05;//0.13279e-4;
144 fAlephParameters[3] = 2.30445734159456084e+00;//1.8631;
145 fAlephParameters[4] = 2.25624744086878559e+00;//1.9479;
146
10d100d4 147 fESDpid = new AliESDpid();
148 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
635c8b48 149
150 // Constructor
151 Printf("*** CONSTRUCTOR CALLED ****");
152 // Define input and output slots here
153 // Input slot #0 works with a TChain
c77fc043 154 //DefineInput(0, TChain::Class()); <-> not needed in AliAnalysisTaskSE
635c8b48 155
156 // Output slot #0 writes into a TList container
c77fc043 157 DefineOutput(1, TList::Class());
635c8b48 158
159
160
161}
162
635c8b48 163
164//________________________________________________________________________
c77fc043 165void AliAnalysisTaskChargedHadronSpectra::UserCreateOutputObjects()
635c8b48 166{
167 // Create histograms
168 // Called once
169 fListHist = new TList();
c77fc043 170 //fListHist->SetOwner(); // Whoever knows how the data handling is ...?
635c8b48 171
172 const Int_t kPtBins = 2*56;
173 const Double_t kPtMax = 16.0;
174 const Int_t kEtaBins = 4;
175 const Double_t kEtaMax = 0.8;
176 const Int_t kDeDxBins = 200;
177 const Double_t kDeDxMax = 1;
178 const Int_t kMultBins = 10;
179 const Int_t kMultMax = 300;
180
181 // sort pT-bins ..
182 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, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.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, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
183 Int_t indexes[kPtBins+1];
184 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
185 Double_t binsPt[kPtBins+1];
186 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
187
188
189 // MC histograms
190 fHistPtMCKaon = new TH3F("HistPtMCKaon", "PtEtaKaon; mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
191 fHistPtMCProton = new TH3F("HistPtMCProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
192 fHistPtMCPion = new TH3F("HistPtMCPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
193 fHistPtMCElectron = new TH3F("HistPtMCElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
194 fHistPtMCMuon = new TH3F("HistPtMCMuon", "PtEtaMuon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
195 //
196 fHistPtMCKaon->GetZaxis()->Set(kPtBins, binsPt);
197 fHistPtMCProton->GetZaxis()->Set(kPtBins, binsPt);
198 fHistPtMCPion->GetZaxis()->Set(kPtBins, binsPt);
199 fHistPtMCElectron->GetZaxis()->Set(kPtBins, binsPt);
200 fHistPtMCMuon->GetZaxis()->Set(kPtBins, binsPt);
201 //
202 fListHist->Add(fHistPtMCKaon);
203 fListHist->Add(fHistPtMCProton);
204 fListHist->Add(fHistPtMCPion);
205 fListHist->Add(fHistPtMCElectron);
206 fListHist->Add(fHistPtMCMuon);
207
208 // reconstructed particle histograms
209 fHistPtEtaKaon = new TH3F("HistPtEtaKaon", "PtEtaKaon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
210 fHistPtEtaKaonNoKink = new TH3F("HistPtEtaKaonNoKink", "PtEtaKaonNoKink;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
211 fHistPtEtaProton = new TH3F("HistPtEtaProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
212 fHistPtEtaProtonDCA = new TH3F("HistPtEtaProtonDCA", "PtEtaProton;DCA (cm); #eta; p_{T} (GeV)",200,0,15,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
213 fHistPtEtaPion = new TH3F("HistPtEtaPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
214 fHistPtEtaElectron = new TH3F("HistPtEtaElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
215 //
216 fHistPtEtaKaon->GetZaxis()->Set(kPtBins, binsPt);
217 fHistPtEtaKaonNoKink->GetZaxis()->Set(kPtBins, binsPt);
218 fHistPtEtaProton->GetZaxis()->Set(kPtBins, binsPt);
219 fHistPtEtaProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
220 fHistPtEtaPion->GetZaxis()->Set(kPtBins, binsPt);
221 fHistPtEtaElectron->GetZaxis()->Set(kPtBins, binsPt);
222 //
223 fListHist->Add(fHistPtEtaKaon);
224 fListHist->Add(fHistPtEtaKaonNoKink);
225 fListHist->Add(fHistPtEtaProton);
226 fListHist->Add(fHistPtEtaProtonDCA);
227 fListHist->Add(fHistPtEtaPion);
228 fListHist->Add(fHistPtEtaElectron);
229
230 // histograms for the classical analysis
231 fHistClassicKaon = new TH3F("HistClassicKaon", "PtEtaKaon;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
232 fHistClassicProton = new TH3F("HistClassicProton", "PtEtaProton;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
233 fHistClassicPion = new TH3F("HistClassicPion", "PtEtaPion;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
234 fHistClassicElectron = new TH3F("HistClassicElectron", "PtEtaElectron;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
235 //
236 fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
237 fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
238 fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
239 fHistClassicElectron->GetXaxis()->Set(kPtBins, binsPt);
240 //
241 fListHist->Add(fHistClassicKaon);
242 fListHist->Add(fHistClassicProton);
243 fListHist->Add(fHistClassicPion);
244 fListHist->Add(fHistClassicElectron);
245
246 // histograms of general interest
247 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,20.,1000,0.,500);
248 BinLogX(fDeDx);
249 fListHist->Add(fDeDx);
250 fHistTrackPerEvent = new TH1F("HistTrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",101,-0.5,100);
251 fListHist->Add(fHistTrackPerEvent);
252 fHistTrackPerEventMC = new TH2F("HistTrackPerEventMC", "Tracks per event;Number of Tracks;Number of Events",10,-0.5,9,101,-0.5,100);
253 fListHist->Add(fHistTrackPerEventMC);
254 fSecProtons = new TH2F("SecProtons", "xy;x;y",100,-10,10,100,-10,10);
255 fListHist->Add(fSecProtons);
256 fVertexZ = new TH1F("VertexZ", "vertex position;z (cm);counts",200,-50,50);
257 fListHist->Add(fVertexZ);
258 //
259 fHistEtaNcls = new TH2F("HistEtaNcls", "EtaNcls;#eta;Ncls",100,-1.5,1.5,80,0,160);
260 fListHist->Add(fHistEtaNcls);
261 fHistEtaPhi = new TH2F("HistEtaPhi", "EtaNcls;#eta;#phi",100,-1.5,1.5,100,0,2*TMath::Pi());
262 fListHist->Add(fHistEtaPhi);
263
264 // histograms for a refined efficiency study
265 fHistEffProton = new TH3F("HistEffProton", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
266 fListHist->Add(fHistEffProton);
267 fHistEffProtonDCA = new TH3F("HistEffProtonDCA", "Eff;p_{T} (GeV); code",10,-0.5,9.5,200,0,15,kPtBins,-kPtMax,kPtMax);
268 fListHist->Add(fHistEffProtonDCA);
269 fHistEffPion = new TH3F("HistEffPion", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
270 fListHist->Add(fHistEffPion);
271 fHistEffKaon = new TH3F("HistEffKaon", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
272 fListHist->Add(fHistEffKaon);
273 //
274 fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
275 fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
276 fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
277 fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
278 //
279
280
281 // histograms for dN/dpT
282 fHighPtElectrons = new TH1F("HistHighPtElectrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
283 fListHist->Add(fHighPtElectrons);
284 fHighPtHadrons = new TH1F("HistHighPtHadrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
285 fListHist->Add(fHighPtHadrons);
286
287 //
288 // Create Objects which are needed for the rest of the analysis
289 //
290 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
291 fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
292 fESDtrackCuts->SetMaxNsigmaToVertex(3);
293 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
c77fc043 294 fESDtrackCuts->SetRequireTPCRefit(kFALSE);
635c8b48 295 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
296 fESDtrackCuts->SetMinNClustersTPC(100);
297 fESDtrackCuts->SetMaxChi2PerClusterTPC(3.5);
298 //fESDtrackCuts->SetMaxDCAToVertex(3);
299 //
300 fListHist->Add(fESDtrackCuts);
301 //
302 //
303
304}
305
306//________________________________________________________________________
c77fc043 307void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *)
635c8b48 308{
309 // main event loop
310
311 AliLog::SetGlobalLogLevel(AliLog::kError);
312 //
313 // Check Monte Carlo information and other access first:
314 //
315 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
316 if (!eventHandler) {
317 //Printf("ERROR: Could not retrieve MC event handler");
318 if (fMCtrue) return;
319 }
320
321 AliMCEvent* mcEvent = eventHandler->MCEvent();
322 if (!mcEvent) {
323 //Printf("ERROR: Could not retrieve MC event");
324 if (fMCtrue) return;
325 }
326
c77fc043 327 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
635c8b48 328 if (!fESD) {
329 //Printf("ERROR: fESD not available");
330 return;
331 }
332
333 if (!fESDtrackCuts) {
334 Printf("ERROR: fESDtrackCuts not available");
335 return;
336 }
337
338 Int_t trackCounter = 0;
339
340 // monitor vertex position
341 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
342 if (!vertex) {
343 fHistTrackPerEvent->Fill(trackCounter);
c77fc043 344 PostData(1, fListHist);
635c8b48 345 return;
346 } else {
347 if (vertex->GetZv() != 0) fVertexZ->Fill(vertex->GetZv());
348 if (TMath::Abs(vertex->GetZv()) > 10) {
349 fHistTrackPerEvent->Fill(trackCounter);
c77fc043 350 PostData(1, fListHist);
635c8b48 351 return;
352 }
353 }
354
355
356
357 // 1st track loop to determine multiplicities
358 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
359 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
360 if (!track) continue;
361 if (fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta())< 0.9) {
362 trackCounter++;
363 //
364 }
365 delete track;
366 }
367
368
369 // MC loop
370 AliStack* stack = mcEvent->Stack();
371
372 if (fMCtrue) {
373 //
374 //
375 //
376 AliHeader * header = mcEvent->Header();
377 Int_t processtype = GetPythiaEventProcessType(header);
378 // non diffractive
379 if (processtype !=92 && processtype !=93 && processtype != 94) fHistTrackPerEventMC->Fill(1., trackCounter);
380 // single diffractive
381 if ((processtype == 92 || processtype == 93)) fHistTrackPerEventMC->Fill(2., trackCounter);
382 // double diffractive
383 if (processtype == 94) fHistTrackPerEventMC->Fill(3., trackCounter);
384 //
385 Int_t mult = trackCounter;//stack->GetNprimary();
386 fHistTrackPerEventMC->Fill(0., trackCounter);
387 //
388 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
389 TParticle * trackMC = stack->Particle(i);
390 Int_t pdg = trackMC->GetPdgCode();
391 //
392 Double_t xv = trackMC->Vx();
393 Double_t yv = trackMC->Vy();
394 Double_t zv = trackMC->Vz();
395 Double_t dxy = 0;
396 dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
397 Double_t dz = 0;
398 dz = TMath::Abs(zv); // so stupid to avoid warnings
399 //
400 // vertex cut - selection of primaries
401 //
402 //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
403 if (pdg == 11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select e- only; before vertex
404 if (pdg == -11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select e+ only; before vertex
405 //
406 if (!stack->IsPhysicalPrimary(i)) continue;
407 //
408 if (pdg == 321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select K+ only
409 if (pdg == 211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select Pi+ only
410 if (pdg == 2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select p+ only
411 if (pdg == -13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select mu+ only
412 //
413 if (pdg == -321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select K- only
414 if (pdg == -211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select Pi- only
415 if (pdg == -2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select p- only
416 if (pdg == 13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select mu- only
417 //
418 // special Kaon checks - those which decay before entering the TPC fiducial volume
419 //
420 if (TMath::Abs(pdg)==321) {
421 TParticle * trackDaughterMC = stack->Particle(TMath::Abs(trackMC->GetFirstDaughter()));
422 Int_t pdgDaughter = trackDaughterMC->GetPdgCode();
423 if (pdgDaughter == TMath::Abs(13) && pdg == 321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // Kaon control hist
424 if (pdgDaughter == TMath::Abs(13) && pdg == -321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // Kaon control hist
425 }
426 }
427 }
428
429 //end MC tracks loop
430
431
432 // ESD loop
433 // definition of PID functions --> to be put to CreateOutputObjects()
434
435 TF1 foProton("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
436 TF1 foPion("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
437 TF1 foElec("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
438 TF1 foKaon("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
439 TF1 foMuon("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
440
441 foProton.SetParameters(fAlephParameters);
442 foPion.SetParameters(fAlephParameters);
443 foElec.SetParameters(fAlephParameters);
444 foKaon.SetParameters(fAlephParameters);
445 foMuon.SetParameters(fAlephParameters);
446
c77fc043 447 const Float_t k2sigmaCorr = 1/0.9545;
635c8b48 448
449 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
450
451 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
452 //
453 // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
454 //
455 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
456 if (!track) continue;
457 AliESDtrack *trackDummy = fESD->GetTrack(i);
458 if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
459 delete track;
460 continue;
461 }
462 AliExternalTrackParam trackIn(*trackDummy->GetInnerParam());
463 Double_t ptot = trackIn.GetP(); // momentum for dEdx determination
464 Double_t pT = track->Pt();
465 Double_t eta = TMath::Abs(track->Eta());
466 //
467 fHistEtaNcls->Fill(track->Eta(),track->GetTPCNcls());
468 fHistEtaPhi->Fill(track->Eta(),track->Phi());
469 // cut for dead regions in the detector
470 // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
471 //
472 // 2.a) apply some standard track cuts according to general recommendations
473 //
474
475 if (!fESDtrackCuts->AcceptTrack(track)) {
476 //
477 if (fMCtrue) {
478 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
479 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
480 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(6,eta,track->GetSign()*pT); // Kaon control hist
481 }
482 //
483 delete track;
484 continue;
485 }
486
487 //
488 // 3. make the PID
489 //
490 Double_t sign = track->GetSign();
491
492 Double_t tpcSignal = track->GetTPCsignal();
493 //
494 fDeDx->Fill(ptot, tpcSignal);
495 //
496 //
497 //
498 fHistClassicKaon->Fill(sign*pT,eta,(tpcSignal-foKaon.Eval(ptot))/foKaon.Eval(ptot));
499 fHistClassicProton->Fill(sign*pT,eta,(tpcSignal-foProton.Eval(ptot))/foProton.Eval(ptot));
500 fHistClassicPion->Fill(sign*pT,eta,(tpcSignal-foPion.Eval(ptot))/foPion.Eval(ptot));
501 fHistClassicElectron->Fill(sign*pT,eta,(tpcSignal-foElec.Eval(ptot))/foElec.Eval(ptot));
502 //
503 // fill histograms for dNdpT
504 //
505 if (eta < 0.15) {
506 fHighPtHadrons->Fill(pT);
507 Double_t delta = (tpcSignal - foElec.Eval(ptot))/foElec.Eval(ptot);
508 if (delta > 0) fHighPtElectrons->Fill(pT);
509 }
510 //
511 /* 2sigma PID with 2sigma eff correction! */
10d100d4 512 Double_t tpcMom = track->GetP();
513 const AliExternalTrackParam *in = track->GetTPCInnerParam();
514 if (in)
515 tpcMom = in->GetP();
635c8b48 516 // PION
10d100d4 517 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kPion)) < 2) {
c77fc043 518 fHistPtEtaPion->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
635c8b48 519 //
520 if (trackCounter < 300 && fMCtrue) {
521 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
522 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
c77fc043 523 if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,eta,sign*pT,k2sigmaCorr);
635c8b48 524 if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
c77fc043 525 fHistEffPion->Fill(1,eta,sign*pT,k2sigmaCorr);
635c8b48 526 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
527 Int_t code = TMath::Abs(trackMother->GetPdgCode());
c77fc043 528 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffPion->Fill(3,eta,sign*pT,k2sigmaCorr);
635c8b48 529 }
c77fc043 530 if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,eta,sign*pT,k2sigmaCorr);
531 if (TMath::Abs(trackMC->GetPdgCode()) == 13) fHistEffPion->Fill(4,eta,sign*pT,k2sigmaCorr);
635c8b48 532 }
533 }
534 // KAON
10d100d4 535 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2) {
c77fc043 536 fHistPtEtaKaon->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
635c8b48 537 //
538 if (trackCounter < 300 && fMCtrue) {
539 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
540 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
c77fc043 541 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,eta,sign*pT,k2sigmaCorr);
635c8b48 542 if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
c77fc043 543 fHistEffKaon->Fill(1,eta,sign*pT,k2sigmaCorr);
635c8b48 544 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
545 Int_t code = TMath::Abs(trackMother->GetPdgCode());
c77fc043 546 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffKaon->Fill(3,eta,sign*pT,k2sigmaCorr);
635c8b48 547 }
c77fc043 548 if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,eta,sign*pT,k2sigmaCorr);
635c8b48 549 }
550 }
551 // KAON NO KINK
10d100d4 552 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2 && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
635c8b48 553 // PROTON
10d100d4 554 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)) < 2) {
c77fc043 555 fHistPtEtaProton->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
635c8b48 556 //
557 track->GetImpactParameters(dca,cov);
558 Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
c77fc043 559 fHistPtEtaProtonDCA->Fill(primVtxDCA, eta, sign*pT, k2sigmaCorr);
635c8b48 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 == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
c77fc043 565 fHistEffProton->Fill(0.,eta,sign*pT,k2sigmaCorr);
566 if (eta < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
635c8b48 567 }
568 if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
c77fc043 569 fHistEffProton->Fill(1,eta,sign*pT,k2sigmaCorr);
570 if (eta < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
635c8b48 571 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
572 Int_t code = TMath::Abs(trackMother->GetPdgCode());
573 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
c77fc043 574 fHistEffProton->Fill(3,eta,sign*pT,k2sigmaCorr);
575 if (eta < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
635c8b48 576 }
577 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());
578 }
c77fc043 579 if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,eta,sign*pT,k2sigmaCorr);
635c8b48 580 }
581 }
582 // ELECTRON
10d100d4 583 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron))) fHistPtEtaElectron->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
635c8b48 584
585 delete track;
586
587 } // end of track loop
588
589 fHistTrackPerEvent->Fill(trackCounter);
590
591 // Post output data
592
593
c77fc043 594 PostData(1, fListHist);
635c8b48 595
596}
597
598
599//________________________________________________________________________
600void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
601{
602 // Draw result to the screen
603 // Called once at the end of the query
604 cout << "*** TERMINATE ***" << endl;
605
606}
607
608
609
610
611//________________________________________________________________________
612void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
613
614 // Method for the correct logarithmic binning of histograms
615
616 TAxis *axis = h->GetXaxis();
617 int bins = axis->GetNbins();
618
619 Double_t from = axis->GetXmin();
620 Double_t to = axis->GetXmax();
621 Double_t *newBins = new Double_t[bins + 1];
622
623 newBins[0] = from;
624 Double_t factor = pow(to/from, 1./bins);
625
626 for (int i = 1; i <= bins; i++) {
627 newBins[i] = factor * newBins[i-1];
628 }
629 axis->Set(bins, newBins);
630 delete newBins;
631
632}
633
634
635//________________________________________________________________________
636Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
637 //
638 // get the process type of the event.
639 //
640
641 // can only read pythia headers, either directly or from cocktalil header
642 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
643
644 if (!pythiaGenHeader) {
645
646 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
647 if (!genCocktailHeader) {
648 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
649 return -1;
650 }
651
652 TList* headerList = genCocktailHeader->GetHeaders();
653 if (!headerList) {
654 return -1;
655 }
656
657 for (Int_t i=0; i<headerList->GetEntries(); i++) {
658 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
659 if (pythiaGenHeader)
660 break;
661 }
662
663 if (!pythiaGenHeader) {
664 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
665 return -1;
666 }
667 }
668
669 if (adebug) {
670 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
671 }
672
673 return pythiaGenHeader->ProcessType();
674}
675
676
677//________________________________________________________________________
678TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBin, const Double_t * AlephParams) {
679 //
680 // The multiple Gauss fitting is happening here...
681 // input histogram is of the form (Pt, eta, difference in dEdx)
682 //
683
684 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
685 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
686 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
687 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
688 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
689 //
690 foProton->SetParameters(AlephParams);
691 foPion->SetParameters(AlephParams);
692 foElec->SetParameters(AlephParams);
693 foKaon->SetParameters(AlephParams);
694 foMuon->SetParameters(AlephParams);
695
696
697 const Double_t kSigma = 0.06; // expected resolution (roughly)
698
699 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
700 const Double_t kPtMax = input->GetXaxis()->GetXmax();
701 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
702
703 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
704 // sort pT-bins ..
705 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, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.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, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
706 Int_t indexes[kPtBins+1];
707 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
708 Double_t binsPt[kPtBins+1];
709 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
710 histPt->GetXaxis()->Set(kPtBins, binsPt);
711 //
712
713 TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
714 canvMany->Print("canvManyProton.ps[");
715
716 for(Int_t x=1; x < kPtBins+1; x++) {
717 Double_t pT = input->GetXaxis()->GetBinCenter(x);
718 for(Int_t y=1; y < kEtaBins+1; y++) {
719 Double_t eta = input->GetYaxis()->GetBinCenter(y);
720 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
721 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
722 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
723 histDeDx->SetTitle(Form("pt_%f",pT));
724 Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
725 //
726 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
727 Double_t parameters[12] = {maxYield/2.,0,kSigma,maxYield/2.,(foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma};
728 gausFit->SetParameters(parameters);
729 gausFit->SetParLimits(0,0.,maxYield);
730 gausFit->SetParLimits(1,-0.05,0.05);
731 gausFit->SetParLimits(2,0.04,0.08);
732 //
733 Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
734 gausFit->SetParLimits(0,0.,maxYield);
735 gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
736 gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
737 //
738 Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
739 gausFit->SetParLimits(0,0.,maxYield);
740 gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
741 gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
742 //
743 Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
744 gausFit->SetParLimits(0,0.,maxYield);
745 gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
746 gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
747 histDeDx->Fit("gausFit","QNR");
748 gausFit->GetParameters(parameters);
749 // visualization
750 if (y == EtaBin) {
751 canvMany->cd(x);
752 gPad->SetLogy();
753 histDeDx->SetMarkerStyle(21);
754 histDeDx->SetMarkerSize(0.7);
755 histDeDx->Draw("E");
756 gausFit->SetLineWidth(2);
757 gausFit->Draw("same");
758 //
759 TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
760 TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
761 TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
762 TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
763 gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
764 gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
765 gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
766 gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
767 //
768 gausFit0.SetLineColor(4);
769 gausFit1.SetLineColor(2);
770 gausFit2.SetLineColor(6);
771 gausFit3.SetLineColor(8);
772 //
773 gausFit0.SetLineWidth(1);
774 gausFit1.SetLineWidth(1);
775 gausFit2.SetLineWidth(1);
776 gausFit3.SetLineWidth(1);
777 //
778 gausFit0.Draw("same");
779 gausFit1.Draw("same");
780 gausFit2.Draw("same");
781 gausFit3.Draw("same");
782 canvMany->Print("canvManyProton.ps");
783 }
784
785 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
786 delete gausFit;
787 //
788 // stupid solution --> remove as soon as possible
789 /*yield = 0;
790 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
791 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
792 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
793 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
794 yield += gausYield->Eval(delta);
795 }*/
796 //
797 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
798 //delete gausYield;
799 }
800 }
801 canvMany->Print("canvManyProton.ps]");
802
803 TCanvas * canvPt = new TCanvas();
804 canvPt->cd();
805 histPt->Draw();
806
807 return histPt;
808}
809
810
811
812//________________________________________________________________________
813TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBin,const Double_t * AlephParams) {
814 //
815 // The multiple Gauss fitting is happening here...
816 //
817
818 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
819 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
820 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
821 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
822 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
823 //
824 foProton->SetParameters(AlephParams);
825 foPion->SetParameters(AlephParams);
826 foElec->SetParameters(AlephParams);
827 foKaon->SetParameters(AlephParams);
828 foMuon->SetParameters(AlephParams);
829
830 // input histogram is of the form (Pt, eta, difference in dEdx)
831
832 const Double_t kSigma = 0.06; // expected resolution (roughly)
833
834 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
835 const Double_t kPtMax = input->GetXaxis()->GetXmax();
836 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
837
838 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
839 // sort pT-bins ..
840 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, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.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, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
841 Int_t indexes[kPtBins+1];
842 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
843 Double_t binsPt[kPtBins+1];
844 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
845 histPt->GetXaxis()->Set(kPtBins, binsPt);
846 //
847
848 TCanvas * canvMany = new TCanvas();
849 canvMany->Divide(8,5);
850
851 for(Int_t x=1; x < kPtBins+1; x++) {
852 Double_t pT = input->GetXaxis()->GetBinCenter(x);
853 for(Int_t y=1; y < kEtaBins+1; y++) {
854 Double_t eta = input->GetYaxis()->GetBinCenter(y);
855 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
856 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
857 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
858 //
859 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
860 //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foKaon->Eval(ptot)/foPion->Eval(ptot)),kSigma};
861 Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma};
862 gausFit->SetParameters(parameters);
863 gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
864 gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
865 //
866 gausFit->FixParameter(4,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
867 gausFit->FixParameter(5,kSigma);
868 gausFit->FixParameter(7,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
869 gausFit->FixParameter(8,kSigma);
870 gausFit->FixParameter(10,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
871 gausFit->FixParameter(11,kSigma);
872 histDeDx->Fit("gausFit","Q0");
873 if (y == EtaBin && pT < 0) {
874 canvMany->cd(x);
875 gPad->SetLogy();
876 histDeDx->Draw();
877 gausFit->Draw("same");
878 }
879 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
880 // stupid solution --> remove as soon as possible
881 yield = 0;
882 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
883 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
884 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
885 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
886 yield += gausYield->Eval(delta);
887 }
888 //
889 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
890 //delete gausFit;
891 delete gausYield;
892 }
893 }
894
895 TCanvas * canvPt = new TCanvas();
896 canvPt->cd();
897 histPt->Draw();
898
899 return histPt;
900
901}
902
903
904//________________________________________________________________________
905TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBin,const Double_t * AlephParams) {
906 //
907 // The multiple Gauss fitting is happening here...
908 //
909
910 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
911 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
912 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
913 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
914 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
915 //
916 foProton->SetParameters(AlephParams);
917 foPion->SetParameters(AlephParams);
918 foElec->SetParameters(AlephParams);
919 foKaon->SetParameters(AlephParams);
920 foMuon->SetParameters(AlephParams);
921
922 const Double_t kSigma = 0.06; // expected resolution (roughly)
923
924 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
925 const Double_t kPtMax = input->GetXaxis()->GetXmax();
926 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
927
928 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
929 // sort pT-bins ..
930 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, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.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, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
931 Int_t indexes[kPtBins+1];
932 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
933 Double_t binsPt[kPtBins+1];
934 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
935 histPt->GetXaxis()->Set(kPtBins, binsPt);
936 //
937
938 TCanvas * canvMany = new TCanvas();
939 canvMany->Divide(8,5);
940
941 for(Int_t x=1; x < kPtBins+1; x++) {
942 Double_t pT = input->GetXaxis()->GetBinCenter(x);
943 for(Int_t y=1; y < kEtaBins+1; y++) {
944 Double_t eta = input->GetYaxis()->GetBinCenter(y);
945 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
946 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
947 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
948 //
949 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
950 //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foPion->Eval(ptot)/foPion->Eval(ptot)),kSigma};
951 Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foPion->Eval(ptot)-foKaon->Eval(ptot))/foPion->Eval(ptot),kSigma};
952 gausFit->SetParameters(parameters);
953 gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
954 gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
955 //
956 gausFit->FixParameter(4,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
957 gausFit->FixParameter(5,kSigma);
958 gausFit->FixParameter(7,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
959 gausFit->FixParameter(8,kSigma);
960 gausFit->FixParameter(10,(foPion->Eval(ptot)- foKaon->Eval(ptot))/foKaon->Eval(ptot));
961 gausFit->FixParameter(11,kSigma);
962 histDeDx->Fit("gausFit","Q0");
963 if (y == EtaBin && pT < 0) {
964 canvMany->cd(x);
965 gPad->SetLogy();
966 histDeDx->Draw();
967 gausFit->Draw("same");
968 }
969 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
970 // stupid solution --> remove as soon as possible
971 yield = 0;
972 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
973 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
974 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
975 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
976 yield += gausYield->Eval(delta);
977 }
978 //
979 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
980 //delete gausFit;
981 delete gausYield;
982 }
983 }
984
985 TCanvas * canvPt = new TCanvas();
986 canvPt->cd();
987 histPt->Draw();
988
989 return histPt;
990
991}
992
993
994//________________________________________________________________________
995void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const TList * ListOfHistogramsData, const Char_t *filename) {
996
997 // Define ranges
998 const Int_t kMultLow = 1;
999 const Int_t kMultHigh = 10;
1000
1001 const Int_t kEtaHigh = 4;
1002
1003 // Extract histograms for the MC efficiency study
1004 TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1005 TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1006 TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
635c8b48 1007
1008 TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
635c8b48 1009 TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1010 TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
635c8b48 1011
1012 TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1013 TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1014 TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1015
1016 // Extract data for the final spectra
1017 TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
635c8b48 1018 TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1019 TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
635c8b48 1020
1021 // "MC" of the real data for consistency checks
1022 TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1023 TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1024 TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
635c8b48 1025
1026 TFile spectraFile(filename,"RECREATE");
1027
1028 // 1. Protons
1029 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1030 TH1D *protonRaw = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1031 protonRaw->Sumw2();
1032 TH1D *protonMC = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1033 protonMC->Sumw2();
1034 TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1035 protonRawPrim->Sumw2();
1036 TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1037 protonSecReac->Sumw2();
1038 TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1039 protonSecWeak->Sumw2();
1040 TH1D *protonMis = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1041 protonMis->Sumw2();
1042 protonSecReac->Add(protonSecWeak,-1);
1043 //
1044 // 1. total (physical) efficiency
1045 //
1046 TH1D *protonEffTot = new TH1D(*protonRaw);
1047 protonEffTot->Divide(protonRaw,protonMC);
1048 protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1049 protonEffTot->Write();
1050 //
1051 // 2. contributions from secondaries created by interactions
1052 //
1053 TH1D *protonEffSecReac = new TH1D(*protonRaw);
1054 TH1D *dummy = new TH1D(*protonRaw);
1055 dummy->Add(protonSecReac,-1);
1056 protonEffSecReac->Divide(protonRaw,dummy);
1057 protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1058 protonEffSecReac->Write();
1059 delete dummy;
1060 //
1061 // 3. contributions from secondaries from weak decays
1062 //
1063 TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1064 TH1D *dummy2 = new TH1D(*protonRaw);
1065 dummy2->Add(protonSecWeak,-1);
1066 protonEffSecWeak->Divide(protonRaw,dummy2);
1067 protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1068 protonEffSecWeak->Write();
1069 delete dummy2;
1070 //
1071 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1072 //
1073 TH1D *protonEffAbs = new TH1D(*protonRaw);
1074 protonEffAbs->Divide(protonRawPrim,protonMC);
1075 protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1076 protonEffAbs->Write();
1077 //
1078 // 5. contamination
1079 //
1080 TH1D *protonEffMis = new TH1D(*protonMis);
1081 protonEffMis->Divide(protonMis,protonRaw);
1082 protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1083 protonEffMis->Write();
1084
1085 //
1086 // PROCESS THE REAL DATA FROM HERE ON
1087 //
1088 TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1089 histPtProtonEta->Write();
1090 //
1091 // Apply corrections as you would like them to have...
1092 //
1093 TH1D *protonFinal = new TH1D(*histPtProtonEta);
1094 protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1095 protonFinal->Sumw2();
1096 protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1097 protonFinal->Write();
1098 //
1099 // if available save also the MC truth for consistency checks
1100 //
1101 TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1102 histPtProtonEtaMcData->Write();
1103 }
1104
1105 // 2. Kaons
1106 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1107 TH1D *kaonRaw = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1108 kaonRaw->Sumw2();
1109 TH1D *kaonMC = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1110 kaonMC->Sumw2();
1111 TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1112 kaonRawPrim->Sumw2();
1113 TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1114 kaonSecReac->Sumw2();
1115 TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1116 kaonSecWeak->Sumw2();
1117 TH1D *kaonMis = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1118 kaonMis->Sumw2();
1119 kaonSecReac->Add(kaonSecWeak,-1);
1120 //
1121 // 1. total (physical) efficiency
1122 //
1123 TH1D *kaonEffTot = new TH1D(*kaonRaw);
1124 kaonEffTot->Divide(kaonRaw,kaonMC);
1125 kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1126 kaonEffTot->Write();
1127 //
1128 // 2. contributions from secondaries created by interactions
1129 //
1130 TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1131 TH1D *dummy = new TH1D(*kaonRaw);
1132 dummy->Add(kaonSecReac,-1);
1133 kaonEffSecReac->Divide(kaonRaw,dummy);
1134 kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1135 kaonEffSecReac->Write();
1136 delete dummy;
1137 //
1138 // 3. contributions from secondaries from weak decays
1139 //
1140 TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1141 TH1D *dummy2 = new TH1D(*kaonRaw);
1142 dummy2->Add(kaonSecWeak,-1);
1143 kaonEffSecWeak->Divide(kaonRaw,dummy2);
1144 kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1145 kaonEffSecWeak->Write();
1146 delete dummy2;
1147 //
1148 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1149 //
1150 TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1151 kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1152 kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1153 kaonEffAbs->Write();
1154 //
1155 // 5. contamination
1156 //
1157 TH1D *kaonEffMis = new TH1D(*kaonMis);
1158 kaonEffMis->Divide(kaonMis,kaonRaw);
1159 kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1160 kaonEffMis->Write();
1161
1162 //
1163 // PROCESS THE REAL DATA FROM HERE ON
1164 //
1165 TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1166 histPtKaonEta->Write();
1167 //
1168 // Apply corrections as you would like them to have...
1169 //
1170 TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1171 kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1172 kaonFinal->Sumw2();
1173 kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1174 kaonFinal->Write();
1175 //
1176 // if available save also the MC truth for consistency checks
1177 //
1178 TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1179 histPtKaonEtaMcData->Write();
1180 }
1181
1182
1183 // 3. Pions
1184 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1185 TH1D *pionRaw = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1186 pionRaw->Sumw2();
1187 TH1D *pionMC = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1188 pionMC->Sumw2();
1189 TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1190 pionRawPrim->Sumw2();
1191 TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1192 pionSecReac->Sumw2();
1193 TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1194 pionSecWeak->Sumw2();
1195 TH1D *pionMis = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1196 pionMis->Sumw2();
1197 pionSecReac->Add(pionSecWeak,-1);
1198 //
1199 // 1. total (physical) efficiency
1200 //
1201 TH1D *pionEffTot = new TH1D(*pionRaw);
1202 pionEffTot->Divide(pionRaw,pionMC);
1203 pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1204 pionEffTot->Write();
1205 //
1206 // 2. contributions from secondaries created by interactions
1207 //
1208 TH1D *pionEffSecReac = new TH1D(*pionRaw);
1209 TH1D *dummy = new TH1D(*pionRaw);
1210 dummy->Add(pionSecReac,-1);
1211 pionEffSecReac->Divide(pionRaw,dummy);
1212 pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1213 pionEffSecReac->Write();
1214 delete dummy;
1215 //
1216 // 3. contributions from secondaries from weak decays
1217 //
1218 TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1219 TH1D *dummy2 = new TH1D(*pionRaw);
1220 dummy2->Add(pionSecWeak,-1);
1221 pionEffSecWeak->Divide(pionRaw,dummy2);
1222 pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1223 pionEffSecWeak->Write();
1224 delete dummy2;
1225 //
1226 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1227 //
1228 TH1D *pionEffAbs = new TH1D(*pionRaw);
1229 pionEffAbs->Divide(pionRawPrim,pionMC);
1230 pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1231 pionEffAbs->Write();
1232 //
1233 // 5. contamination
1234 //
1235 TH1D *pionEffMis = new TH1D(*pionMis);
1236 pionEffMis->Divide(pionMis,pionRaw);
1237 pionEffMis->SetName(Form("PionContam_%i",iEta));
1238 pionEffMis->Write();
1239
1240 //
1241 // PROCESS THE REAL DATA FROM HERE ON
1242 //
1243 TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1244 histPtPionEta->Write();
1245 //
1246 // Apply corrections as you would like them to have...
1247 //
1248 TH1D *pionFinal = new TH1D(*histPtPionEta);
1249 pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1250 pionFinal->Sumw2();
1251 pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1252 pionFinal->Write();
1253 //
1254 // if available save also the MC truth for consistency checks
1255 //
1256 TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1257 histPtPionEtaMcData->Write();
1258
1259 }
1260
1261 //
1262 // Close the file after everthing is done
1263 //
1264 spectraFile.Close();
1265
1266}
1267
1268
1269
1270