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