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