]>
Commit | Line | Data |
---|---|---|
88b71b9f | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de> * | |
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 | // The task: | |
17 | // stores TPC PID quantities in a THnSparse | |
18 | // output can then be used for e.g. dEdx calibration | |
19 | // | |
20 | // Author: | |
21 | // Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de> | |
22 | // | |
23 | ||
24 | ||
25 | #include "AliTPCcalibResidualPID.h" | |
26 | #include "TChain.h" | |
27 | #include "AliAnalysisManager.h" | |
28 | #include "AliESDEvent.h" | |
29 | #include "AliAODEvent.h" | |
30 | #include "AliMCEvent.h" | |
31 | #include "AliMCParticle.h" | |
32 | //#include "AliStack.h" | |
33 | #include "AliPID.h" | |
34 | #include "AliTPCcalibResidualPID.h" | |
35 | #include "AliESDtrack.h" | |
36 | #include "AliESDtrackCuts.h" | |
37 | #include "AliESDpid.h" | |
38 | #include "AliESDInputHandler.h" | |
39 | #include "AliESDv0KineCuts.h" | |
40 | #include "AliESDv0.h" | |
41 | #include "AliCentrality.h" | |
42 | #include "THnSparse.h" | |
43 | #include "TH2D.h" | |
44 | #include "TCanvas.h" | |
45 | #include "TGraphErrors.h" | |
46 | #include "TMultiGraph.h" | |
47 | #include "TAxis.h" | |
48 | #include "TFile.h" | |
49 | #include "TSpline.h" | |
50 | #include "TStyle.h" | |
51 | #include "AliTPCdEdxInfo.h" | |
52 | #include "TString.h" | |
53 | #include "TFitResult.h" | |
54 | #include "TH1F.h" | |
55 | #include "TLegend.h" | |
56 | #include "TVirtualFitter.h" | |
57 | ||
58 | using namespace std; | |
59 | ||
60 | ClassImp(AliTPCcalibResidualPID) | |
61 | ||
62 | Double_t AliTPCcalibResidualPID::fgCutGeo = 1.; | |
63 | Double_t AliTPCcalibResidualPID::fgCutNcr = 0.85; | |
64 | Double_t AliTPCcalibResidualPID::fgCutNcl = 0.7; | |
65 | ||
66 | //________________________________________________________________________ | |
67 | AliTPCcalibResidualPID::AliTPCcalibResidualPID() | |
68 | : AliAnalysisTaskSE(), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0), | |
69 | fUseTPCCutMIGeo(kFALSE), | |
70 | fUseMCinfo(kTRUE), | |
71 | fIsPbpOrpPb(kFALSE), | |
72 | fZvtxCutEvent(9999.0), | |
73 | fV0KineCuts(0x0), | |
74 | fNumTagsStored(0), | |
75 | fV0tags(0x0), | |
76 | fV0motherIndex(0x0), | |
77 | fV0motherPDG(0x0), | |
78 | fProduceAllPadTypes(0), | |
79 | fProduceGlobal(0), | |
80 | fProduceShortPads(0), | |
81 | fProduceMediumPads(0), | |
82 | fProduceLongPads(0), | |
83 | fProduceOroc(0), | |
84 | fHistPidQA(0), | |
85 | fHistPidQAshort(0), | |
86 | fHistPidQAmedium(0), | |
87 | fHistPidQAlong(0), | |
88 | fHistPidQAoroc(0), | |
89 | fProduceTPCSignalSparse(0), | |
90 | fCorrectdEdxEtaDependence(0), | |
91 | fCorrectdEdxMultiplicityDependence(0), | |
92 | fThnspTpc(0), | |
93 | fQAList(0x0), | |
94 | fhInvMassGamma(0x0), | |
95 | fhInvMassK0s(0x0), | |
96 | fhInvMassLambda(0x0), | |
97 | fhInvMassAntiLambda(0x0), | |
98 | fhArmenterosAll(0x0), | |
99 | fhArmenterosGamma(0x0), | |
100 | fhArmenterosK0s(0x0), | |
101 | fhArmenterosLambda(0x0), | |
102 | fhArmenterosAntiLambda(0x0), | |
103 | fHistSharedClusQAV0Pi(0x0), | |
104 | fHistSharedClusQAV0Pr(0x0), | |
105 | fHistSharedClusQAV0El(0x0) | |
106 | { | |
107 | // default Constructor | |
108 | /* fast compilation test | |
109 | gSystem->Load("libANALYSIS"); | |
110 | gSystem->Load("libANALYSISalice"); | |
111 | .L /lustre/alice/akalweit/train/trunk/util/statsQA/AliTPCcalibResidualPID.cxx++ | |
112 | */ | |
113 | ||
114 | } | |
115 | ||
116 | ||
117 | //________________________________________________________________________ | |
118 | AliTPCcalibResidualPID::AliTPCcalibResidualPID(const char *name) | |
119 | : AliAnalysisTaskSE(name), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0), | |
120 | fUseTPCCutMIGeo(kFALSE), | |
121 | fUseMCinfo(kTRUE), | |
122 | fIsPbpOrpPb(kFALSE), | |
123 | fZvtxCutEvent(9999.0), | |
124 | fV0KineCuts(0x0), | |
125 | fNumTagsStored(0), | |
126 | fV0tags(0x0), | |
127 | fV0motherIndex(0x0), | |
128 | fV0motherPDG(0x0), | |
129 | fProduceAllPadTypes(0), | |
130 | fProduceGlobal(0), | |
131 | fProduceShortPads(0), | |
132 | fProduceMediumPads(0), | |
133 | fProduceLongPads(0), | |
134 | fProduceOroc(0), | |
135 | fHistPidQA(0), | |
136 | fHistPidQAshort(0), | |
137 | fHistPidQAmedium(0), | |
138 | fHistPidQAlong(0), | |
139 | fHistPidQAoroc(0), | |
140 | fProduceTPCSignalSparse(0), | |
141 | fCorrectdEdxEtaDependence(0), | |
142 | fCorrectdEdxMultiplicityDependence(0), | |
143 | fThnspTpc(0), | |
144 | fQAList(0x0), | |
145 | fhInvMassGamma(0x0), | |
146 | fhInvMassK0s(0x0), | |
147 | fhInvMassLambda(0x0), | |
148 | fhInvMassAntiLambda(0x0), | |
149 | fhArmenterosAll(0x0), | |
150 | fhArmenterosGamma(0x0), | |
151 | fhArmenterosK0s(0x0), | |
152 | fhArmenterosLambda(0x0), | |
153 | fhArmenterosAntiLambda(0x0), | |
154 | fHistSharedClusQAV0Pi(0x0), | |
155 | fHistSharedClusQAV0Pr(0x0), | |
156 | fHistSharedClusQAV0El(0x0) | |
157 | { | |
158 | ||
159 | //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE); | |
160 | //fESDtrackCutsV0 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE); | |
161 | ||
162 | // Constructor | |
163 | DefineInput(0, TChain::Class()); | |
164 | DefineOutput(1, TObjArray::Class()); | |
165 | DefineOutput(2, TObjArray::Class()); | |
166 | ||
167 | } | |
168 | ||
169 | ||
170 | //_________________________________________________ | |
171 | AliTPCcalibResidualPID::~AliTPCcalibResidualPID() | |
172 | { | |
173 | delete fOutputContainer; | |
174 | fOutputContainer = 0; | |
175 | ||
176 | delete fQAList; | |
177 | fQAList = 0; | |
178 | ||
179 | delete fESDtrackCuts; | |
180 | fESDtrackCuts = 0; | |
181 | ||
182 | delete fESDtrackCutsV0; | |
183 | fESDtrackCutsV0 = 0; | |
184 | ||
185 | delete fV0KineCuts; | |
186 | fV0KineCuts = 0; | |
187 | ||
188 | delete fV0tags; | |
189 | fV0tags = 0; | |
190 | fNumTagsStored = 0; | |
191 | ||
192 | delete fV0motherIndex; | |
193 | fV0motherIndex = 0; | |
194 | ||
195 | delete fV0motherPDG; | |
196 | fV0motherPDG = 0; | |
197 | } | |
198 | ||
199 | ||
200 | //________________________________________________________________________ | |
201 | void AliTPCcalibResidualPID::UserCreateOutputObjects() | |
202 | { | |
203 | ||
204 | // THnSparse binning | |
205 | const Int_t kNdim = 9; | |
206 | //v0id, dEdx, TPCsigele, TPCsigpion, TOFbit, eta, TPCclus, centr, p | |
207 | Int_t bins[kNdim] = { 4, 250, 200, 200, 8, 20, 50, 20, 100}; | |
208 | Double_t xmin[kNdim] = { 0.5, 30, -10., -10., -1.5, -1., 60., 0., 0.1}; | |
209 | Double_t xmax[kNdim] = { 4.5, 500., 10., 10., 6.5, 1., 160., 100, 4}; | |
210 | fThnspTpc= new THnSparseF("tpcsignals", "TPC signal;v0id;tpc signal;tpcsigele;tpcsigpion;tofbit;eta;tpcclus;centr;p (GeV/c)", kNdim, bins, xmin, xmax); | |
211 | BinLogAxis(fThnspTpc, 8); | |
212 | ||
213 | // | |
214 | // 0.ptot, 1.tpcSig, 2.particle ID, 3. assumed particle, 4. nSigmaTPC (4x), 5. nSigmaTOF (4x), 6. centrality | |
215 | // Concerning 2 (particle ID): | |
216 | // - in case of MC: Contains MC ID. Bin 1-4 (<=> Slot 0-3): el, pi, ka, pr | |
217 | // - in case of data: Contains V0 particles + bin with all others: Bin 1-4 (<=> Slot 0-3): All non-V0s, V0-el, V0-pi, V0-pr | |
218 | // | |
219 | ||
220 | Int_t binsHistQA[7] = {135, 1980, 4, 5, 40, 10, 40}; | |
221 | Double_t xminHistQA[7] = {0.1, 20, -0.5, -0.5, -10, -5, 0.}; | |
222 | Double_t xmaxHistQA[7] = {50., 2000, 3.5, 4.5, 10, 5, 20000}; | |
223 | fHistPidQA = new THnSparseF("fHistPidQA","PID QA",7,binsHistQA,xminHistQA,xmaxHistQA); | |
224 | BinLogAxis(fHistPidQA, 0); | |
225 | ||
226 | // | |
227 | fHistPidQAshort = new THnSparseF("fHistPidQAshort" ,"PID QA -- short pads",7,binsHistQA,xminHistQA,xmaxHistQA); | |
228 | fHistPidQAmedium = new THnSparseF("fHistPidQAmedium","PID QA -- med pads",7,binsHistQA,xminHistQA,xmaxHistQA); | |
229 | fHistPidQAlong = new THnSparseF("fHistPidQAlong" ,"PID QA -- long pads",7,binsHistQA,xminHistQA,xmaxHistQA); | |
230 | fHistPidQAoroc = new THnSparseF("fHistPidQAoroc" ,"PID QA -- oroc",7,binsHistQA,xminHistQA,xmaxHistQA); | |
231 | BinLogAxis(fHistPidQAshort, 0); | |
232 | BinLogAxis(fHistPidQAmedium, 0); | |
233 | BinLogAxis(fHistPidQAlong, 0); | |
234 | BinLogAxis(fHistPidQAoroc, 0); | |
235 | // | |
236 | fOutputContainer = new TObjArray(2); | |
237 | fOutputContainer->SetName(GetName()); | |
238 | fOutputContainer->SetOwner(); | |
239 | ||
240 | if(fProduceTPCSignalSparse)fOutputContainer->Add(fThnspTpc); | |
241 | if(fProduceGlobal)fOutputContainer->Add(fHistPidQA); | |
242 | if(fProduceAllPadTypes || fProduceShortPads)fOutputContainer->Add(fHistPidQAshort); | |
243 | if(fProduceAllPadTypes || fProduceMediumPads)fOutputContainer->Add(fHistPidQAmedium); | |
244 | if(fProduceAllPadTypes || fProduceLongPads)fOutputContainer->Add(fHistPidQAlong); | |
245 | if(fProduceAllPadTypes || fProduceOroc)fOutputContainer->Add(fHistPidQAoroc); | |
246 | ||
247 | // V0 Kine cuts | |
248 | fV0KineCuts = new AliESDv0KineCuts; | |
249 | fV0KineCuts->SetGammaCutChi2NDF(5.); | |
250 | // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values! | |
251 | Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. }; | |
252 | fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]); | |
253 | ||
254 | ||
255 | fQAList = new TObjArray(4); | |
256 | fQAList->SetName(GetName()); | |
257 | fQAList->SetOwner(); | |
258 | ||
259 | fhInvMassGamma = new TH1F("fhInvMassGamma", "Invariant Mass of gammas; m_{ee} (GeV/#it{c}^{2}); Entries", 200, 0., 0.2); | |
260 | fhInvMassK0s = new TH1F("fhInvMassK0s", "Invariant Mass of K0s; m_{#pi#pi} (GeV/#it{c}^{2}); Entries;", 200, 0.45, 0.55); | |
261 | fhInvMassLambda = new TH1F("fhInvMassLambda", "Invariant Mass of lambdas; m_{p#pi^{-}} (GeV/#it{c}^{2}); Entries", 200, 1.05, 1.15); | |
262 | fhInvMassAntiLambda = new TH1F("fhInvMassAntiLambda", "Invariant Mass of anti-lambdas; m_{#pi^{+}#bar{p}} (GeV/#it{c}^{2}); Entries", | |
263 | 200, 1.05, 1.15); | |
264 | ||
265 | fQAList->Add(fhInvMassGamma); | |
266 | fQAList->Add(fhInvMassK0s); | |
267 | fQAList->Add(fhInvMassLambda); | |
268 | fQAList->Add(fhInvMassAntiLambda); | |
269 | ||
270 | ||
271 | fhArmenterosAll = new TH2F("fhArmenterosAll", | |
272 | "Armenteros plot all V0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", | |
273 | 200, -1., 1., 200, 0., 0.4); | |
274 | fhArmenterosGamma = new TH2F("fhArmenterosGamma", | |
275 | "Armenteros plot Gamma;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", | |
276 | 200, -1., 1., 200, 0., 0.4); | |
277 | fhArmenterosK0s = new TH2F("fhArmenterosK0s", | |
278 | "Armenteros plot K0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", | |
279 | 200, -1., 1., 200, 0., 0.4); | |
280 | fhArmenterosLambda = new TH2F("fhArmenterosLambda", | |
281 | "Armenteros plot lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", | |
282 | 200, -1., 1., 200, 0., 0.4); | |
283 | fhArmenterosAntiLambda = new TH2F("fhArmenterosAntiLambda", | |
284 | "Armenteros plot anti-lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", | |
285 | 200, -1., 1., 200, 0., 0.4); | |
286 | ||
287 | fQAList->Add(fhArmenterosAll); | |
288 | fQAList->Add(fhArmenterosGamma); | |
289 | fQAList->Add(fhArmenterosK0s); | |
290 | fQAList->Add(fhArmenterosLambda); | |
291 | fQAList->Add(fhArmenterosAntiLambda); | |
292 | ||
293 | ||
294 | const Int_t dimQASharedClusters = 4; | |
295 | const TString axisTitles[dimQASharedClusters] = { "#it{p}_{TPC} (GeV/#it{c})", "#it{#Delta'}", "#it{N}_{shared cl}", "Pad row"}; | |
296 | Int_t binsHistQASharedClusters[dimQASharedClusters] = { 100, 100, 160, 160}; | |
297 | Double_t xminHistQASharedClusters[dimQASharedClusters] = { 0.3, 0.5, 0, -1}; | |
298 | Double_t xmaxHistQASharedClusters[dimQASharedClusters] = { 20, 1.5, 160, 159}; | |
299 | ||
300 | fHistSharedClusQAV0Pi = new THnSparseF("fHistSharedClusQAV0Pi" ,"PID QA shared clusters - V0 pi", dimQASharedClusters, | |
301 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); | |
302 | BinLogAxis(fHistSharedClusQAV0Pi, 0); | |
303 | for (Int_t i = 0; i < dimQASharedClusters; i++) | |
304 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); | |
305 | fQAList->Add(fHistSharedClusQAV0Pi); | |
306 | ||
307 | fHistSharedClusQAV0Pr = new THnSparseF("fHistSharedClusQAV0Pr" ,"PID QA shared clusters - V0 pr", dimQASharedClusters, | |
308 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); | |
309 | BinLogAxis(fHistSharedClusQAV0Pr, 0); | |
310 | for (Int_t i = 0; i < dimQASharedClusters; i++) | |
311 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); | |
312 | fQAList->Add(fHistSharedClusQAV0Pr); | |
313 | ||
314 | fHistSharedClusQAV0El = new THnSparseF("fHistSharedClusQAV0El" ,"PID QA shared clusters - V0 el", dimQASharedClusters, | |
315 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); | |
316 | BinLogAxis(fHistSharedClusQAV0El, 0); | |
317 | for (Int_t i = 0; i < dimQASharedClusters; i++) | |
318 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); | |
319 | fQAList->Add(fHistSharedClusQAV0El); | |
320 | ||
321 | PostData(1,fOutputContainer); | |
322 | PostData(2,fQAList); | |
323 | } | |
324 | ||
325 | ||
326 | //_____________________________________________________________________________ | |
327 | void AliTPCcalibResidualPID::UserExec(Option_t *) | |
328 | { | |
329 | //calls the Process function | |
330 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
331 | ||
332 | if (!esdH) { | |
333 | printf("ERROR: Could not get ESDInputHandler \n"); | |
334 | } | |
335 | else fESD = esdH->GetEvent(); | |
336 | ||
337 | // If MC, forward MCevent | |
338 | fMC = dynamic_cast<AliMCEvent*>(MCEvent()); | |
339 | // | |
340 | Process(fESD, fMC); | |
341 | // | |
342 | PostData(1,fOutputContainer); | |
343 | PostData(2,fQAList); | |
344 | } | |
345 | ||
346 | ||
347 | ||
348 | //________________________________________________________________________ | |
349 | void AliTPCcalibResidualPID::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent) | |
350 | { | |
351 | ||
352 | //called for each event | |
353 | if (!esdEvent) { | |
354 | Printf("ERROR: esdEvent not available"); | |
355 | return; | |
356 | } | |
357 | ||
358 | if(!((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))) printf("ESD inputhandler not available \n"); | |
359 | ||
360 | if (!fESDpid) { | |
361 | fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid(); | |
362 | } | |
363 | ||
364 | if (!fESDpid || !fV0KineCuts) | |
365 | return; | |
366 | ||
367 | ||
368 | Float_t centralityFper=99; | |
369 | ||
370 | AliCentrality *esdCentrality = esdEvent->GetCentrality(); | |
371 | centralityFper = esdCentrality->GetCentralityPercentile("V0M"); | |
372 | ||
373 | if (!GetVertexIsOk(esdEvent)) | |
374 | return; | |
375 | ||
376 | const AliESDVertex* fESDvertex = esdEvent->GetPrimaryVertexTracks(); | |
377 | if (!fESDvertex) | |
378 | return; | |
379 | ||
380 | Int_t ncontr = fESDvertex->GetNContributors(); | |
381 | ||
382 | if (ncontr <= 0) | |
383 | return; | |
384 | ||
385 | // Fill V0 arrays with V0s | |
386 | FillV0PIDlist(esdEvent); | |
387 | ||
388 | // Array with flags wheter QA for this V0 was already done or not | |
389 | const Int_t numV0s = esdEvent->GetNumberOfV0s(); | |
390 | Bool_t v0QAadded[numV0s]; | |
391 | for (Int_t i = 0; i < numV0s; i++) | |
392 | v0QAadded[i] = kFALSE; | |
393 | ||
394 | Int_t nTotTracks = esdEvent->GetNumberOfTracks(); | |
395 | for (Int_t iTracks = 0; iTracks < nTotTracks; iTracks++){//begin track loop | |
396 | AliESDtrack *trackESD = esdEvent->GetTrack(iTracks); | |
397 | if(!trackESD) { | |
398 | Printf("ERROR: Could not receive track %d (esd loop)", iTracks); | |
399 | continue; | |
400 | } | |
401 | if((TMath::Abs(trackESD->Eta())) > 0.9) | |
402 | continue; | |
403 | ||
404 | // Do not distinguish positively and negatively charged V0's | |
405 | Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks)); | |
406 | if (v0tagAllCharges == -99) { | |
407 | AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, esdEvent->GetNumberOfTracks(), | |
408 | fV0tags != 0x0)); | |
409 | continue; | |
410 | } | |
411 | ||
412 | Bool_t isV0el = v0tagAllCharges == 14; | |
413 | Bool_t isV0pi = v0tagAllCharges == 15; | |
414 | Bool_t isV0pr = v0tagAllCharges == 16; | |
415 | Bool_t isV0 = isV0el || isV0pi || isV0pr; | |
416 | ||
417 | if (mcEvent && fUseMCinfo) { | |
418 | // For MC, do not look for V0s, i.e. demand the non-V0 track cuts | |
419 | if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue; | |
420 | ||
421 | if (fUseTPCCutMIGeo) { | |
422 | // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included. | |
423 | if (!TPCCutMIGeo(trackESD, esdEvent)) | |
424 | continue; | |
425 | } | |
426 | else { | |
427 | // If cut on geometry is not active, always cut on num clusters | |
428 | if (trackESD->GetTPCsignalN() < 60) | |
429 | continue; | |
430 | } | |
431 | } | |
432 | else { | |
433 | // For data, take V0 AND non-V0 with separate cuts | |
434 | //if (!isV0 && fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue; | |
435 | if (isV0) { | |
436 | if (fESDtrackCutsV0 && !fESDtrackCutsV0->AcceptTrack(trackESD)) | |
437 | continue; | |
438 | } | |
439 | else { | |
440 | if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) | |
441 | continue; | |
442 | } | |
443 | ||
444 | if (fUseTPCCutMIGeo) { | |
445 | // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included. | |
446 | if (!isV0) { | |
447 | if (!TPCCutMIGeo(trackESD, esdEvent)) | |
448 | continue; | |
449 | } | |
450 | else { | |
451 | // One should not cut on geometry for V0's since they have different topology. (Loosely) cut on num of clusters instead. | |
452 | if (trackESD->GetTPCsignalN() < 60) | |
453 | continue; | |
454 | } | |
455 | } | |
456 | else { | |
457 | // If cut on geometry is not active, always cut on num clusters | |
458 | if (trackESD->GetTPCsignalN() < 60) | |
459 | continue; | |
460 | } | |
461 | } | |
462 | ||
463 | const AliExternalTrackParam *paramIn = trackESD->GetInnerParam(); | |
464 | Float_t precin=-1; | |
465 | if(paramIn) { | |
466 | precin=paramIn->GetP(); | |
467 | } else { | |
468 | continue; | |
469 | } | |
470 | Int_t precdefault=CompareFloat(precin,-1); | |
471 | if(precdefault==1) continue; // momentum cut | |
472 | // | |
473 | ||
474 | AliMCParticle* mcTrack = 0x0; | |
475 | Int_t particleID = -1; | |
476 | Int_t v0id = 0; | |
477 | ||
478 | if (mcEvent && fUseMCinfo) { | |
479 | // MC - particle ID = MC ID | |
480 | Int_t label = trackESD->GetLabel(); | |
481 | ||
482 | if (label < 0) | |
483 | continue; // If MC is available, reject tracks with negative ESD label | |
484 | ||
485 | mcTrack = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(TMath::Abs(label))); | |
486 | if (!mcTrack) { | |
487 | Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks); | |
488 | continue; | |
489 | } | |
490 | ||
491 | /*// Only accept MC primaries | |
492 | if (!mcEvent->Stack()->IsPhysicalPrimary(TMath::Abs(label))) { | |
493 | continue; | |
494 | }*/ | |
495 | ||
496 | Int_t pdgAbs = TMath::Abs(mcTrack->PdgCode()); | |
497 | ||
498 | if (pdgAbs == 11) { // electron | |
499 | particleID = 0; | |
500 | } | |
501 | else if (pdgAbs == 211) { // pion | |
502 | particleID = 1; | |
503 | } | |
504 | else if (pdgAbs == 321) { // kaon | |
505 | particleID = 2; | |
506 | } | |
507 | else if (pdgAbs == 2212) { // proton | |
508 | particleID = 3; | |
509 | } | |
510 | else | |
511 | continue; // Reject all other particles | |
512 | } | |
513 | else { | |
514 | // Data - particle ID = V0 ID (if V0) | |
515 | if (isV0pi) { // pion | |
516 | particleID = 2; | |
517 | v0id = 2; | |
518 | } | |
519 | else if (isV0el) { // electron | |
520 | particleID = 1; | |
521 | v0id = 1; | |
522 | } | |
523 | else if (isV0pr) { // proton | |
524 | particleID = 3; | |
525 | v0id = 3; | |
526 | } | |
527 | else { // No V0-ID available (species must be determined via TPC, TOF, ...) | |
528 | particleID = 0; | |
529 | } | |
530 | } | |
531 | ||
532 | ||
533 | Float_t tpcsignal=trackESD->GetTPCsignal(); | |
534 | Int_t tofbit=-1; | |
535 | Int_t iTOFpid=0; | |
536 | Int_t ikTIME=0; | |
537 | Double_t tpcnsigmaele=-10; | |
538 | Double_t tpcnsigmapion=-10; | |
539 | // | |
540 | if ((trackESD->GetStatus() & AliESDtrack::kTOFpid)) iTOFpid = 1; | |
541 | if ((trackESD->GetStatus() & AliESDtrack::kTIME)) ikTIME = 1; | |
542 | Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero(); | |
543 | // | |
544 | if((iTOFpid==1) &&(ikTIME==1)){ | |
545 | Double_t tofnsigmaele= fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kElectron, time0); | |
546 | Double_t tofnsigmapion=fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kPion, time0); | |
547 | if (TMath::Abs(tofnsigmapion)<3) tofbit = 1; | |
548 | if (TMath::Abs(tofnsigmaele)<3) tofbit = 0; | |
549 | tpcnsigmaele=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kElectron); | |
550 | tpcnsigmapion=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kPion); | |
551 | } | |
552 | // | |
553 | Int_t tpcnclusN=trackESD->GetTPCsignalN(); | |
554 | Double_t eta=trackESD->Eta(); | |
555 | ||
556 | // | |
557 | if(fProduceTPCSignalSparse){ | |
558 | Double_t contentSignal[9]; | |
559 | contentSignal[0]=v0id; | |
560 | contentSignal[1]=tpcsignal; | |
561 | contentSignal[2]=tpcnsigmaele; | |
562 | contentSignal[3]=tpcnsigmapion; | |
563 | contentSignal[4]=tofbit; | |
564 | contentSignal[5]=eta; | |
565 | contentSignal[6]=tpcnclusN; | |
566 | contentSignal[7]=centralityFper; | |
567 | contentSignal[8]=precin; | |
568 | // | |
569 | if (fThnspTpc->GetEntries() < 1e6) fThnspTpc->Fill(contentSignal); | |
570 | }//should it be created or not | |
571 | ||
572 | ||
573 | ||
574 | // | |
575 | // 2nd THnSparse | |
576 | // | |
577 | Double_t tpcQA[5] = {fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kElectron), | |
578 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kPion), | |
579 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kKaon), | |
580 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kProton), | |
581 | 0}; | |
582 | Double_t tofQA[5] = {fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kElectron, time0), | |
583 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kPion, time0), | |
584 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kKaon, time0), | |
585 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kProton, time0), | |
586 | 0 }; | |
587 | ||
588 | // id 5 is just again Kaons in restricted eta range | |
589 | tpcQA[4] = tpcQA[2]; | |
590 | tofQA[4] = tofQA[2]; | |
591 | ||
592 | // | |
593 | // dE/dx in different pad regions | |
594 | // | |
595 | AliTPCdEdxInfo * infoTpcPid = trackESD->GetTPCdEdxInfo(); | |
596 | Double32_t signal[4]; Char_t ncl[3]; Char_t nrows[3]; | |
597 | if (infoTpcPid) { | |
598 | infoTpcPid->GetTPCSignalRegionInfo(signal, ncl, nrows); | |
599 | } else { | |
600 | for(Int_t iarr = 0; iarr < 3; iarr++) { | |
601 | signal[iarr] = 0; | |
602 | ncl[iarr] = 0; | |
603 | nrows[iarr] = 0; | |
604 | } | |
605 | signal[3] = 0; | |
606 | } | |
607 | // | |
608 | UInt_t status = trackESD->GetStatus(); | |
609 | Bool_t hasTOFout = status&AliESDtrack::kTOFout; | |
610 | Bool_t hasTOFtime = status&AliESDtrack::kTIME; | |
611 | Bool_t hasTOF = kFALSE; | |
612 | if (hasTOFout && hasTOFtime) hasTOF = kTRUE; | |
613 | Float_t length = trackESD->GetIntegratedLength(); | |
614 | // Length check only for primaries! | |
615 | if (length < 350 && !isV0) hasTOF = kFALSE; | |
616 | // | |
617 | ||
618 | if (!hasTOF) { | |
619 | // Make sure that number of sigmas is large in this case, so that the track will be rejected if a TOF cut is applied | |
620 | for (Int_t i = 0; i < 5; i++) { | |
621 | tofQA[i] = 999; | |
622 | } | |
623 | } | |
624 | ||
625 | ||
626 | Double_t processedTPCsignal[5] = { tpcsignal, tpcsignal, tpcsignal, tpcsignal, tpcsignal }; | |
627 | ||
628 | ||
629 | if (fCorrectdEdxEtaDependence && fCorrectdEdxMultiplicityDependence) { | |
630 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron, | |
631 | AliTPCPIDResponse::kdEdxDefault); | |
632 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion, | |
633 | AliTPCPIDResponse::kdEdxDefault); | |
634 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon, | |
635 | AliTPCPIDResponse::kdEdxDefault); | |
636 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton, | |
637 | AliTPCPIDResponse::kdEdxDefault); | |
638 | } | |
639 | else if (fCorrectdEdxEtaDependence) { | |
640 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault); | |
641 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault); | |
642 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault); | |
643 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault); | |
644 | } | |
645 | else if (fCorrectdEdxMultiplicityDependence) { | |
646 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault); | |
647 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault); | |
648 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault); | |
649 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault); | |
650 | } | |
651 | ||
652 | // id 5 is just again Kaons in restricted eta range | |
653 | processedTPCsignal[4] = processedTPCsignal[2]; | |
654 | ||
655 | for(Int_t iPart = 0; iPart < 5; iPart++) { | |
656 | // Only accept "Kaons" within |eta| < 0.2 for index 4 in case of data (no contamination in case of MC, so this index is not used) | |
657 | if (iPart == 4 && ((mcEvent && fUseMCinfo) || abs(trackESD->Eta()) > 0.2)) | |
658 | continue; | |
659 | ||
0e60aeb1 | 660 | Double_t vecHistQA[7] = {precin, processedTPCsignal[iPart], (Double_t)particleID, (Double_t)iPart, tpcQA[iPart], tofQA[iPart], |
661 | (Double_t)nTotTracks}; | |
88b71b9f | 662 | if (fProduceGlobal) fHistPidQA->Fill(vecHistQA); |
663 | vecHistQA[1] = signal[0]; vecHistQA[4] = ncl[0]; | |
664 | if ((fProduceAllPadTypes || fProduceShortPads)) fHistPidQAshort->Fill(vecHistQA); | |
665 | // | |
666 | vecHistQA[1] = signal[1]; vecHistQA[4] = ncl[1]; | |
667 | if ((fProduceAllPadTypes || fProduceMediumPads)) fHistPidQAmedium->Fill(vecHistQA); | |
668 | // | |
669 | vecHistQA[1] = signal[2]; vecHistQA[4] = ncl[2]; | |
670 | if ((fProduceAllPadTypes || fProduceLongPads)) fHistPidQAlong->Fill(vecHistQA); | |
671 | // | |
672 | vecHistQA[1] = signal[3]; vecHistQA[4] = nrows[1] + nrows[2]; | |
673 | if ((fProduceAllPadTypes || fProduceOroc)) fHistPidQAoroc->Fill(vecHistQA); | |
674 | } | |
675 | ||
676 | if (isV0) { | |
677 | // Fill QA hists for shared clusters dE/dx | |
0e60aeb1 | 678 | Double_t vecHistQA[4] = { precin, -1., (Double_t)trackESD->GetTPCSharedMap().CountBits(), -1. }; |
88b71b9f | 679 | THnSparse* currHist = fHistSharedClusQAV0Pi; |
680 | ||
681 | if (isV0pi) { | |
682 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, | |
683 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); | |
684 | if (expectedDeDx > 0 ) { | |
685 | vecHistQA[1] = tpcsignal / expectedDeDx; | |
686 | currHist = fHistSharedClusQAV0Pi; | |
687 | } | |
688 | } | |
689 | else if (isV0pr) { | |
690 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, | |
691 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); | |
692 | if (expectedDeDx > 0 ) { | |
693 | vecHistQA[1] = tpcsignal / expectedDeDx; | |
694 | currHist = fHistSharedClusQAV0Pr; | |
695 | } | |
696 | } | |
697 | else if (isV0el) { | |
698 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, | |
699 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); | |
700 | if (expectedDeDx > 0 ) { | |
701 | vecHistQA[1] = tpcsignal / expectedDeDx; | |
702 | currHist = fHistSharedClusQAV0El; | |
703 | } | |
704 | } | |
705 | ||
706 | Int_t iRowInd = -1; | |
707 | // iRowInd == -1 for "all rows w/o multiple counting" | |
708 | currHist->Fill(vecHistQA); | |
709 | ||
710 | // Fill hist for every pad row with shared cluster | |
711 | for (iRowInd = 0; iRowInd < 159; iRowInd++) { | |
712 | if (trackESD->GetTPCSharedMap().TestBitNumber(iRowInd)) { | |
713 | vecHistQA[3] = iRowInd; | |
714 | currHist->Fill(vecHistQA); | |
715 | } | |
716 | } | |
717 | ||
718 | ||
719 | ||
720 | // Check, whether the QA for this V0 mother was already filled into the histograms (is the case if the other daughter has already | |
721 | // been processed). If not, set flag to kTRUE and fill the QA. | |
722 | const Int_t iV0 = GetV0motherIndex(iTracks); | |
723 | if (!v0QAadded[iV0]) { | |
724 | v0QAadded[iV0] = kTRUE; | |
725 | ||
726 | AliESDv0* esdV0 = (AliESDv0*)esdEvent->GetV0(iV0); | |
727 | ||
728 | if (!esdV0) { | |
729 | printf("Error: V0 tagged, but does not exist!\n"); | |
730 | continue; | |
731 | } | |
732 | ||
733 | Float_t armVar[2] = {0.0,0.0}; | |
734 | fV0KineCuts->Armenteros(esdV0, armVar); | |
735 | ||
736 | const Int_t motherPDG = GetV0motherPDG(iTracks); | |
737 | ||
738 | // Mother and daughter match the requirements, otherwise it wouldn't be tagged as a V0. Now just fill the QA histos. | |
739 | fhArmenterosAll->Fill(armVar[0], armVar[1]); | |
740 | //if (esdV0->TestBit(BIT(14))) { | |
741 | if (TMath::Abs(motherPDG) == 22) { | |
742 | fhInvMassGamma->Fill(esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron)); | |
743 | fhArmenterosGamma->Fill(armVar[0], armVar[1]); | |
744 | } | |
745 | else if (TMath::Abs(motherPDG) == 310) { | |
746 | //else if (esdV0->TestBit(BIT(15))) { | |
747 | fhInvMassK0s->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kPion)); | |
748 | fhArmenterosK0s->Fill(armVar[0], armVar[1]); | |
749 | } | |
750 | else if (motherPDG == 3122) { | |
751 | //else if (esdV0->TestBit(BIT(16))) { | |
752 | fhInvMassLambda->Fill(esdV0->GetEffMass(AliPID::kProton, AliPID::kPion)); | |
753 | fhArmenterosLambda->Fill(armVar[0], armVar[1]); | |
754 | } | |
755 | else if (motherPDG == -3122) { | |
756 | //else if (esdV0->TestBit(BIT(17))) { | |
757 | fhInvMassAntiLambda->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kProton)); | |
758 | fhArmenterosAntiLambda->Fill(armVar[0], armVar[1]); | |
759 | } | |
760 | } | |
761 | } | |
762 | } //end track loop | |
763 | ||
764 | ||
765 | // Clear the V0 PID arrays | |
766 | ClearV0PIDlist(); | |
767 | } | |
768 | ||
769 | ||
770 | //________________________________________________________________________ | |
771 | Int_t AliTPCcalibResidualPID::CompareFloat(Float_t f1, Float_t f2) const | |
772 | { | |
773 | //compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false | |
774 | Float_t precision = 0.00001; | |
775 | if (((f1 - precision) < f2) && | |
776 | ((f1 + precision) > f2)) | |
777 | { | |
778 | return 1; | |
779 | } | |
780 | else | |
781 | { | |
782 | return 0; | |
783 | } | |
784 | } | |
785 | ||
786 | ||
787 | //________________________________________________________________________ | |
788 | void AliTPCcalibResidualPID::Terminate(const Option_t *) | |
789 | { | |
790 | ||
791 | ||
792 | } | |
793 | ||
794 | ||
795 | ||
796 | //________________________________________________________________________ | |
797 | void AliTPCcalibResidualPID::BinLogAxis(const THnSparseF *h, Int_t axisNumber) { | |
798 | // | |
799 | // Method for the correct logarithmic binning of histograms | |
800 | // | |
801 | TAxis *axis = h->GetAxis(axisNumber); | |
802 | int bins = axis->GetNbins(); | |
803 | ||
804 | Double_t from = axis->GetXmin(); | |
805 | Double_t to = axis->GetXmax(); | |
806 | Double_t *newBins = new Double_t[bins + 1]; | |
807 | ||
808 | newBins[0] = from; | |
809 | Double_t factor = pow(to/from, 1./bins); | |
810 | ||
811 | for (int i = 1; i <= bins; i++) { | |
812 | newBins[i] = factor * newBins[i-1]; | |
813 | } | |
814 | axis->Set(bins, newBins); | |
815 | delete [] newBins; | |
816 | ||
817 | } | |
818 | ||
819 | ||
820 | //________________________________________________________________________ | |
821 | Double_t* AliTPCcalibResidualPID::ExtractResidualPID(THnSparseF * histPidQA, const Bool_t useV0s, const Char_t * outFile, | |
822 | const Char_t * type, const Char_t * period, const Char_t * pass, | |
823 | const Char_t * system, const Double_t * initialParameters, | |
824 | const Char_t *dedxtype, | |
825 | AliTPCcalibResidualPID::FitType fitType) { | |
826 | // | |
827 | // (1.) obtain residual graphs | |
828 | // | |
829 | TObjArray * arr = 0x0; | |
830 | ||
831 | Bool_t isMC = kFALSE; | |
832 | if (strcmp(type, "MC") == 0) { | |
833 | arr = GetResidualGraphsMC(histPidQA, system); | |
834 | isMC = kTRUE; | |
835 | } | |
836 | else if (strcmp(type, "DATA") == 0) { | |
837 | arr = GetResidualGraphs(histPidQA, system, useV0s); | |
838 | } | |
839 | else { | |
840 | Printf("ERROR - ExtractResidualPID: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type); | |
841 | ||
842 | return 0x0; | |
843 | } | |
844 | ||
845 | Bool_t isPPb = kFALSE; | |
846 | if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { | |
847 | Printf("p-Pb/Pb-p detected - Applying special handling!"); | |
848 | isPPb = kTRUE; | |
849 | } | |
850 | // | |
851 | // (2.) get old-style Bethe-Bloch parameters | |
852 | // | |
853 | TF1* parametrisation = FitBB(arr, isMC, isPPb, useV0s, initialParameters, fitType); | |
854 | // | |
855 | // (3.) obtain response functions to OADB | |
856 | // | |
857 | TObjArray* arrResponse = GetResponseFunctions(parametrisation, arr, type, period, pass, system, dedxtype); | |
858 | // | |
859 | // (4.) write results to file and exit | |
860 | // | |
861 | if (arrResponse) { | |
862 | TFile outputFile(outFile,"RECREATE"); | |
863 | arrResponse->Write(); | |
864 | outputFile.Close(); | |
865 | } | |
866 | // | |
867 | ||
868 | return parametrisation->GetParameters(); | |
869 | } | |
870 | ||
871 | ||
872 | //________________________________________________________________________ | |
873 | TObjArray * AliTPCcalibResidualPID::GetSeparation(THnSparseF * histPidQA, Int_t kParticle1, Int_t kParticle2) { | |
874 | ||
875 | //get THnSparse | |
876 | THnSparse * hist = histPidQA; | |
877 | ||
878 | ||
879 | Float_t nSigmaPlus1 = 0, nSigmaMinus1 = 0; //sigma of first particle in the TPC | |
880 | Float_t nSigmaPlus2 = 0, nSigmaMinus2 = 0; //sigma of second particle in the TPC | |
881 | ||
882 | if(kParticle1 == 0){ // electron | |
883 | nSigmaMinus1 = -1.999; | |
884 | nSigmaPlus1 = 2.999; | |
885 | } else if(kParticle1 == 1){ //pion | |
886 | nSigmaMinus1 = -1.999; | |
887 | nSigmaPlus1 = 2.999; | |
888 | } else if(kParticle1 == 2){ //kaons | |
889 | nSigmaMinus1 = -2.999; | |
890 | nSigmaPlus1 = 2.999; | |
891 | } else if(kParticle1 == 3){ //protons | |
892 | nSigmaMinus1 = -2.999; | |
893 | nSigmaPlus1 = 2.999; | |
894 | } | |
895 | ||
896 | // | |
897 | // 1. select and fit first particle | |
898 | // | |
899 | hist->GetAxis(3)->SetRangeUser(kParticle1,kParticle1); | |
900 | hist->GetAxis(4)->SetRangeUser(nSigmaMinus1,nSigmaPlus1); | |
901 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
902 | TH2D * histPart1 = (TH2D*) hist->Projection(1,0)->Clone("histPart1"); | |
903 | histPart1->SetTitle(";p /(GeV/c) ; TPCSignal /(a.u.)"); | |
904 | histPart1->RebinX(4); | |
905 | ||
906 | TObjArray arr; | |
907 | histPart1->FitSlicesY(0,0,-1,10,"QNR",&arr); | |
908 | TH1D * PartMean1 = (TH1D *) arr.At(1)->Clone("PartMean1"); | |
909 | TH1D * PartSigma1 = (TH1D *) arr.At(2)->Clone("PartSigma1"); | |
910 | PartMean1->SetTitle(";p /(GeV/c) ; Mean (Gauss)"); | |
911 | PartSigma1->SetTitle(";p /(GeV/c) ; Sigma (Gauss)"); | |
912 | ||
913 | histPart1->SetMarkerStyle(22); | |
914 | PartMean1->SetMarkerStyle(21); | |
915 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
916 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
917 | ||
918 | if(kParticle2==0){ // electron | |
919 | nSigmaMinus2 = -1.999; | |
920 | nSigmaPlus2 = 2.999; | |
921 | } else if(kParticle2 == 1){ //pion | |
922 | nSigmaMinus2 = -1.999; | |
923 | nSigmaPlus2 = 2.999; | |
924 | } else if(kParticle2 == 2){ //kaons | |
925 | nSigmaMinus2 = -2.999; | |
926 | nSigmaPlus2 = 2.999; | |
927 | } else if(kParticle2 == 3){ //protons | |
928 | nSigmaMinus2 = -2.999; | |
929 | nSigmaPlus2 = 2.999; | |
930 | } | |
931 | ||
932 | // | |
933 | // 2. select and fit second particle | |
934 | // | |
935 | hist->GetAxis(3)->SetRangeUser(kParticle2,kParticle2); | |
936 | hist->GetAxis(4)->SetRangeUser(nSigmaMinus2,nSigmaPlus2); | |
937 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
938 | TH2D * histPart2 = (TH2D*)hist->Projection(1,0)->Clone("histPart2"); | |
939 | histPart2->RebinX(4); | |
940 | histPart2->FitSlicesY(0,0,-1,10,"QNR",&arr); | |
941 | TH1D * PartMean2 = (TH1D *) arr.At(1)->Clone("PartMean2"); | |
942 | TH1D * PartSigma2 = (TH1D *) arr.At(2)->Clone("PartSigma2"); | |
943 | PartMean2->SetTitle(";p /(GeV/c) ; Mean (Gauss)"); | |
944 | PartSigma2->SetTitle(";p /(GeV/c) ; Sigma (Gauss)"); | |
945 | PartMean2->SetMarkerStyle(20); | |
946 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
947 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
948 | ||
949 | // | |
950 | // 3. separation | |
951 | // | |
952 | ||
953 | TH1F *fHistSeparation=(TH1F*)PartMean1->Clone("fHistSeparation"); //to get same binning | |
954 | fHistSeparation->SetMarkerStyle(22); | |
955 | const Int_t Nbins = PartMean1->GetNbinsX(); | |
956 | fHistSeparation->SetTitle(";p /(GeV/c) ; Separation"); | |
957 | ||
958 | Float_t DeltaMean[Nbins] ; | |
959 | Float_t DeltaSigma[Nbins]; | |
960 | ||
961 | for(Int_t i=0 ; i<Nbins; i++){ | |
962 | ||
963 | DeltaMean[i] = TMath::Abs(PartMean1->GetBinContent(i) - PartMean2->GetBinContent(i)); | |
964 | DeltaSigma[i] = TMath::Abs((PartSigma1->GetBinContent(i) + PartSigma2->GetBinContent(i)))/2.; | |
965 | ||
966 | if(!(TMath::Abs(DeltaSigma[i])<0.000001))fHistSeparation->SetBinContent(i,DeltaMean[i]/DeltaSigma[i]); | |
967 | }//for(Int_t i=0 ; i<Nbins ; i++) | |
968 | ||
969 | TObjArray *array = new TObjArray(); | |
970 | array->Add(histPart1); | |
971 | array->Add(histPart2); | |
972 | array->Add(PartMean1); | |
973 | array->Add(PartSigma1); | |
974 | array->Add(PartMean2); | |
975 | array->Add(PartSigma2); | |
976 | array->Add(fHistSeparation); | |
977 | return array; | |
978 | ||
979 | } | |
980 | ||
981 | ||
982 | //________________________________________________________________________ | |
983 | TObjArray * AliTPCcalibResidualPID::GetResidualGraphs(THnSparseF * histPidQA, const Char_t * system, Bool_t useV0s) { | |
984 | // | |
985 | // Extracts the residual graphs from THnSparse created from data (NOT MC!) | |
986 | // | |
987 | ||
988 | const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})"; | |
989 | const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)"; | |
990 | ||
991 | Bool_t isPPb = kFALSE; | |
992 | if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { | |
993 | Printf("p-Pb/Pb-p detected - Applying special handling!"); | |
994 | isPPb = kTRUE; | |
995 | } | |
996 | ||
997 | // the following parameters have to be extracted from the data itself | |
998 | // | |
999 | ||
1000 | Float_t nSigmaPionMinus = -3.999; | |
1001 | Float_t nSigmaPionPlus = 3.999; | |
1002 | Float_t nSigmaKaonMinus = -4.199; | |
1003 | Float_t nSigmaKaonPlus = 4.7999; | |
1004 | Float_t nSigmaElectronMinus = -1.999; | |
1005 | Float_t nSigmaElectronPlus = 4.999; | |
1006 | ||
1007 | Int_t cutForFitting = 10; | |
1008 | Double_t heightFractionForFittingRange = 0.1; | |
1009 | // | |
1010 | THnSparse * hist = histPidQA; | |
1011 | // | |
1012 | TCanvas * canvasQAtpc = new TCanvas("canvasQAtpcResGraph","Control canvas for residual graphs (TPC)",100,10,1380,800); | |
1013 | canvasQAtpc->Divide(2,2); | |
1014 | ||
1015 | TCanvas * canvasQAtof = new TCanvas("canvasQAtofResGraph","Control canvas for residual graphs (TOF)",100,10,1380,800); | |
1016 | canvasQAtof->Divide(2,2); | |
1017 | ||
1018 | TCanvas * canvasQAv0 = new TCanvas("canvasQAv0ResGraph","Control canvas for residual graphs (V0)",100,10,1380,800); | |
1019 | canvasQAv0->Divide(2,2); | |
1020 | ||
1021 | TCanvas * canvasQAv0plusTOF = new TCanvas("canvasQAv0plusTOFResGraph","Control canvas for residual graphs (V0+TOF)",100,10,1380,800); | |
1022 | canvasQAv0plusTOF->Divide(2,2); | |
1023 | ||
1024 | ||
1025 | TCanvas * canvasQAv0DeDxPurityEl = new TCanvas("canvasQAv0DeDxPurityEl","Control canvas for residual graphs (V0 el)",100,10,600,400); | |
1026 | TCanvas * canvasQAv0DeDxPurityPi = new TCanvas("canvasQAv0DeDxPurityPi","Control canvas for residual graphs (V0 pi)",100,10,600,400); | |
1027 | TCanvas * canvasQAv0DeDxPurityPr = new TCanvas("canvasQAv0DeDxPurityPr","Control canvas for residual graphs (V0 pr)",100,10,600,400); | |
1028 | ||
1029 | canvasQAv0DeDxPurityEl->SetLogx(); | |
1030 | canvasQAv0DeDxPurityPi->SetLogx(); | |
1031 | canvasQAv0DeDxPurityPr->SetLogx(); | |
1032 | ||
1033 | canvasQAv0DeDxPurityEl->SetLogz(); | |
1034 | canvasQAv0DeDxPurityPi->SetLogz(); | |
1035 | canvasQAv0DeDxPurityPr->SetLogz(); | |
1036 | ||
1037 | canvasQAv0DeDxPurityEl->SetTopMargin(0.03); | |
1038 | canvasQAv0DeDxPurityPi->SetTopMargin(0.03); | |
1039 | canvasQAv0DeDxPurityPr->SetTopMargin(0.03); | |
1040 | ||
1041 | for (Int_t i = 1; i <= 4; i++) { | |
1042 | canvasQAtpc->GetPad(i)->SetGrid(1, 1); | |
1043 | canvasQAtof->GetPad(i)->SetGrid(1, 1); | |
1044 | canvasQAv0->GetPad(i)->SetGrid(1, 1); | |
1045 | canvasQAv0plusTOF->GetPad(i)->SetGrid(1, 1); | |
1046 | ||
1047 | canvasQAtpc->GetPad(i)->SetLogz(); | |
1048 | canvasQAtof->GetPad(i)->SetLogz(); | |
1049 | canvasQAv0->GetPad(i)->SetLogz(); | |
1050 | canvasQAv0plusTOF->GetPad(i)->SetLogz(); | |
1051 | ||
1052 | canvasQAtpc->GetPad(i)->SetLogx(); | |
1053 | canvasQAtof->GetPad(i)->SetLogx(); | |
1054 | canvasQAv0->GetPad(i)->SetLogx(); | |
1055 | canvasQAv0plusTOF->GetPad(i)->SetLogx(); | |
1056 | } | |
1057 | // | |
1058 | // 1. select and fit electrons | |
1059 | // | |
1060 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s | |
1061 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons | |
1062 | hist->GetAxis(4)->SetRangeUser(nSigmaElectronMinus,nSigmaElectronPlus); // 3sigma in TPC | |
1063 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1064 | TH2D * histElectron = hist->Projection(1,0); | |
1065 | histElectron->SetName("histElectron"); | |
1066 | histElectron->RebinX(4); | |
1067 | histElectron->GetXaxis()->SetRangeUser(0.2,3.0); | |
1068 | TObjArray arr; | |
1069 | FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1070 | TH1D * electronPoints = (TH1D *) arr.At(1)->Clone(); | |
1071 | // | |
1072 | TGraphErrors * electronGraph = new TGraphErrors(electronPoints); | |
1073 | electronGraph->SetName("electronGraph"); | |
1074 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { | |
1075 | Bool_t removePoint = electronGraph->GetY()[ip] < 10 || | |
1076 | electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.05 || | |
1077 | electronGraph->GetX()[ip] > 2.0 || | |
1078 | electronGraph->GetX()[ip] < 0.5; | |
1079 | if (removePoint) { | |
1080 | electronGraph->RemovePoint(ip); | |
1081 | ip--; | |
1082 | continue; | |
1083 | } | |
1084 | } | |
1085 | // | |
1086 | canvasQAtof->cd(1); | |
1087 | histElectron->SetTitle("electrons"); | |
1088 | histElectron->GetYaxis()->SetRangeUser(50, 120); | |
1089 | histElectron->GetXaxis()->SetTitle(momTitle.Data()); | |
1090 | histElectron->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1091 | histElectron->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1092 | histElectron->GetXaxis()->SetNoExponent(kTRUE); | |
1093 | histElectron->Draw("colz"); | |
1094 | electronPoints->SetMarkerStyle(24); | |
1095 | electronPoints->Draw("same"); | |
1096 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
1097 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
1098 | // | |
1099 | // 2. protons | |
1100 | // | |
1101 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s | |
1102 | hist->GetAxis(3)->SetRangeUser(3,3); // protons | |
1103 | // | |
1104 | //hist->GetAxis(4)->SetRangeUser(nSigmaProtonMinus,nSigmaProtonPlus); // protons --> not reliable, use PATTERN RECOGNITION | |
1105 | TH2D * histProtonTPC = hist->Projection(1,0); | |
1106 | histProtonTPC->SetName("histProtonTPC"); | |
1107 | histProtonTPC->GetXaxis()->SetRangeUser(0.15,0.7); | |
1108 | histProtonTPC->GetYaxis()->SetRangeUser(50, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins())); | |
1109 | ||
1110 | // PATTERN RECOGNITION | |
1111 | //OLD TF1 betaSq("betaSq","50./TMath::Power(x,1.3)",0.1,10); | |
1112 | // Add some additional Erf - cuts away kaons for pPb and PbPb and does not hurt in pp | |
1113 | TF1 betaSq("betaSq","50./TMath::Power(x,1.3) + (1-TMath::Erf((x-0.2) / 0.075)) * 100",0.1,10); | |
1114 | for(Int_t ix = 1; ix <= histProtonTPC->GetXaxis()->GetNbins(); ix++) { | |
1115 | for(Int_t jy = 1; jy <= histProtonTPC->GetYaxis()->GetNbins(); jy++) { | |
1116 | Float_t yPos = histProtonTPC->GetYaxis()->GetBinCenter(jy); | |
1117 | Float_t xPos = histProtonTPC->GetXaxis()->GetBinCenter(ix); | |
1118 | Int_t bin = histProtonTPC->GetBin(ix,jy); | |
1119 | if (yPos < betaSq.Eval(xPos)) histProtonTPC->SetBinContent(bin,0); | |
1120 | } | |
1121 | } | |
1122 | // | |
1123 | //TODO Use likelihood option -> Inside fitting range ~no contamination, but little statistics at low momenta requires L | |
1124 | FitSlicesY(histProtonTPC, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); | |
1125 | TH1D * protonPointsTPC = (TH1D *) arr.At(1)->Clone(); | |
1126 | // | |
1127 | TGraphErrors * protonTpcGraph = new TGraphErrors(protonPointsTPC); | |
1128 | protonTpcGraph->SetName("protonTpcGraph"); | |
1129 | for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) { | |
1130 | Bool_t removePoint = protonTpcGraph->GetY()[ip] < 10 || protonTpcGraph->GetEY()[ip]/protonTpcGraph->GetY()[ip] > 0.05 // TODO Larger tolerance, since this is only for the correction function, where 5% are still better than no data point | |
1131 | || protonTpcGraph->GetX()[ip] > (useV0s ? (isPPb ? 0.499 : 0.349) : 0.649); // If V0s are to be used, don't use TPC only protons to too high momenta; for pPb low statistics for the moment, therefore, go a bit further with TPC protons | |
1132 | //TODO NOW -> Changed range of TPC-signal -> Hopefully, sufficient statistics now for very low momenta!|| protonTpcGraph->GetX()[ip] < 0.19; // Added this cut because almost always bad statistics below 0.19 | |
1133 | if (removePoint) { | |
1134 | protonTpcGraph->RemovePoint(ip); | |
1135 | ip--; | |
1136 | continue; | |
1137 | } | |
1138 | } | |
1139 | // | |
1140 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1141 | TH2D * histProtonTOF = hist->Projection(1,0); | |
1142 | histProtonTOF->SetName("histProtonTOF"); | |
1143 | histProtonTOF->GetXaxis()->SetRangeUser(0.6,2.5); | |
1144 | ||
1145 | // Same pattern recognition as before | |
1146 | for(Int_t ix = 1; ix <= histProtonTOF->GetXaxis()->GetNbins(); ix++) { | |
1147 | for(Int_t jy = 1; jy <= histProtonTOF->GetYaxis()->GetNbins(); jy++) { | |
1148 | Float_t yPos = histProtonTOF->GetYaxis()->GetBinCenter(jy); | |
1149 | Float_t xPos = histProtonTOF->GetXaxis()->GetBinCenter(ix); | |
1150 | Int_t bin = histProtonTOF->GetBin(ix,jy); | |
1151 | if (yPos < betaSq.Eval(xPos)) histProtonTOF->SetBinContent(bin,0); | |
1152 | } | |
1153 | } | |
1154 | ||
1155 | FitSlicesY(histProtonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1156 | ||
1157 | TH1D * protonPointsTOF = (TH1D *)arr.At(1)->Clone(); | |
1158 | // | |
1159 | TGraphErrors * protonTofGraph = new TGraphErrors(protonPointsTOF); | |
1160 | protonTofGraph->SetName("protonTofGraph"); | |
1161 | for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) { | |
1162 | // If V0s are to be used, TPC+TOF protons are only for reference and can be plotted to higher momenta | |
1163 | Bool_t removePoint = protonTofGraph->GetY()[ip] < 10 || protonTofGraph->GetEY()[ip]/protonTofGraph->GetY()[ip] > 0.02 || protonTofGraph->GetX()[ip] > (useV0s ? 3 : 2) || protonTofGraph->GetX()[ip] < 0.65; | |
1164 | if (removePoint) { | |
1165 | protonTofGraph->RemovePoint(ip); | |
1166 | ip--; | |
1167 | continue; | |
1168 | } | |
1169 | } | |
1170 | // | |
1171 | canvasQAtpc->cd(2); | |
1172 | histProtonTPC->SetTitle("protons"); | |
1173 | histProtonTPC->GetXaxis()->SetTitle(momTitle.Data()); | |
1174 | histProtonTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1175 | histProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1176 | histProtonTPC->GetXaxis()->SetNoExponent(kTRUE); | |
1177 | histProtonTPC->Draw("colz"); | |
1178 | betaSq.DrawCopy("same"); | |
1179 | protonPointsTPC->SetMarkerStyle(20); | |
1180 | protonPointsTOF->SetMarkerStyle(24); | |
1181 | protonPointsTOF->Draw("same"); | |
1182 | protonPointsTPC->Draw("same"); | |
1183 | // | |
1184 | protonTpcGraph->SetMarkerStyle(26); | |
1185 | protonTpcGraph->SetMarkerColor(kMagenta); | |
1186 | protonTpcGraph->DrawClone("p"); | |
1187 | protonTofGraph->SetMarkerStyle(25); | |
1188 | protonTofGraph->SetMarkerColor(kMagenta); | |
1189 | protonTofGraph->DrawClone("p"); | |
1190 | ||
1191 | canvasQAtof->cd(2); | |
1192 | histProtonTOF->SetTitle("protons"); | |
1193 | histProtonTOF->GetYaxis()->SetRangeUser(30, 250); | |
1194 | histProtonTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1195 | histProtonTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1196 | histProtonTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1197 | histProtonTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1198 | histProtonTOF->Draw("colz"); | |
1199 | betaSq.DrawCopy("same"); | |
1200 | protonPointsTOF->Draw("same"); | |
1201 | protonPointsTPC->Draw("same"); | |
1202 | // | |
1203 | protonTpcGraph->DrawClone("p"); | |
1204 | protonTofGraph->DrawClone("p"); | |
1205 | ||
1206 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
1207 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
1208 | // | |
1209 | // 3. pions | |
1210 | // | |
1211 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s | |
1212 | hist->GetAxis(3)->SetRangeUser(1,1); // pions | |
1213 | // | |
1214 | Double_t lowerPionThreshold = 0.3; | |
1215 | ||
1216 | hist->GetAxis(4)->SetRangeUser(nSigmaPionMinus,nSigmaPionPlus); // pions | |
1217 | TH2D * histPionTPC = hist->Projection(1,0); | |
1218 | histPionTPC->SetName("histPionTPC"); | |
1219 | histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7); | |
1220 | FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1221 | ||
1222 | // In case of really bad splines comment last 5 lines and use the following 6 lines: | |
1223 | //TH2D * histPionTPC = hist->Projection(1,0); | |
1224 | //histPionTPC->SetName("histPionTPC"); | |
1225 | //histPionTPC->GetYaxis()->SetRangeUser(30,100); | |
1226 | //histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7); | |
1227 | //FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1228 | //lowerPionThreshold = 0.15; | |
1229 | ||
1230 | ||
1231 | TH1D * pionPointsTPC = (TH1D *) arr.At(1)->Clone(); | |
1232 | // | |
1233 | TGraphErrors * pionTpcGraph = new TGraphErrors(pionPointsTPC); | |
1234 | pionTpcGraph->SetName("pionTpcGraph"); | |
1235 | for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) { | |
1236 | Bool_t removePoint = pionTpcGraph->GetY()[ip] < 10 || pionTpcGraph->GetEY()[ip]/pionTpcGraph->GetY()[ip] > 0.02 || pionTpcGraph->GetX()[ip] > 0.5 | |
1237 | || pionTpcGraph->GetX()[ip] < lowerPionThreshold; // Exclude contamination by electrons | |
1238 | if (removePoint) { | |
1239 | pionTpcGraph->RemovePoint(ip); | |
1240 | ip--; | |
1241 | continue; | |
1242 | } | |
1243 | } | |
1244 | // | |
1245 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1246 | TH2D * histPionTOF = hist->Projection(1,0); | |
1247 | histPionTOF->SetName("histPionTOF"); | |
1248 | histPionTOF->GetXaxis()->SetRangeUser(0.5,1.1); | |
1249 | FitSlicesY(histPionTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1250 | TH1D * pionPointsTOF = (TH1D *) arr.At(1)->Clone(); | |
1251 | TGraphErrors * pionTofGraph = new TGraphErrors(pionPointsTOF); | |
1252 | pionTofGraph->SetName("pionTofGraph"); | |
1253 | for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) { | |
1254 | Bool_t removePoint = pionTofGraph->GetY()[ip] < 10 || pionTofGraph->GetEY()[ip]/pionTofGraph->GetY()[ip] > 0.02 || pionTofGraph->GetX()[ip] > 1.1 || pionTofGraph->GetX()[ip] < 0.5; // Changed by Ben (was 0.35) to avoid problems with TOF efficiency/systematics | |
1255 | if (removePoint) { | |
1256 | pionTofGraph->RemovePoint(ip); | |
1257 | ip--; | |
1258 | continue; | |
1259 | } | |
1260 | } | |
1261 | canvasQAtpc->cd(3); | |
1262 | histPionTPC->GetYaxis()->SetRangeUser(30, 90); | |
1263 | histPionTPC->SetTitle("pions"); | |
1264 | histPionTPC->GetXaxis()->SetTitle(momTitle.Data()); | |
1265 | histPionTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1266 | histPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1267 | histPionTPC->GetXaxis()->SetNoExponent(kTRUE); | |
1268 | histPionTPC->Draw("colz"); | |
1269 | pionPointsTPC->SetMarkerStyle(20); | |
1270 | pionPointsTOF->SetMarkerStyle(24); | |
1271 | pionPointsTOF->Draw("same"); | |
1272 | pionPointsTPC->Draw("same"); | |
1273 | // | |
1274 | pionTpcGraph->SetMarkerStyle(26); | |
1275 | pionTpcGraph->SetMarkerColor(kMagenta); | |
1276 | pionTpcGraph->DrawClone("p"); | |
1277 | pionTofGraph->SetMarkerStyle(25); | |
1278 | pionTofGraph->SetMarkerColor(kMagenta); | |
1279 | pionTofGraph->DrawClone("p"); | |
1280 | ||
1281 | canvasQAtof->cd(3); | |
1282 | histPionTOF->GetYaxis()->SetRangeUser(30, 80); | |
1283 | histPionTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1284 | histPionTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1285 | histPionTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1286 | histPionTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1287 | histPionTOF->Draw("colz"); | |
1288 | histPionTOF->SetTitle("pions"); | |
1289 | pionPointsTOF->Draw("same"); | |
1290 | pionPointsTPC->Draw("same"); | |
1291 | // | |
1292 | pionTpcGraph->DrawClone("p"); | |
1293 | pionTofGraph->DrawClone("p"); | |
1294 | ||
1295 | ||
1296 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
1297 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
1298 | // | |
1299 | // 4. kaons | |
1300 | // | |
1301 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s | |
1302 | hist->GetAxis(3)->SetRangeUser(2,2); // kaons | |
1303 | // | |
1304 | hist->GetAxis(4)->SetRangeUser(nSigmaKaonMinus,nSigmaKaonPlus); // kaons | |
1305 | TH2D * histKaonTPC = hist->Projection(1,0); | |
1306 | histKaonTPC->SetName("histKaonTPC"); | |
1307 | histKaonTPC->GetXaxis()->SetRangeUser(0.12,0.6999); | |
1308 | FitSlicesY(histKaonTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1309 | TH1D * kaonPointsTPC = (TH1D*) arr.At(1)->Clone(); | |
1310 | // | |
1311 | TGraphErrors * kaonTpcGraph = new TGraphErrors(kaonPointsTPC); | |
1312 | kaonTpcGraph->SetName("kaonTpcGraph"); | |
1313 | for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) { | |
1314 | Bool_t removePoint = kaonTpcGraph->GetY()[ip] < 10 || kaonTpcGraph->GetEY()[ip]/kaonTpcGraph->GetY()[ip] > 0.02 | |
1315 | || kaonTpcGraph->GetX()[ip] < 0.12 || kaonTpcGraph->GetX()[ip] > 0.319; //TODO BEN was > 0.399 | |
1316 | if (removePoint) { | |
1317 | kaonTpcGraph->RemovePoint(ip); | |
1318 | ip--; | |
1319 | continue; | |
1320 | } | |
1321 | } | |
1322 | // | |
1323 | ||
1324 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // sigma in TOF | |
1325 | ||
1326 | /* | |
1327 | hist->GetAxis(5)->SetRangeUser(-1.999,1.999); // sigma in TOF | |
1328 | ||
1329 | if (hist->GetAxis(3)->GetNbins() >= 5) { | |
1330 | printf("Using restricted eta range for TPC+TOF Kaons!\n"); | |
1331 | hist->GetAxis(3)->SetRangeUser(4,4); | |
1332 | } | |
1333 | else { | |
1334 | printf("WARNING: Data points for restricted eta range for TPC+TOF Kaons not available! Using full eta range (worse w.r.t. contamination)\n"); | |
1335 | }*/ | |
1336 | ||
1337 | TH2D * histKaonTOF = hist->Projection(1,0); | |
1338 | histKaonTOF->SetName("histKaonTOF"); | |
1339 | histKaonTOF->GetXaxis()->SetRangeUser(0.3,1.5); | |
1340 | FitSlicesY(histKaonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1341 | TH1D * kaonPointsTOF = (TH1D*) arr.At(1)->Clone(); | |
1342 | // | |
1343 | TGraphErrors * kaonTofGraph = new TGraphErrors(kaonPointsTOF); | |
1344 | kaonTofGraph->SetName("kaonTofGraph"); | |
1345 | // | |
1346 | for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) { | |
1347 | Bool_t removePoint = kaonTofGraph->GetY()[ip] < 10 || kaonTofGraph->GetEY()[ip]/kaonTofGraph->GetY()[ip] > 0.02 || kaonTofGraph->GetX()[ip] > 1.0 | |
1348 | || kaonTofGraph->GetX()[ip] < 0.36; //TODO BEN NOW was 0.5 | |
1349 | if (removePoint) { | |
1350 | kaonTofGraph->RemovePoint(ip); | |
1351 | ip--; | |
1352 | continue; | |
1353 | } | |
1354 | } | |
1355 | // | |
1356 | canvasQAtpc->cd(4); | |
1357 | histKaonTPC->SetTitle("kaons"); | |
1358 | histKaonTPC->GetYaxis()->SetRangeUser(30, 500); | |
1359 | histKaonTPC->GetXaxis()->SetTitle(momTitle.Data()); | |
1360 | histKaonTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1361 | histKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1362 | histKaonTPC->GetXaxis()->SetNoExponent(kTRUE); | |
1363 | histKaonTPC->Draw("colz"); | |
1364 | kaonPointsTPC->SetMarkerStyle(20); | |
1365 | kaonPointsTOF->SetMarkerStyle(24); | |
1366 | kaonPointsTOF->Draw("same"); | |
1367 | kaonPointsTPC->Draw("same"); | |
1368 | // | |
1369 | kaonTpcGraph->SetMarkerStyle(26); | |
1370 | kaonTpcGraph->SetMarkerColor(kMagenta); | |
1371 | kaonTpcGraph->DrawClone("p"); | |
1372 | kaonTofGraph->SetMarkerStyle(25); | |
1373 | kaonTofGraph->SetMarkerColor(kMagenta); | |
1374 | kaonTofGraph->DrawClone("p"); | |
1375 | ||
1376 | ||
1377 | canvasQAtof->cd(4); | |
1378 | histKaonTOF->GetYaxis()->SetRangeUser(30, 250); | |
1379 | histKaonTOF->SetTitle("kaons"); | |
1380 | histKaonTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1381 | histKaonTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1382 | histKaonTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1383 | histKaonTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1384 | histKaonTOF->Draw("colz"); | |
1385 | kaonPointsTOF->Draw("same"); | |
1386 | kaonPointsTPC->Draw("same"); | |
1387 | kaonTpcGraph->DrawClone("p"); | |
1388 | kaonTofGraph->DrawClone("p"); | |
1389 | ||
1390 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES | |
1391 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
1392 | ||
1393 | ||
1394 | // | |
1395 | // 5. V0 electrons | |
1396 | // | |
1397 | hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons | |
1398 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons | |
1399 | ||
1400 | // | |
1401 | TH2D * histElectronV0 = hist->Projection(1,0); | |
1402 | histElectronV0->SetName("histElectronV0"); | |
1403 | //histElectronV0->RebinX(2); | |
1404 | histElectronV0->GetXaxis()->SetRangeUser(0.15,5.0); | |
1405 | FitSlicesY(histElectronV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1406 | TH1D * electronPointsV0 = (TH1D*) arr.At(1)->Clone(); | |
1407 | // | |
1408 | TGraphErrors * electronV0Graph = new TGraphErrors(electronPointsV0); | |
1409 | electronV0Graph->SetName("electronV0Graph"); | |
1410 | for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) { | |
1411 | Bool_t removePoint = electronV0Graph->GetY()[ip] < 10 || electronV0Graph->GetEY()[ip]/electronV0Graph->GetY()[ip] > 0.02; | |
1412 | ||
1413 | if (removePoint) { | |
1414 | electronV0Graph->RemovePoint(ip); | |
1415 | ip--; | |
1416 | continue; | |
1417 | } | |
1418 | } | |
1419 | // | |
1420 | canvasQAv0->cd(1); | |
1421 | histElectronV0->SetTitle("V0 electrons"); | |
1422 | histElectronV0->GetYaxis()->SetRangeUser(50, 120); | |
1423 | histElectronV0->GetXaxis()->SetTitle(momTitle.Data()); | |
1424 | histElectronV0->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1425 | histElectronV0->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1426 | histElectronV0->GetXaxis()->SetNoExponent(kTRUE); | |
1427 | histElectronV0->SetStats(0); | |
1428 | histElectronV0->Draw("colz"); | |
1429 | electronPointsV0->SetMarkerStyle(34); | |
1430 | electronPointsV0->Draw("same"); | |
1431 | electronGraph->SetMarkerStyle(25); | |
1432 | electronGraph->SetMarkerColor(kMagenta); | |
1433 | electronGraph->DrawClone("p"); | |
1434 | ||
1435 | canvasQAv0DeDxPurityEl->cd(); | |
1436 | gStyle->SetOptTitle(0); | |
1437 | TH1D* hElNew = (TH1D*)histElectronV0->DrawClone("colz"); | |
1438 | hElNew->GetXaxis()->SetTitleOffset(1.0); | |
1439 | hElNew->SetTitle(""); | |
1440 | electronPointsV0->DrawClone("same"); | |
1441 | gStyle->SetOptTitle(1); | |
1442 | ||
1443 | ||
1444 | ||
1445 | canvasQAtof->cd(1); | |
1446 | electronV0Graph->SetMarkerStyle(24); | |
1447 | electronV0Graph->SetMarkerColor(kMagenta); | |
1448 | electronV0Graph->DrawClone("p"); | |
1449 | ||
1450 | // | |
1451 | // 6. V0 pions | |
1452 | // | |
1453 | hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions | |
1454 | hist->GetAxis(3)->SetRangeUser(1,1); // pions | |
1455 | // | |
1456 | TH2D * histPionV0 = hist->Projection(1,0); | |
1457 | histPionV0->SetName("histPionV0"); | |
1458 | histPionV0->GetXaxis()->SetRangeUser(0.15,10.); | |
1459 | FitSlicesY(histPionV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1460 | TH1D * pionPointsV0 = (TH1D*) arr.At(1)->Clone(); | |
1461 | // | |
1462 | TGraphErrors * pionV0Graph = new TGraphErrors(pionPointsV0); | |
1463 | pionV0Graph->SetName("pionV0Graph"); | |
1464 | for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) { | |
1465 | Bool_t removePoint = pionV0Graph->GetY()[ip] < 10 || pionV0Graph->GetEY()[ip]/pionV0Graph->GetY()[ip] > 0.02; | |
1466 | ||
1467 | if (removePoint) { | |
1468 | pionV0Graph->RemovePoint(ip); | |
1469 | ip--; | |
1470 | continue; | |
1471 | } | |
1472 | } | |
1473 | // | |
1474 | canvasQAv0->cd(3); | |
1475 | histPionV0->SetTitle("V0 pions"); | |
1476 | histPionV0->GetYaxis()->SetRangeUser(30, 100); | |
1477 | histPionV0->GetXaxis()->SetTitle(momTitle.Data()); | |
1478 | histPionV0->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1479 | histPionV0->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1480 | histPionV0->GetXaxis()->SetNoExponent(kTRUE); | |
1481 | histPionV0->Draw("colz"); | |
1482 | pionPointsV0->SetMarkerStyle(34); | |
1483 | pionPointsV0->Draw("same"); | |
1484 | pionTofGraph->DrawClone("p"); | |
1485 | pionTpcGraph->DrawClone("p"); | |
1486 | ||
1487 | canvasQAv0DeDxPurityPi->cd(); | |
1488 | gStyle->SetOptTitle(0); | |
1489 | TH1D* hPiNew = (TH1D*)histPionV0->DrawClone("colz"); | |
1490 | hPiNew->GetXaxis()->SetTitleOffset(1.0); | |
1491 | hPiNew->SetTitle(""); | |
1492 | pionPointsV0->DrawClone("same"); | |
1493 | gStyle->SetOptTitle(1); | |
1494 | ||
1495 | canvasQAtof->cd(3); | |
1496 | pionV0Graph->SetMarkerStyle(24); | |
1497 | pionV0Graph->SetMarkerColor(kMagenta); | |
1498 | pionV0Graph->DrawClone("p"); | |
1499 | ||
1500 | canvasQAtpc->cd(3); | |
1501 | pionV0Graph->DrawClone("p"); | |
1502 | ||
1503 | // | |
1504 | // 6. V0 protons | |
1505 | // | |
1506 | hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons | |
1507 | hist->GetAxis(3)->SetRangeUser(3,3); // protons | |
1508 | // | |
1509 | TH2D * histProtonV0 = hist->Projection(1,0); | |
1510 | histProtonV0->SetName("histProtonV0"); | |
1511 | histProtonV0->GetXaxis()->SetRangeUser(0.15,10.); | |
1512 | ||
1513 | // Remove contamination due to el and pions and low momenta because there the statistics | |
1514 | // for the V0 protons goes down and becomes comparable with the contamination | |
1515 | // -> This spoils the fit, if the contamination is not removed | |
1516 | ||
1517 | // PATTERN RECOGNITION | |
1518 | Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.6); | |
1519 | ||
1520 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { | |
1521 | for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) { | |
1522 | Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy); | |
1523 | Float_t xPos = histProtonV0->GetXaxis()->GetBinCenter(ix); | |
1524 | Int_t bin = histProtonV0->GetBin(ix,jy); | |
1525 | if (yPos < betaSq.Eval(xPos)) histProtonV0->SetBinContent(bin,0); | |
1526 | } | |
1527 | } | |
1528 | ||
1529 | /* | |
1530 | Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.4); | |
1531 | Double_t threshold = 120.; | |
1532 | ||
1533 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { | |
1534 | for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) { | |
1535 | Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy); | |
1536 | Int_t bin = histProtonV0->GetBin(ix,jy); | |
1537 | if (yPos < threshold) histProtonV0->SetBinContent(bin,0); | |
1538 | else break; | |
1539 | } | |
1540 | } | |
1541 | */ | |
1542 | ||
1543 | FitSlicesY(histProtonV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1544 | TH1D * protonPointsV0 = (TH1D*) arr.At(1)->Clone(); | |
1545 | // | |
1546 | TGraphErrors * protonV0Graph = new TGraphErrors(protonPointsV0); | |
1547 | protonV0Graph->SetName("protonV0Graph"); | |
1548 | for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) { | |
1549 | Bool_t removePoint = protonV0Graph->GetY()[ip] < 10 || protonV0Graph->GetEY()[ip]/protonV0Graph->GetY()[ip] > 0.02; | |
1550 | ||
1551 | if (removePoint) { | |
1552 | protonV0Graph->RemovePoint(ip); | |
1553 | ip--; | |
1554 | continue; | |
1555 | } | |
1556 | } | |
1557 | // | |
1558 | canvasQAv0->cd(2); | |
1559 | histProtonV0->SetTitle("V0 protons"); | |
1560 | histProtonV0->GetXaxis()->SetTitle(momTitle.Data()); | |
1561 | histProtonV0->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1562 | histProtonV0->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1563 | histProtonV0->GetXaxis()->SetNoExponent(kTRUE); | |
1564 | histProtonV0->Draw("colz"); | |
1565 | protonPointsV0->SetMarkerStyle(34); | |
1566 | protonPointsV0->Draw("same"); | |
1567 | protonTofGraph->DrawClone("p"); | |
1568 | protonTpcGraph->DrawClone("p"); | |
1569 | ||
1570 | canvasQAv0DeDxPurityPr->cd(); | |
1571 | gStyle->SetOptTitle(0); | |
1572 | TH1D* hPrNew = (TH1D*)histProtonV0->DrawClone("colz"); | |
1573 | hPrNew->GetXaxis()->SetTitleOffset(1.0); | |
1574 | hPrNew->SetTitle(""); | |
1575 | protonPointsV0->DrawClone("same"); | |
1576 | gStyle->SetOptTitle(1); | |
1577 | ||
1578 | canvasQAtof->cd(2); | |
1579 | protonV0Graph->SetMarkerStyle(24); | |
1580 | protonV0Graph->SetMarkerColor(kMagenta); | |
1581 | protonV0Graph->DrawClone("p"); | |
1582 | ||
1583 | canvasQAtpc->cd(2); | |
1584 | protonV0Graph->DrawClone("p"); | |
1585 | ||
1586 | ||
1587 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES | |
1588 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES | |
1589 | ||
1590 | ||
1591 | // | |
1592 | // 5. V0 electrons + TOF | |
1593 | // | |
1594 | hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons | |
1595 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons | |
1596 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1597 | // | |
1598 | TH2D * histElectronV0plusTOF = hist->Projection(1,0); | |
1599 | histElectronV0plusTOF->SetName("histElectronV0plusTOF"); | |
1600 | histElectronV0plusTOF->RebinX(2); | |
1601 | histElectronV0plusTOF->GetXaxis()->SetRangeUser(0.2,5.0); | |
1602 | FitSlicesY(histElectronV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1603 | TH1D * electronPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); | |
1604 | // | |
1605 | TGraphErrors * electronV0plusTOFGraph = new TGraphErrors(electronPointsV0plusTOF); | |
1606 | electronV0plusTOFGraph->SetName("electronV0plusTOFGraph"); | |
1607 | for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) { | |
1608 | Bool_t removePoint = electronV0plusTOFGraph->GetY()[ip] < 10 || electronV0plusTOFGraph->GetEY()[ip]/electronV0plusTOFGraph->GetY()[ip] > 0.02; | |
1609 | ||
1610 | if (removePoint) { | |
1611 | electronV0plusTOFGraph->RemovePoint(ip); | |
1612 | ip--; | |
1613 | continue; | |
1614 | } | |
1615 | } | |
1616 | // | |
1617 | canvasQAv0plusTOF->cd(1); | |
1618 | histElectronV0plusTOF->SetTitle("V0+TOF electrons"); | |
1619 | histElectronV0plusTOF->GetYaxis()->SetRangeUser(50, 120); | |
1620 | histElectronV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1621 | histElectronV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1622 | histElectronV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1623 | histElectronV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1624 | histElectronV0plusTOF->Draw("colz"); | |
1625 | electronPointsV0plusTOF->SetMarkerStyle(29); | |
1626 | electronPointsV0plusTOF->Draw("same"); | |
1627 | electronGraph->DrawClone("p"); | |
1628 | electronV0Graph->DrawClone("p"); | |
1629 | ||
1630 | canvasQAv0->cd(1); | |
1631 | electronV0plusTOFGraph->SetMarkerStyle(30); | |
1632 | electronV0plusTOFGraph->SetMarkerColor(kMagenta); | |
1633 | electronV0plusTOFGraph->DrawClone("p"); | |
1634 | ||
1635 | canvasQAtof->cd(1); | |
1636 | electronV0plusTOFGraph->DrawClone("p"); | |
1637 | ||
1638 | // | |
1639 | // 6. V0 pions + TOF | |
1640 | // | |
1641 | hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions | |
1642 | hist->GetAxis(3)->SetRangeUser(1,1); // pions | |
1643 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1644 | // | |
1645 | TH2D * histPionV0plusTOF = hist->Projection(1,0); | |
1646 | histPionV0plusTOF->SetName("histPionV0plusTOF"); | |
1647 | histPionV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.); | |
1648 | FitSlicesY(histPionV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1649 | TH1D * pionPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); | |
1650 | // | |
1651 | TGraphErrors * pionV0plusTOFGraph = new TGraphErrors(pionPointsV0plusTOF); | |
1652 | pionV0plusTOFGraph->SetName("pionV0plusTOFGraph"); | |
1653 | for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) { | |
1654 | Bool_t removePoint = pionV0plusTOFGraph->GetY()[ip] < 10 || pionV0plusTOFGraph->GetEY()[ip]/pionV0plusTOFGraph->GetY()[ip] > 0.02; | |
1655 | ||
1656 | if (removePoint) { | |
1657 | pionV0plusTOFGraph->RemovePoint(ip); | |
1658 | ip--; | |
1659 | continue; | |
1660 | } | |
1661 | } | |
1662 | // | |
1663 | canvasQAv0plusTOF->cd(3); | |
1664 | histPionV0plusTOF->SetTitle("V0+TOF pions"); | |
1665 | histPionV0plusTOF->GetYaxis()->SetRangeUser(30, 100); | |
1666 | histPionV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1667 | histPionV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1668 | histPionV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1669 | histPionV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1670 | histPionV0plusTOF->Draw("colz"); | |
1671 | pionPointsV0plusTOF->SetMarkerStyle(29); | |
1672 | pionPointsV0plusTOF->Draw("same"); | |
1673 | pionV0Graph->DrawClone("p"); | |
1674 | pionTofGraph->DrawClone("p"); | |
1675 | pionTpcGraph->DrawClone("p"); | |
1676 | ||
1677 | canvasQAv0->cd(3); | |
1678 | pionV0plusTOFGraph->SetMarkerStyle(30); | |
1679 | pionV0plusTOFGraph->SetMarkerColor(kMagenta); | |
1680 | pionV0plusTOFGraph->DrawClone("p"); | |
1681 | ||
1682 | canvasQAtof->cd(3); | |
1683 | pionV0plusTOFGraph->DrawClone("p"); | |
1684 | ||
1685 | canvasQAtpc->cd(3); | |
1686 | pionV0plusTOFGraph->DrawClone("p"); | |
1687 | ||
1688 | // | |
1689 | // 6. V0 protons + TOF | |
1690 | // | |
1691 | hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons | |
1692 | hist->GetAxis(3)->SetRangeUser(3,3); // protons | |
1693 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF | |
1694 | // | |
1695 | TH2D * histProtonV0plusTOF = hist->Projection(1,0); | |
1696 | histProtonV0plusTOF->SetName("histProtonV0plusTOF"); | |
1697 | histProtonV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.); | |
1698 | ||
1699 | // Remove contamination due to el and pions and low momenta because there the statistics | |
1700 | // for the V0 protons goes down and becomes comparable with the contamination | |
1701 | // -> This spoils the fit, if the contamination is not removed | |
1702 | ||
1703 | // PATTERN RECOGNITION | |
1704 | correctionUpToBin = histProtonV0plusTOF->GetXaxis()->FindBin(0.6); | |
1705 | ||
1706 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { | |
1707 | for(Int_t jy = 1; jy <= histProtonV0plusTOF->GetYaxis()->GetNbins(); jy++) { | |
1708 | Float_t yPos = histProtonV0plusTOF->GetYaxis()->GetBinCenter(jy); | |
1709 | Float_t xPos = histProtonV0plusTOF->GetXaxis()->GetBinCenter(ix); | |
1710 | Int_t bin = histProtonV0plusTOF->GetBin(ix,jy); | |
1711 | if (yPos < betaSq.Eval(xPos)) histProtonV0plusTOF->SetBinContent(bin,0); | |
1712 | } | |
1713 | } | |
1714 | ||
1715 | FitSlicesY(histProtonV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); | |
1716 | TH1D * protonPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); | |
1717 | // | |
1718 | TGraphErrors * protonV0plusTOFGraph = new TGraphErrors(protonPointsV0plusTOF); | |
1719 | protonV0plusTOFGraph->SetName("protonV0plusTOFGraph"); | |
1720 | for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) { | |
1721 | Bool_t removePoint = protonV0plusTOFGraph->GetY()[ip] < 10 || protonV0plusTOFGraph->GetEY()[ip]/protonV0plusTOFGraph->GetY()[ip] > 0.02; | |
1722 | ||
1723 | if (removePoint) { | |
1724 | protonV0plusTOFGraph->RemovePoint(ip); | |
1725 | ip--; | |
1726 | continue; | |
1727 | } | |
1728 | } | |
1729 | // | |
1730 | canvasQAv0plusTOF->cd(2); | |
1731 | histProtonV0plusTOF->SetTitle("V0+TOF protons"); | |
1732 | histProtonV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); | |
1733 | histProtonV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
1734 | histProtonV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1735 | histProtonV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); | |
1736 | histProtonV0plusTOF->Draw("colz"); | |
1737 | protonPointsV0plusTOF->SetMarkerStyle(29); | |
1738 | protonPointsV0plusTOF->Draw("same"); | |
1739 | protonV0Graph->DrawClone("p"); | |
1740 | protonTofGraph->DrawClone("p"); | |
1741 | protonTpcGraph->DrawClone("p"); | |
1742 | ||
1743 | canvasQAv0->cd(2); | |
1744 | protonV0plusTOFGraph->SetMarkerStyle(30); | |
1745 | protonV0plusTOFGraph->SetMarkerColor(kMagenta); | |
1746 | protonV0plusTOFGraph->DrawClone("p"); | |
1747 | ||
1748 | canvasQAtof->cd(2); | |
1749 | protonV0plusTOFGraph->DrawClone("p"); | |
1750 | ||
1751 | canvasQAtpc->cd(2); | |
1752 | protonV0plusTOFGraph->DrawClone("p"); | |
1753 | ||
1754 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES | |
1755 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES | |
1756 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES | |
1757 | ||
1758 | ||
1759 | // Clone (V0, V0+TOF) graphs: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that | |
1760 | // these graphs can be used for the BB fit. The original graph will keep all point in order to determine the correction and response | |
1761 | // functions | |
1762 | TGraphErrors * electronV0GraphForBBfit = new TGraphErrors(*electronV0Graph); | |
1763 | electronV0GraphForBBfit->SetName("electronV0GraphForBBfit"); | |
1764 | ||
1765 | for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) { | |
1766 | Bool_t removePoint = electronV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.4);// || electronV0GraphForBBfit->GetX()[ip] > 2.2; | |
1767 | ||
1768 | if (removePoint) { | |
1769 | electronV0GraphForBBfit->RemovePoint(ip); | |
1770 | ip--; | |
1771 | continue; | |
1772 | } | |
1773 | } | |
1774 | ||
1775 | TGraphErrors * pionV0GraphForBBfit = new TGraphErrors(*pionV0Graph); | |
1776 | pionV0GraphForBBfit->SetName("pionV0GraphForBBfit"); | |
1777 | ||
1778 | for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) { | |
1779 | Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.0 : 0.76);// || pionV0GraphForBBfit->GetX()[ip] > 6.; //TODO NOW was < 0.4 before | |
1780 | //Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < 0.4 || pionV0GraphForBBfit->GetX()[ip] > 0.399;// In this case NOT used for BB fit | |
1781 | if (removePoint) { | |
1782 | pionV0GraphForBBfit->RemovePoint(ip); | |
1783 | ip--; | |
1784 | continue; | |
1785 | } | |
1786 | } | |
1787 | ||
1788 | TGraphErrors * protonV0GraphForBBfit = new TGraphErrors(*protonV0Graph); | |
1789 | protonV0GraphForBBfit->SetName("protonV0GraphForBBfit"); | |
1790 | ||
1791 | for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) { | |
1792 | Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 5.0; //TODO NOW was 2.5 before | |
1793 | //Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 0.599;// In this case NOT used for BB fit | |
1794 | if (removePoint) { | |
1795 | protonV0GraphForBBfit->RemovePoint(ip); | |
1796 | ip--; | |
1797 | continue; | |
1798 | } | |
1799 | } | |
1800 | ||
1801 | ||
1802 | // V0 plus TOF | |
1803 | TGraphErrors * electronV0plusTOFGraphForBBfit = new TGraphErrors(*electronV0plusTOFGraph); | |
1804 | electronV0plusTOFGraphForBBfit->SetName("electronV0plusTOFGraphForBBfit"); | |
1805 | ||
1806 | for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1807 | Bool_t removePoint = electronV0plusTOFGraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.6);//0.4 || electronV0plusTOFGraphForBBfit->GetX()[ip] > 1.3799; | |
1808 | if (removePoint) { | |
1809 | electronV0plusTOFGraphForBBfit->RemovePoint(ip); | |
1810 | ip--; | |
1811 | continue; | |
1812 | } | |
1813 | } | |
1814 | ||
1815 | TGraphErrors * pionV0plusTOFGraphForBBfit = new TGraphErrors(*pionV0plusTOFGraph); | |
1816 | pionV0plusTOFGraphForBBfit->SetName("pionV0plusTOFGraphForBBfit"); | |
1817 | ||
1818 | for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1819 | Bool_t removePoint = pionV0plusTOFGraphForBBfit->GetX()[ip] < 0.4;// || pionV0plusTOFGraphForBBfit->GetX()[ip] > 6.; | |
1820 | if (removePoint) { | |
1821 | pionV0plusTOFGraphForBBfit->RemovePoint(ip); | |
1822 | ip--; | |
1823 | continue; | |
1824 | } | |
1825 | } | |
1826 | ||
1827 | TGraphErrors * protonV0plusTOFGraphForBBfit = new TGraphErrors(*protonV0plusTOFGraph); | |
1828 | protonV0plusTOFGraphForBBfit->SetName("protonV0plusTOFGraphForBBfit"); | |
1829 | ||
1830 | for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1831 | Bool_t removePoint = protonV0plusTOFGraphForBBfit->GetX()[ip] < 0.9 || protonV0plusTOFGraphForBBfit->GetX()[ip] > 2.5; | |
1832 | if (removePoint) { | |
1833 | protonV0plusTOFGraphForBBfit->RemovePoint(ip); | |
1834 | ip--; | |
1835 | continue; | |
1836 | } | |
1837 | } | |
1838 | ||
1839 | ||
1840 | // | |
1841 | // rescale graphs to betaGamma | |
1842 | // | |
1843 | kaonTofGraph->SetMarkerStyle(21); | |
1844 | kaonTpcGraph->SetMarkerStyle(22); | |
1845 | for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) { | |
1846 | kaonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
1847 | kaonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
1848 | } | |
1849 | for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) { | |
1850 | kaonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
1851 | kaonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
1852 | } | |
1853 | // | |
1854 | electronGraph->SetMarkerStyle(22); | |
1855 | electronV0Graph->SetMarkerStyle(34); | |
1856 | electronV0GraphForBBfit->SetMarkerStyle(34); | |
1857 | electronV0plusTOFGraph->SetMarkerStyle(30); | |
1858 | electronV0plusTOFGraphForBBfit->SetMarkerStyle(30); | |
1859 | electronGraph->SetMarkerColor(2); | |
1860 | electronV0Graph->SetMarkerColor(2); | |
1861 | electronV0GraphForBBfit->SetMarkerColor(2); | |
1862 | electronV0plusTOFGraph->SetMarkerColor(2); | |
1863 | electronV0plusTOFGraphForBBfit->SetMarkerColor(2); | |
1864 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { | |
1865 | electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1866 | electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1867 | } | |
1868 | for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) { | |
1869 | electronV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1870 | electronV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1871 | } | |
1872 | for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) { | |
1873 | electronV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1874 | electronV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1875 | } | |
1876 | for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) { | |
1877 | electronV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1878 | electronV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1879 | } | |
1880 | for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1881 | electronV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1882 | electronV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
1883 | } | |
1884 | // | |
1885 | pionTofGraph->SetMarkerStyle(21); | |
1886 | pionTpcGraph->SetMarkerStyle(22); | |
1887 | pionV0Graph->SetMarkerStyle(34); | |
1888 | pionV0GraphForBBfit->SetMarkerStyle(34); | |
1889 | pionV0plusTOFGraph->SetMarkerStyle(30); | |
1890 | pionV0plusTOFGraphForBBfit->SetMarkerStyle(30); | |
1891 | pionTofGraph->SetMarkerColor(3); | |
1892 | pionTpcGraph->SetMarkerColor(3); | |
1893 | pionV0Graph->SetMarkerColor(3); | |
1894 | pionV0GraphForBBfit->SetMarkerColor(3); | |
1895 | pionV0plusTOFGraph->SetMarkerColor(3); | |
1896 | pionV0plusTOFGraphForBBfit->SetMarkerColor(3); | |
1897 | for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) { | |
1898 | pionTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1899 | pionTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1900 | } | |
1901 | for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) { | |
1902 | pionTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1903 | pionTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1904 | } | |
1905 | for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) { | |
1906 | pionV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1907 | pionV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1908 | } | |
1909 | for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) { | |
1910 | pionV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1911 | pionV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1912 | } | |
1913 | for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) { | |
1914 | pionV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1915 | pionV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1916 | } | |
1917 | for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1918 | pionV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1919 | pionV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
1920 | } | |
1921 | // | |
1922 | protonTofGraph->SetMarkerStyle(21); | |
1923 | protonTpcGraph->SetMarkerStyle(22); | |
1924 | protonV0Graph->SetMarkerStyle(34); | |
1925 | protonV0GraphForBBfit->SetMarkerStyle(34); | |
1926 | protonV0plusTOFGraph->SetMarkerStyle(30); | |
1927 | protonV0plusTOFGraphForBBfit->SetMarkerStyle(30); | |
1928 | protonTofGraph->SetMarkerColor(4); | |
1929 | protonTpcGraph->SetMarkerColor(4); | |
1930 | protonV0Graph->SetMarkerColor(4); | |
1931 | protonV0GraphForBBfit->SetMarkerColor(4); | |
1932 | protonV0plusTOFGraph->SetMarkerColor(4); | |
1933 | protonV0plusTOFGraphForBBfit->SetMarkerColor(4); | |
1934 | for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) { | |
1935 | protonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1936 | protonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1937 | } | |
1938 | for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) { | |
1939 | protonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1940 | protonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1941 | } | |
1942 | for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) { | |
1943 | protonV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1944 | protonV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1945 | } | |
1946 | for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) { | |
1947 | protonV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1948 | protonV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1949 | } | |
1950 | for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) { | |
1951 | protonV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1952 | protonV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1953 | } | |
1954 | for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) { | |
1955 | protonV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1956 | protonV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
1957 | } | |
1958 | // | |
1959 | // | |
1960 | TGraphErrors * graphAll = new TGraphErrors(); | |
1961 | TList * listColl = new TList(); | |
1962 | if (!useV0s) { | |
1963 | //listColl->Add(protonTpcGraph); | |
1964 | //listColl->Add(kaonTpcGraph); | |
1965 | listColl->Add(pionTpcGraph); | |
1966 | listColl->Add(electronGraph); | |
1967 | listColl->Add(protonTofGraph); | |
1968 | listColl->Add(kaonTofGraph); | |
1969 | listColl->Add(pionTofGraph); | |
1970 | } | |
1971 | else { | |
1972 | listColl->Add(pionV0GraphForBBfit); | |
1973 | //listColl->Add(electronV0GraphForBBfit); | |
1974 | listColl->Add(protonV0GraphForBBfit); | |
1975 | ||
1976 | //listColl->Add(pionV0plusTOFGraphForBBfit); | |
1977 | listColl->Add(electronV0plusTOFGraphForBBfit); | |
1978 | //listColl->Add(protonV0plusTOFGraphForBBfit); | |
1979 | } | |
1980 | MergeGraphErrors(graphAll, listColl); | |
1981 | graphAll->SetNameTitle("beamDataPoints","beamDataPoints"); | |
1982 | ||
1983 | // Create graph out of TPC only, V0, V0+TOF and TPC+TOF to determine correction functions | |
1984 | TList * listCollPr = new TList(); | |
1985 | if (useV0s) { | |
1986 | listCollPr->Add(protonTpcGraph); | |
1987 | //listCollPr->Add(protonTofGraph); | |
1988 | ||
1989 | TGraphErrors * protonV0GraphCut = new TGraphErrors(*protonV0Graph); | |
1990 | protonV0GraphCut->SetName("protonV0GraphCut"); | |
1991 | for(Int_t ip = 0; ip < protonV0GraphCut->GetN(); ip ++) { | |
1992 | Double_t mom = protonV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton); | |
1993 | //Bool_t removePoint = mom >= 0.6 || mom < 0.35; // Will take TPC protons instead (statistics) for very low momenta | |
1994 | Bool_t removePoint = mom < (isPPb ? 0.6 : 0.35); // Will take TPC protons instead (statistics) for very low momenta; | |
1995 | // for pPb low statistics for the moment, therefore, go a bit further with TPC protons | |
1996 | ||
1997 | ||
1998 | if (removePoint) { | |
1999 | protonV0GraphCut->RemovePoint(ip); | |
2000 | ip--; | |
2001 | continue; | |
2002 | } | |
2003 | } | |
2004 | /* | |
2005 | TGraphErrors * protonV0plusTOFGraphCut = new TGraphErrors(*protonV0plusTOFGraph); | |
2006 | protonV0plusTOFGraphCut->SetName("protonV0plusTOFGraphCut"); | |
2007 | for(Int_t ip = 0; ip < protonV0plusTOFGraphCut->GetN(); ip ++) { | |
2008 | Double_t mom = protonV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton); | |
2009 | Bool_t removePoint = mom < 0.6; | |
2010 | ||
2011 | if (removePoint) { | |
2012 | protonV0plusTOFGraphCut->RemovePoint(ip); | |
2013 | ip--; | |
2014 | continue; | |
2015 | } | |
2016 | }*/ | |
2017 | ||
2018 | listCollPr->Add(protonV0GraphCut); | |
2019 | //listCollPr->Add(protonV0plusTOFGraphCut); | |
2020 | } | |
2021 | else { | |
2022 | listCollPr->Add(protonTpcGraph); | |
2023 | listCollPr->Add(protonTofGraph); | |
2024 | } | |
2025 | TGraphErrors * graphPrCombined = new TGraphErrors(); | |
2026 | MergeGraphErrors(graphPrCombined, listCollPr); | |
2027 | graphPrCombined->SetMarkerStyle(22); | |
2028 | graphPrCombined->SetMarkerColor(4); | |
2029 | graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined"); | |
2030 | ||
2031 | TList * listCollPi = new TList(); | |
2032 | if (useV0s) { | |
2033 | //listCollPi->Add(pionTpcGraph); | |
2034 | //listCollPi->Add(pionTofGraph); | |
2035 | ||
2036 | TGraphErrors * pionV0GraphCut = new TGraphErrors(*pionV0Graph); | |
2037 | pionV0GraphCut->SetName("pionV0GraphCut"); | |
2038 | for(Int_t ip = 0; ip < pionV0GraphCut->GetN(); ip ++) { | |
2039 | //Double_t mom = pionV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion); | |
2040 | Bool_t removePoint = kFALSE;//mom >= 0.4; | |
2041 | //|| mom < 0.2;//TODO NOW NEEDED? | |
2042 | if (removePoint) { | |
2043 | pionV0GraphCut->RemovePoint(ip); | |
2044 | ip--; | |
2045 | continue; | |
2046 | } | |
2047 | } | |
2048 | /* | |
2049 | TGraphErrors * pionV0plusTOFGraphCut = new TGraphErrors(*pionV0plusTOFGraph); | |
2050 | pionV0plusTOFGraphCut->SetName("pionV0plusTOFGraphCut"); | |
2051 | for(Int_t ip = 0; ip < pionV0plusTOFGraphCut->GetN(); ip ++) { | |
2052 | Double_t mom = pionV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion); | |
2053 | Bool_t removePoint = mom < 0.4; | |
2054 | ||
2055 | if (removePoint) { | |
2056 | pionV0plusTOFGraphCut->RemovePoint(ip); | |
2057 | ip--; | |
2058 | continue; | |
2059 | } | |
2060 | }*/ | |
2061 | ||
2062 | listCollPi->Add(pionV0GraphCut); | |
2063 | //listCollPi->Add(pionV0plusTOFGraphCut); | |
2064 | } | |
2065 | else { | |
2066 | listCollPi->Add(pionTpcGraph); | |
2067 | listCollPi->Add(pionTofGraph); | |
2068 | } | |
2069 | TGraphErrors * graphPiCombined = new TGraphErrors(); | |
2070 | MergeGraphErrors(graphPiCombined, listCollPi); | |
2071 | graphPiCombined->SetMarkerStyle(22); | |
2072 | graphPiCombined->SetMarkerColor(3); | |
2073 | graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined"); | |
2074 | ||
2075 | TList * listCollKa = new TList(); | |
2076 | listCollKa->Add(kaonTpcGraph); | |
2077 | listCollKa->Add(kaonTofGraph); | |
2078 | TGraphErrors * graphKaCombined = new TGraphErrors(); | |
2079 | MergeGraphErrors(graphKaCombined, listCollKa); | |
2080 | graphKaCombined->SetMarkerStyle(22); | |
2081 | graphKaCombined->SetMarkerColor(5); | |
2082 | graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined"); | |
2083 | ||
2084 | TList * listCollEl = new TList(); | |
2085 | if (useV0s) { | |
2086 | //listCollEl->Add(electronGraph); | |
2087 | ||
2088 | TGraphErrors * electronV0GraphCut = new TGraphErrors(*electronV0Graph); | |
2089 | electronV0GraphCut->SetName("electronV0GraphCut"); | |
2090 | /* | |
2091 | for(Int_t ip = 0; ip < electronV0GraphCut->GetN(); ip ++) { | |
2092 | Double_t mom = electronV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron); | |
2093 | Bool_t removePoint = mom > 0.3; | |
2094 | if (removePoint) { | |
2095 | electronV0GraphCut->RemovePoint(ip); | |
2096 | ip--; | |
2097 | continue; | |
2098 | } | |
2099 | } | |
2100 | ||
2101 | TGraphErrors * electronV0plusTOFGraphCut = new TGraphErrors(*electronV0plusTOFGraph); | |
2102 | electronV0plusTOFGraphCut->SetName("electronV0plusTOFGraphCut"); | |
2103 | for(Int_t ip = 0; ip < electronV0plusTOFGraphCut->GetN(); ip ++) { | |
2104 | Double_t mom = electronV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron); | |
2105 | Bool_t removePoint = mom <= 0.4; | |
2106 | ||
2107 | if (removePoint) { | |
2108 | electronV0plusTOFGraphCut->RemovePoint(ip); | |
2109 | ip--; | |
2110 | continue; | |
2111 | } | |
2112 | } | |
2113 | */ | |
2114 | listCollEl->Add(electronV0GraphCut); | |
2115 | //listCollEl->Add(electronV0plusTOFGraphCut); | |
2116 | } | |
2117 | else { | |
2118 | listCollEl->Add(electronGraph); | |
2119 | } | |
2120 | TGraphErrors * graphElCombined = new TGraphErrors(); | |
2121 | MergeGraphErrors(graphElCombined, listCollEl); | |
2122 | graphElCombined->SetMarkerStyle(22); | |
2123 | graphElCombined->SetMarkerColor(2); | |
2124 | graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined"); | |
2125 | ||
2126 | ||
2127 | ||
2128 | // | |
2129 | // return array with all graphs | |
2130 | // | |
2131 | TObjArray * arrGraphs = new TObjArray(); | |
2132 | arrGraphs->Add(graphAll); | |
2133 | // | |
2134 | arrGraphs->Add(electronGraph); | |
2135 | arrGraphs->Add(electronV0Graph); | |
2136 | arrGraphs->Add(electronV0plusTOFGraph); | |
2137 | arrGraphs->Add(electronV0GraphForBBfit); | |
2138 | arrGraphs->Add(pionTpcGraph); | |
2139 | arrGraphs->Add(pionTofGraph); | |
2140 | arrGraphs->Add(pionV0Graph); | |
2141 | arrGraphs->Add(pionV0plusTOFGraph); | |
2142 | arrGraphs->Add(pionV0GraphForBBfit); | |
2143 | arrGraphs->Add(kaonTpcGraph); | |
2144 | arrGraphs->Add(kaonTofGraph); | |
2145 | arrGraphs->Add(protonTpcGraph); | |
2146 | arrGraphs->Add(protonTofGraph); | |
2147 | arrGraphs->Add(protonV0Graph); | |
2148 | arrGraphs->Add(protonV0plusTOFGraph); | |
2149 | arrGraphs->Add(protonV0GraphForBBfit); | |
2150 | // | |
2151 | arrGraphs->Add(graphElCombined); | |
2152 | arrGraphs->Add(graphPiCombined); | |
2153 | arrGraphs->Add(graphKaCombined); | |
2154 | arrGraphs->Add(graphPrCombined); | |
2155 | // | |
2156 | // | |
2157 | ||
2158 | canvasQAtpc->SaveAs("splines_QA_ResidualGraphsTPC.root"); | |
2159 | canvasQAtof->SaveAs("splines_QA_ResidualGraphsTOF.root"); | |
2160 | canvasQAv0->SaveAs("splines_QA_ResidualGraphsV0.root"); | |
2161 | canvasQAv0plusTOF->SaveAs("splines_QA_ResidualGraphsV0plusTOF.root"); | |
2162 | ||
2163 | ||
2164 | canvasQAv0DeDxPurityEl->SaveAs("V0_dEdx_purity_el.root"); | |
2165 | canvasQAv0DeDxPurityPi->SaveAs("V0_dEdx_purity_Pi.root"); | |
2166 | canvasQAv0DeDxPurityPr->SaveAs("V0_dEdx_purity_Pr.root"); | |
2167 | ||
2168 | return arrGraphs; | |
2169 | ||
2170 | } | |
2171 | ||
2172 | ||
2173 | //________________________________________________________________________ | |
2174 | TObjArray * AliTPCcalibResidualPID::GetResidualGraphsMC(THnSparseF * histPidQA, const Char_t * /*system*/) { | |
2175 | // | |
2176 | // Extracts the residual graphs from THnSparse created from MC (NOT DATA!) | |
2177 | // | |
2178 | ||
2179 | const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})"; | |
2180 | const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)"; | |
2181 | ||
2182 | Int_t cutForFitting = 10; | |
2183 | Double_t heightFractionForFittingRange = 0.1; | |
2184 | ||
2185 | // | |
2186 | THnSparse * hist = histPidQA; | |
2187 | // | |
2188 | TCanvas * canvasQAmc = new TCanvas("canvasQAmcResGraph","Control canvas for residual graphs (MC)",100,10,1380,800); | |
2189 | canvasQAmc->Divide(2,2); | |
2190 | ||
2191 | for (Int_t i = 1; i <= 4; i++) { | |
2192 | canvasQAmc->GetPad(i)->SetGrid(1, 1); | |
2193 | canvasQAmc->GetPad(i)->SetLogz(); | |
2194 | canvasQAmc->GetPad(i)->SetLogx(); | |
2195 | } | |
2196 | // | |
2197 | // 1. select and fit electrons | |
2198 | // | |
2199 | // In the following, axis 3 will also be limited in range, although nsigma itself is not used. This is to avoid multiple counting of entries | |
2200 | hist->GetAxis(2)->SetRangeUser(0,0); // electrons | |
2201 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons (used for nsigma and for the eta correction of the signal) | |
2202 | TH2D * histElectron = hist->Projection(1,0); | |
2203 | histElectron->SetName("histElectron"); | |
2204 | histElectron->RebinX(2); | |
2205 | histElectron->GetXaxis()->SetRangeUser(0.15,8.0); | |
2206 | TObjArray arr; | |
2207 | FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); | |
2208 | TH1D * electronPoints = (TH1D *) arr.At(1)->Clone(); | |
2209 | // | |
2210 | TGraphErrors * electronGraph = new TGraphErrors(electronPoints); | |
2211 | electronGraph->SetName("electronGraphMC"); | |
2212 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { | |
2213 | Bool_t removePoint = electronGraph->GetY()[ip] < 10 || electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.03 | |
2214 | || electronGraph->GetX()[ip] < 0.15;// || electronGraph->GetX()[ip] > 3.1; | |
2215 | if (removePoint) { | |
2216 | electronGraph->RemovePoint(ip); | |
2217 | ip--; | |
2218 | continue; | |
2219 | } | |
2220 | } | |
2221 | // | |
2222 | canvasQAmc->cd(1); | |
2223 | histElectron->SetTitle("electrons"); | |
2224 | histElectron->GetYaxis()->SetRangeUser(50, 120); | |
2225 | histElectron->GetXaxis()->SetTitle(momTitle.Data()); | |
2226 | histElectron->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
2227 | histElectron->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2228 | histElectron->GetXaxis()->SetNoExponent(kTRUE); | |
2229 | histElectron->Draw("colz"); | |
2230 | electronPoints->SetMarkerStyle(24); | |
2231 | electronPoints->Draw("same"); | |
2232 | // | |
2233 | // 2. protons | |
2234 | // | |
2235 | hist->GetAxis(2)->SetRangeUser(3,3); // protons | |
2236 | hist->GetAxis(3)->SetRangeUser(3,3); // protons (used for nsigma and for the eta correction of the signal) | |
2237 | // | |
2238 | TH2D * histProton = hist->Projection(1,0); | |
2239 | histProton->SetName("histProton"); | |
2240 | histProton->GetXaxis()->SetRangeUser(0.15,8); | |
2241 | histProton->GetYaxis()->SetRangeUser(10, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins())); | |
2242 | FitSlicesY(histProton, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); | |
2243 | TH1D * protonPoints = (TH1D *) arr.At(1)->Clone(); | |
2244 | // | |
2245 | TGraphErrors * protonGraph = new TGraphErrors(protonPoints); | |
2246 | protonGraph->SetName("protonGraphMC"); | |
2247 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { | |
2248 | Bool_t removePoint = protonGraph->GetY()[ip] < 10 || protonGraph->GetEY()[ip]/protonGraph->GetY()[ip] > 0.02 || | |
2249 | protonGraph->GetX()[ip] < 0.15 || protonGraph->GetX()[ip] > 6; | |
2250 | if (removePoint) { | |
2251 | protonGraph->RemovePoint(ip); | |
2252 | ip--; | |
2253 | continue; | |
2254 | } | |
2255 | } | |
2256 | // | |
2257 | canvasQAmc->cd(2); | |
2258 | histProton->SetTitle("protons"); | |
2259 | histProton->GetXaxis()->SetTitle(momTitle.Data()); | |
2260 | histProton->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
2261 | histProton->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2262 | histProton->GetXaxis()->SetNoExponent(kTRUE); | |
2263 | histProton->Draw("colz"); | |
2264 | protonPoints->SetMarkerStyle(20); | |
2265 | protonPoints->Draw("same"); | |
2266 | ||
2267 | // | |
2268 | // 3. pions | |
2269 | // | |
2270 | hist->GetAxis(2)->SetRangeUser(1,1); // pions | |
2271 | hist->GetAxis(3)->SetRangeUser(1,1); // pions (used for nsigma and for the eta correction of the signal) | |
2272 | // | |
2273 | TH2D * histPion = hist->Projection(1,0); | |
2274 | histPion->SetName("histPion"); | |
2275 | histPion->GetXaxis()->SetRangeUser(0.15,50); | |
2276 | FitSlicesY(histPion, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); | |
2277 | TH1D * pionPoints = (TH1D *) arr.At(1)->Clone(); | |
2278 | // | |
2279 | TGraphErrors * pionGraph = new TGraphErrors(pionPoints); | |
2280 | pionGraph->SetName("pionGraphMC"); | |
2281 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { | |
2282 | Bool_t removePoint = pionGraph->GetY()[ip] < 10 || pionGraph->GetEY()[ip]/pionGraph->GetY()[ip] > 0.005 | |
2283 | || pionGraph->GetX()[ip] < 0.15 || pionGraph->GetX()[ip] > 30; | |
2284 | if (removePoint) { | |
2285 | pionGraph->RemovePoint(ip); | |
2286 | ip--; | |
2287 | continue; | |
2288 | } | |
2289 | } | |
2290 | // | |
2291 | canvasQAmc->cd(3); | |
2292 | histPion->GetYaxis()->SetRangeUser(30, 90); | |
2293 | histPion->GetXaxis()->SetTitle(momTitle.Data()); | |
2294 | histPion->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
2295 | histPion->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2296 | histPion->GetXaxis()->SetNoExponent(kTRUE); | |
2297 | histPion->Draw("colz"); | |
2298 | histPion->SetTitle("pions"); | |
2299 | pionPoints->SetMarkerStyle(20); | |
2300 | pionPoints->Draw("same"); | |
2301 | ||
2302 | // | |
2303 | // 4. kaons | |
2304 | // | |
2305 | hist->GetAxis(2)->SetRangeUser(2,2); // kaons | |
2306 | hist->GetAxis(3)->SetRangeUser(2,2); // kaons (used for nsigma and for the eta correction of the signal) | |
2307 | // | |
2308 | TH2D * histKaon = hist->Projection(1,0); | |
2309 | histKaon->SetName("histKaon"); | |
2310 | histKaon->GetXaxis()->SetRangeUser(0.15,8); | |
2311 | FitSlicesY(histKaon, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); | |
2312 | TH1D * kaonPoints = (TH1D*) arr.At(1)->Clone(); | |
2313 | // | |
2314 | TGraphErrors * kaonGraph = new TGraphErrors(kaonPoints); | |
2315 | kaonGraph->SetName("kaonGraphMC"); | |
2316 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { | |
2317 | Bool_t removePoint = kaonGraph->GetY()[ip] < 10 || kaonGraph->GetEY()[ip]/kaonGraph->GetY()[ip] > 0.02 | |
2318 | || kaonGraph->GetX()[ip] < 0.15; | |
2319 | if (removePoint) { | |
2320 | kaonGraph->RemovePoint(ip); | |
2321 | ip--; | |
2322 | continue; | |
2323 | } | |
2324 | } | |
2325 | // | |
2326 | canvasQAmc->cd(4); | |
2327 | histKaon->SetTitle("kaons"); | |
2328 | histKaon->GetYaxis()->SetRangeUser(30, 500); | |
2329 | histKaon->GetXaxis()->SetTitle(momTitle.Data()); | |
2330 | histKaon->GetYaxis()->SetTitle(dEdxTitle.Data()); | |
2331 | histKaon->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2332 | histKaon->GetXaxis()->SetNoExponent(kTRUE); | |
2333 | histKaon->Draw("colz"); | |
2334 | kaonPoints->SetMarkerStyle(20); | |
2335 | kaonPoints->Draw("same"); | |
2336 | ||
2337 | ||
2338 | ||
2339 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES | |
2340 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES | |
2341 | ||
2342 | ||
2343 | // Clone graphs to have graphs with the same names as for the data case -> same subsequent functions can be used to determine | |
2344 | // the correction functions and, finally, the response functions. | |
2345 | // In addition: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that | |
2346 | // these graphs can be used for the BB fit | |
2347 | TGraphErrors * graphPrCombined = new TGraphErrors(*protonGraph); | |
2348 | graphPrCombined->SetMarkerStyle(22); | |
2349 | graphPrCombined->SetMarkerColor(4); | |
2350 | graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined"); | |
2351 | ||
2352 | TGraphErrors * graphPiCombined = new TGraphErrors(*pionGraph); | |
2353 | graphPiCombined->SetMarkerStyle(22); | |
2354 | graphPiCombined->SetMarkerColor(3); | |
2355 | graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined"); | |
2356 | ||
2357 | TGraphErrors * graphKaCombined = new TGraphErrors(*kaonGraph); | |
2358 | graphKaCombined->SetMarkerStyle(22); | |
2359 | graphKaCombined->SetMarkerColor(5); | |
2360 | graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined"); | |
2361 | ||
2362 | TGraphErrors * graphElCombined = new TGraphErrors(*electronGraph); | |
2363 | graphElCombined->SetMarkerStyle(22); | |
2364 | graphElCombined->SetMarkerColor(2); | |
2365 | graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined"); | |
2366 | ||
2367 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { | |
2368 | Bool_t removePoint = electronGraph->GetX()[ip] < 2.5;// || electronGraph->GetX()[ip] > 5.0; | |
2369 | //Bool_t removePoint = electronGraph->GetX()[ip] < 1.2 || electronGraph->GetX()[ip] > 5.0; | |
2370 | if (removePoint) { | |
2371 | electronGraph->RemovePoint(ip); | |
2372 | ip--; | |
2373 | continue; | |
2374 | } | |
2375 | } | |
2376 | ||
2377 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { | |
2378 | Bool_t removePoint = protonGraph->GetX()[ip] < 0.9 | |
2379 | || protonGraph->GetX()[ip] > 2.5;//3.5; //TODO NEW in order not to get into conflict with pions/kaons | |
2380 | if (removePoint) { | |
2381 | protonGraph->RemovePoint(ip); | |
2382 | ip--; | |
2383 | continue; | |
2384 | } | |
2385 | } | |
2386 | ||
2387 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { | |
2388 | Bool_t removePoint = pionGraph->GetX()[ip] < 2.1;//TODO APR18 0.8; | |
2389 | if (removePoint) { | |
2390 | pionGraph->RemovePoint(ip); | |
2391 | ip--; | |
2392 | continue; | |
2393 | } | |
2394 | } | |
2395 | ||
2396 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { | |
2397 | Bool_t removePoint = kaonGraph->GetX()[ip] < 1.25 || kaonGraph->GetX()[ip] > 4.0;// < 0.7;// || kaonGraph->GetX()[ip] > 4; | |
2398 | if (removePoint) { | |
2399 | kaonGraph->RemovePoint(ip); | |
2400 | ip--; | |
2401 | continue; | |
2402 | } | |
2403 | } | |
2404 | // | |
2405 | // rescale graphs to betaGamma | |
2406 | // | |
2407 | kaonGraph->SetMarkerStyle(22); | |
2408 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { | |
2409 | kaonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
2410 | kaonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
2411 | } | |
2412 | for(Int_t ip = 0; ip < graphKaCombined->GetN(); ip ++) { | |
2413 | graphKaCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
2414 | graphKaCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); | |
2415 | } | |
2416 | // | |
2417 | electronGraph->SetMarkerStyle(22); | |
2418 | electronGraph->SetMarkerColor(2); | |
2419 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { | |
2420 | electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
2421 | electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
2422 | } | |
2423 | for(Int_t ip = 0; ip < graphElCombined->GetN(); ip ++) { | |
2424 | graphElCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
2425 | graphElCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); | |
2426 | } | |
2427 | // | |
2428 | pionGraph->SetMarkerStyle(22); | |
2429 | pionGraph->SetMarkerColor(3); | |
2430 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { | |
2431 | pionGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
2432 | pionGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
2433 | } | |
2434 | for(Int_t ip = 0; ip < graphPiCombined->GetN(); ip ++) { | |
2435 | graphPiCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
2436 | graphPiCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); | |
2437 | } | |
2438 | // | |
2439 | protonGraph->SetMarkerStyle(22); | |
2440 | protonGraph->SetMarkerColor(4); | |
2441 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { | |
2442 | protonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
2443 | protonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
2444 | } | |
2445 | for(Int_t ip = 0; ip < graphPrCombined->GetN(); ip ++) { | |
2446 | graphPrCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
2447 | graphPrCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); | |
2448 | } | |
2449 | // | |
2450 | // | |
2451 | ||
2452 | // Store graphs without further restriction | |
2453 | TGraphErrors * graphPrAll = new TGraphErrors(*graphPrCombined); | |
2454 | graphPrAll->SetNameTitle("Protons_MC_all", "Protons_MC_all"); | |
2455 | ||
2456 | TGraphErrors * graphPiAll = new TGraphErrors(*graphPiCombined); | |
2457 | graphPiAll->SetNameTitle("Pions_MC_all", "Pions_MC_all"); | |
2458 | ||
2459 | TGraphErrors * graphKaAll = new TGraphErrors(*graphKaCombined); | |
2460 | graphKaAll->SetNameTitle("Kaons_MC_all", "Kaons_MC_all"); | |
2461 | ||
2462 | TGraphErrors * graphElAll = new TGraphErrors(*graphElCombined); | |
2463 | graphElAll->SetNameTitle("Electrons_MC_all", "Electrons_MC_all"); | |
2464 | ||
2465 | // | |
2466 | // | |
2467 | TGraphErrors * graphAll = new TGraphErrors(); | |
2468 | TList * listColl = new TList(); | |
2469 | listColl->Add(electronGraph); | |
2470 | listColl->Add(protonGraph); | |
2471 | listColl->Add(kaonGraph); | |
2472 | listColl->Add(pionGraph); | |
2473 | MergeGraphErrors(graphAll, listColl); | |
2474 | graphAll->SetNameTitle("beamDataPoints","beamDataPoints"); | |
2475 | ||
2476 | // | |
2477 | // return array with all graphs | |
2478 | // | |
2479 | TObjArray * arrGraphs = new TObjArray(); | |
2480 | arrGraphs->Add(graphAll); | |
2481 | // | |
2482 | arrGraphs->Add(electronGraph); | |
2483 | arrGraphs->Add(pionGraph); | |
2484 | arrGraphs->Add(kaonGraph); | |
2485 | arrGraphs->Add(protonGraph); | |
2486 | // | |
2487 | arrGraphs->Add(graphElAll); | |
2488 | arrGraphs->Add(graphPiAll); | |
2489 | arrGraphs->Add(graphKaAll); | |
2490 | arrGraphs->Add(graphPrAll); | |
2491 | // | |
2492 | arrGraphs->Add(graphElCombined); | |
2493 | arrGraphs->Add(graphPiCombined); | |
2494 | arrGraphs->Add(graphKaCombined); | |
2495 | arrGraphs->Add(graphPrCombined); | |
2496 | // | |
2497 | // | |
2498 | ||
2499 | canvasQAmc->SaveAs("splines_QA_ResidualGraphsTPC.root"); | |
2500 | ||
2501 | return arrGraphs; | |
2502 | ||
2503 | } | |
2504 | ||
2505 | ||
2506 | //________________________________________________________________________ | |
2507 | TObjArray * AliTPCcalibResidualPID::GetResponseFunctions(TF1* parametrisation, TObjArray* inputGraphs, const Char_t * type, const Char_t * period, const Char_t * pass, const Char_t * system, const Char_t* dedxtype) { | |
2508 | // | |
2509 | // | |
2510 | ||
2511 | Bool_t isMC = kFALSE; | |
2512 | if (strcmp(type, "MC") == 0) { | |
2513 | isMC = kTRUE; | |
2514 | } | |
2515 | else if (strcmp(type, "DATA") == 0) { | |
2516 | isMC = kFALSE; | |
2517 | } | |
2518 | else { | |
2519 | Printf("ERROR - GetResponseFunctions: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type); | |
2520 | ||
2521 | return 0x0; | |
2522 | } | |
2523 | ||
2524 | Bool_t isPbPb = kFALSE; | |
2525 | Bool_t isPPb = kFALSE; | |
2526 | if (strcmp(system, "PBPB") == 0) { | |
2527 | Printf("PbPb detected - Applying special handling!"); | |
2528 | isPbPb = kTRUE; | |
2529 | } | |
2530 | else if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { | |
2531 | Printf("p-Pb/Pb-p detected - Applying special handling!"); | |
2532 | isPPb = kTRUE; | |
2533 | } | |
2534 | ||
2535 | TCanvas* canvasQA[4] = { 0x0, }; | |
2536 | for (Int_t i = 0; i < 4; i++) { | |
2537 | canvasQA[i] = new TCanvas(Form("canvasQAres_%d", i), "Control canvas for residual polynomial",100,10,1380,800); | |
2538 | canvasQA[i]->SetGrid(1, 1); | |
2539 | canvasQA[i]->SetLogx(); | |
2540 | } | |
2541 | /* | |
2542 | TCanvas * canvasQA = new TCanvas("canvasQAres","Control canvas for residual polynomial",100,10,1380,800); | |
2543 | canvasQA->Divide(2,2); | |
2544 | ||
2545 | for (Int_t i = 1; i <= 4; i++) { | |
2546 | canvasQA->GetPad(i)->SetGrid(1, 1); | |
2547 | } | |
2548 | */ | |
2549 | // | |
2550 | // smallest needed (100MeV proton) bg = 0.1/0.938(=AliPID::ParticleMass(AliPID::kProton)) = 0.1, | |
2551 | // largest needed 100 GeV electron, bg = 100/0.000511(=AliPID::ParticleMass(AliPID::kElectron)) = 200 000 | |
2552 | // | |
2553 | Double_t from = 0.1; | |
2554 | Double_t to = 200000; | |
2555 | ||
2556 | TF1* funcBB = new TF1(*parametrisation); | |
2557 | funcBB->SetLineColor(2); | |
2558 | funcBB->SetLineWidth(1); | |
2559 | ||
2560 | TString fitFuncString = "[0] + [1] / x + [2] / (x**2) + [3] / (x**3) + [4] / (x**4) + [5] / (x**5) + [6] * x"; | |
2561 | TString fitFuncString2 = "pol6"; | |
2562 | // | |
2563 | // 1. extract proton corrections | |
2564 | // | |
2565 | TGraphErrors * graphProtonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Protons_Combined"); | |
2566 | TGraphErrors * graphProtonTPCfinal = new TGraphErrors(*graphProtonTPC); | |
2567 | graphProtonTPCfinal->SetName("graphProtonTPCfinal"); | |
2568 | for(Int_t ip = 0; ip < graphProtonTPC->GetN(); ip ++) { | |
2569 | graphProtonTPC->GetY()[ip] -= funcBB->Eval(graphProtonTPC->GetX()[ip]); | |
2570 | graphProtonTPC->GetY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]); | |
2571 | graphProtonTPC->GetEY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);; | |
2572 | } | |
2573 | canvasQA[0]->cd(); | |
2574 | //canvasQA->cd(1); | |
2575 | //gPad->SetLogx(); | |
2576 | graphProtonTPC->GetXaxis()->SetTitle("#beta#gamma"); | |
2577 | graphProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2578 | graphProtonTPC->GetXaxis()->SetNoExponent(kTRUE); | |
2579 | graphProtonTPC->GetYaxis()->SetLabelSize(0.04); | |
2580 | graphProtonTPC->GetYaxis()->SetTitleOffset(0.85); | |
2581 | graphProtonTPC->GetXaxis()->SetTitleOffset(1.0); | |
2582 | graphProtonTPC->GetYaxis()->SetTitle("(data - fit) / fit"); | |
2583 | graphProtonTPC->Draw("ap"); | |
2584 | ||
2585 | if (graphProtonTPC->GetN() <= 0) { | |
2586 | Printf("ERROR - GetResponseFunctions: Proton graph has no data points!"); | |
2587 | ||
2588 | return 0x0; | |
2589 | } | |
2590 | graphProtonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function | |
2591 | TF1 * funcCorrProton = new TF1("funcCorrProton", fitFuncString2.Data(), TMath::Max(graphProtonTPC->GetX()[0], 0.15),1.0); // TODO BEN was 0.18 - 0.85 | |
2592 | graphProtonTPC->Fit(funcCorrProton, "QREX0M"); | |
2593 | // | |
2594 | // 2. extract kaon corrections | |
2595 | // | |
2596 | TGraphErrors * graphKaonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Kaons_Combined"); | |
2597 | TGraphErrors * graphKaonTPCfinal = new TGraphErrors(*graphKaonTPC); | |
2598 | graphKaonTPCfinal->SetName("graphKaonTPCfinal"); | |
2599 | for(Int_t ip = 0; ip < graphKaonTPC->GetN(); ip ++) { | |
2600 | graphKaonTPC->GetY()[ip] -= funcBB->Eval(graphKaonTPC->GetX()[ip]); | |
2601 | graphKaonTPC->GetY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]); | |
2602 | graphKaonTPC->GetEY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);; | |
2603 | } | |
2604 | canvasQA[1]->cd(); | |
2605 | //canvasQA->cd(2); | |
2606 | //gPad->SetLogx(); | |
2607 | graphKaonTPC->GetXaxis()->SetTitle("#beta#gamma"); | |
2608 | graphKaonTPC->GetYaxis()->SetTitle("(data - fit) / fit"); | |
2609 | graphKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2610 | graphKaonTPC->GetXaxis()->SetNoExponent(kTRUE); | |
2611 | graphKaonTPC->GetYaxis()->SetLabelSize(0.04); | |
2612 | graphKaonTPC->GetYaxis()->SetTitleOffset(0.85); | |
2613 | graphKaonTPC->GetXaxis()->SetTitleOffset(1.0); | |
2614 | graphKaonTPC->Draw("ap"); | |
2615 | ||
2616 | if (graphKaonTPC->GetN() <= 0) { | |
2617 | Printf("ERROR - GetResponseFunctions: Kaon graph has no data points!"); | |
2618 | ||
2619 | return 0x0; | |
2620 | } | |
2621 | graphKaonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function | |
2622 | ||
2623 | TF1 * funcCorrKaon = 0x0; | |
2624 | ||
2625 | if (isMC) { | |
2626 | funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), isPPb ? 3.44 : 1.981); | |
2627 | graphKaonTPC->Fit(funcCorrKaon, "QREX0M"); | |
2628 | } | |
2629 | else { | |
2630 | // In case of data there are sometimes problems to fit the shape with one function (could describe the overall deviation, | |
2631 | // but will not cover all the details). | |
2632 | // Nevertheless, this shape (including most of the "details") can be fitted with the following approach with two functions | |
2633 | TF1 * funcCorrKaon1 = new TF1("funcCorrKaon1", fitFuncString.Data(), | |
2634 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); | |
2635 | graphKaonTPC->Fit(funcCorrKaon1, "QREX0M", "same", TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.0); | |
2636 | ||
2637 | TF1 * funcCorrKaon2 = new TF1("funcCorrKaon2", fitFuncString2.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); | |
2638 | graphKaonTPC->Fit(funcCorrKaon2, "QREX0M", "same", (isMC ? 1.981 : 1.0), 1.981); | |
2639 | ||
2640 | funcCorrKaon = new TF1("funcCorrKaon", | |
2641 | "funcCorrKaon1 * 0.5*(1.+TMath::Erf((1 - x) / 0.1)) + ([0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x+[6]*x*x*x*x*x*x) * 0.5*(1.+TMath::Erf((x - 1) / 0.1))", | |
2642 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); | |
2643 | ||
2644 | for (Int_t i = funcCorrKaon1->GetNpar(), j = 0; i < funcCorrKaon1->GetNpar() + funcCorrKaon2->GetNpar(); i++, j++) { | |
2645 | funcCorrKaon->SetParameter(j, funcCorrKaon2->GetParameter(j)); | |
2646 | } | |
2647 | funcCorrKaon->SetLineColor(kRed); | |
2648 | funcCorrKaon->GetHistogram()->DrawClone("csame"); | |
2649 | //funcCorrKaon->Draw("same"); | |
2650 | } | |
2651 | /*TODO | |
2652 | TF1 * funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(),//TODO BEN was fitFuncString2 | |
2653 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), (isMC ? 1.981 : 1.0)); //TODO BEN was 0.79 for data and 1.45 for MC | |
2654 | graphKaonTPC->Fit(funcCorrKaon, "QREX0M"); | |
2655 | */ | |
2656 | // | |
2657 | // 3. extract pion corrections | |
2658 | // | |
2659 | TGraphErrors * graphPionTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Pions_Combined"); | |
2660 | TGraphErrors * graphPionTPCfinal = new TGraphErrors(*graphPionTPC); | |
2661 | graphPionTPCfinal->SetName("graphPionTPCfinal"); | |
2662 | for(Int_t ip = 0; ip < graphPionTPC->GetN(); ip ++) { | |
2663 | graphPionTPC->GetY()[ip] -= funcBB->Eval(graphPionTPC->GetX()[ip]); | |
2664 | graphPionTPC->GetY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]); | |
2665 | graphPionTPC->GetEY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]); | |
2666 | } | |
2667 | canvasQA[2]->cd(); | |
2668 | //canvasQA->cd(3); | |
2669 | //gPad->SetLogx(); | |
2670 | graphPionTPC->GetXaxis()->SetTitle("#beta#gamma"); | |
2671 | graphPionTPC->GetYaxis()->SetTitle("(data - fit) / fit"); | |
2672 | graphPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2673 | graphPionTPC->GetXaxis()->SetNoExponent(kTRUE); | |
2674 | graphPionTPC->GetYaxis()->SetLabelSize(0.04); | |
2675 | graphPionTPC->GetYaxis()->SetTitleOffset(0.85); | |
2676 | graphPionTPC->GetXaxis()->SetTitleOffset(1.0); | |
2677 | graphPionTPC->Draw("ap"); | |
2678 | ||
2679 | if (graphPionTPC->GetN() <= 0) { | |
2680 | Printf("ERROR - GetResponseFunctions: Pion graph has no data points!"); | |
2681 | ||
2682 | return 0x0; | |
2683 | } | |
2684 | graphPionTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function | |
2685 | // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to electrons | |
2686 | TF1 * funcCorrPion = new TF1("funcCorrPion", fitFuncString.Data(), ((isPbPb && !isMC) ? (0.2 / AliPID::ParticleMass(AliPID::kPion)) : graphPionTPC->GetX()[0]), isMC ? (isPPb ? 12.8 : (isPbPb ? 12.8 : 7.1)) : (isPPb ? 7.68 : 7.1)); | |
2687 | graphPionTPC->Fit(funcCorrPion, "QREX0M"); | |
2688 | // | |
2689 | // 4. extract electron corrections | |
2690 | // | |
2691 | TGraphErrors * graphElectronTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Electrons_Combined"); | |
2692 | TGraphErrors * graphElectronTPCfinal = new TGraphErrors(*graphElectronTPC); | |
2693 | graphElectronTPCfinal->SetName("graphElectronTPCfinal"); | |
2694 | for(Int_t ip = 0; ip < graphElectronTPC->GetN(); ip ++) { | |
2695 | graphElectronTPC->GetY()[ip] -= funcBB->Eval(graphElectronTPC->GetX()[ip]); | |
2696 | graphElectronTPC->GetY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]); | |
2697 | graphElectronTPC->GetEY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);; | |
2698 | } | |
2699 | canvasQA[3]->cd(); | |
2700 | //canvasQA->cd(4); | |
2701 | //gPad->SetLogx(); | |
2702 | graphElectronTPC->GetXaxis()->SetTitle("#beta#gamma"); | |
2703 | graphElectronTPC->GetYaxis()->SetTitle("(data - fit) / fit"); | |
2704 | graphElectronTPC->GetXaxis()->SetMoreLogLabels(kTRUE); | |
2705 | graphElectronTPC->GetXaxis()->SetNoExponent(kTRUE); | |
2706 | graphElectronTPC->GetYaxis()->SetLabelSize(0.04); | |
2707 | graphElectronTPC->GetYaxis()->SetTitleOffset(0.85); | |
2708 | graphElectronTPC->GetXaxis()->SetTitleOffset(1.0); | |
2709 | graphElectronTPC->Draw("ap"); | |
2710 | ||
2711 | if (graphElectronTPC->GetN() <= 0) { | |
2712 | Printf("ERROR - GetResponseFunctions: Electron graph has no data points!"); | |
2713 | ||
2714 | return 0x0; | |
2715 | } | |
2716 | graphElectronTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function | |
2717 | // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to pions | |
2718 | TF1 * funcCorrElectron = new TF1("funcCorrElectron", fitFuncString.Data(), (!isMC && isPbPb ? (0.2 / AliPID::ParticleMass(AliPID::kElectron)) :graphElectronTPC->GetX()[0]), (isMC ? 3565 : (isPPb ? 2900 : 1920/*970*/)));// TODO was 1800 for pp data | |
2719 | // NOTE: For data, the results are almost the same for fitFuncString and fitFuncString2. Maybe, fitFuncString2 is slightly better. | |
2720 | graphElectronTPC->Fit(funcCorrElectron, "QREX0M"); | |
2721 | // | |
2722 | // EXTRACT GRAPHS AND PUT THEM TO THE TOBJARRAY | |
2723 | // | |
2724 | const Int_t nBins = 500; | |
2725 | Double_t xBetaGamma[nBins]; | |
2726 | Double_t yProton[nBins]; | |
2727 | Double_t yKaon[nBins]; | |
2728 | Double_t yPion[nBins]; | |
2729 | Double_t yElectron[nBins]; | |
2730 | Double_t yDefault[nBins]; | |
2731 | // | |
2732 | // | |
2733 | // | |
2734 | xBetaGamma[0] = from; | |
2735 | Double_t factor = pow(to/from, 1./nBins); | |
2736 | ||
2737 | // | |
2738 | for(Int_t kk = 0; kk < nBins; kk++) { | |
2739 | if (kk > 0) xBetaGamma[kk] = factor * xBetaGamma[kk-1]; | |
2740 | yProton[kk] = funcBB->Eval(xBetaGamma[kk])/50.; | |
2741 | yPion[kk] = funcBB->Eval(xBetaGamma[kk])/50.; | |
2742 | yKaon[kk] = funcBB->Eval(xBetaGamma[kk])/50.; | |
2743 | yElectron[kk] = funcBB->Eval(xBetaGamma[kk])/50.; | |
2744 | yDefault[kk] = funcBB->Eval(xBetaGamma[kk])/50.; | |
2745 | ||
2746 | // Added by Ben | |
2747 | Double_t widthFactor = 0.020; | |
2748 | Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kProton)/widthFactor) + 1); | |
2749 | Double_t smoothKaon = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kKaon)/widthFactor) + 1); | |
2750 | Double_t smoothPion = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kPion)/widthFactor) + 1); | |
2751 | Double_t smoothElectron = 0.5*(TMath::Erf((funcCorrElectron->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kElectron)/widthFactor) + 1); | |
2752 | ||
2753 | if (xBetaGamma[kk] > funcCorrProton->GetXmax()) | |
2754 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmax())); | |
2755 | // Correction is so large that one runs into trouble at low bg for Protons. | |
2756 | // The deviation is smaller if one takes the lower bound of the correction function | |
2757 | else if (xBetaGamma[kk] < funcCorrProton->GetXmin()) | |
2758 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmin())); | |
2759 | else | |
2760 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk])); | |
2761 | ||
2762 | if (xBetaGamma[kk] > funcCorrKaon->GetXmax()) | |
2763 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmax())); | |
2764 | // Correction is so large that one runs into trouble at low bg for Kaons. | |
2765 | // The deviation is smaller if one takes the lower bound of the correction function | |
2766 | else if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) | |
2767 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmin())); | |
2768 | else | |
2769 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk])); | |
2770 | ||
2771 | if (xBetaGamma[kk] > funcCorrElectron->GetXmax()) | |
2772 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmax())); | |
2773 | else if (xBetaGamma[kk] < funcCorrElectron->GetXmin()) | |
2774 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmin())); | |
2775 | else | |
2776 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(xBetaGamma[kk])); | |
2777 | ||
2778 | // Only true for LHC10d.pass?: Seems not to be needed because any (very small) deviations are most likely due to electron contamination(BEN) | |
2779 | if (xBetaGamma[kk] > funcCorrPion->GetXmax()) | |
2780 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmax())); | |
2781 | else if (xBetaGamma[kk] < funcCorrPion->GetXmin()) | |
2782 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmin())); | |
2783 | else | |
2784 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk])); | |
2785 | ||
2786 | ||
2787 | /* Removed by Ben | |
2788 | Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])/0.002) + 1); | |
2789 | Double_t smoothKaon = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])/0.002) + 1); | |
2790 | Double_t smoothPion = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])/0.002) + 1); | |
2791 | ||
2792 | if (xBetaGamma[kk] < funcCorrProton->GetXmax() && xBetaGamma[kk] > funcCorrProton->GetXmin()) yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk])); | |
2793 | if (xBetaGamma[kk] < funcCorrKaon->GetXmax() && xBetaGamma[kk] > funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk])); | |
2794 | if (xBetaGamma[kk] < funcCorrPion->GetXmax() && xBetaGamma[kk] > funcCorrPion->GetXmin()) yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk])); | |
2795 | //if (xBetaGamma[kk] < funcCorrElectron->GetXmax()) yElectron[kk] *= (1 + funcCorrElectron->Eval(xBetaGamma[kk])); | |
2796 | // | |
2797 | if (xBetaGamma[kk] < funcCorrProton->GetXmin()) yProton[kk] *= (1 + funcCorrProton->Eval(funcCorrProton->GetXmin())); | |
2798 | if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + funcCorrKaon->Eval(funcCorrKaon->GetXmin())); | |
2799 | if (xBetaGamma[kk] < funcCorrPion->GetXmin()) yPion[kk] *= (1 + funcCorrPion->Eval(funcCorrPion->GetXmin())); | |
2800 | */ | |
2801 | } | |
2802 | // | |
2803 | TGraph * graphProton = new TGraph(nBins, xBetaGamma, yProton); | |
2804 | TGraph * graphElectron = new TGraph(nBins, xBetaGamma, yElectron); | |
2805 | TGraph * graphKaon = new TGraph(nBins, xBetaGamma, yKaon); | |
2806 | TGraph * graphPion = new TGraph(nBins, xBetaGamma, yPion); | |
2807 | TGraph * graphDefault = new TGraph(nBins, xBetaGamma, yDefault); | |
2808 | // | |
2809 | // | |
2810 | TSpline3 * splineProton = new TSpline3("splineProton", graphProton); | |
2811 | TSpline3 * splineElectron = new TSpline3("splineElectron", graphElectron); | |
2812 | TSpline3 * splineKaon = new TSpline3("splineKaon", graphKaon); | |
2813 | TSpline3 * splinePion = new TSpline3("splinePion", graphPion); | |
2814 | TSpline3 * splineDefault = new TSpline3("splineDefault", graphDefault); | |
2815 | // | |
2816 | // | |
2817 | splineProton->SetNameTitle(Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), | |
2818 | Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); | |
2819 | splineElectron->SetNameTitle(Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), | |
2820 | Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); | |
2821 | splineKaon->SetNameTitle(Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), | |
2822 | Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); | |
2823 | splinePion->SetNameTitle(Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), | |
2824 | Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); | |
2825 | splineDefault->SetNameTitle(Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), | |
2826 | Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); | |
2827 | // | |
2828 | TObjArray * arrResponse = new TObjArray(); | |
2829 | arrResponse->SetName("TPCpidResponseFunctions"); | |
2830 | arrResponse->AddLast(splineProton); | |
2831 | arrResponse->AddLast(splineElectron); | |
2832 | arrResponse->AddLast(splineKaon); | |
2833 | arrResponse->AddLast(splinePion); | |
2834 | arrResponse->AddLast(splineDefault); | |
2835 | // | |
2836 | ||
2837 | // Draw deviation from final results | |
2838 | for(Int_t ip = 0; ip < graphProtonTPCfinal->GetN(); ip ++) { | |
2839 | graphProtonTPCfinal->GetY()[ip] -= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]); | |
2840 | graphProtonTPCfinal->GetY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]); | |
2841 | graphProtonTPCfinal->GetEY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);; | |
2842 | } | |
2843 | canvasQA[0]->cd(); | |
2844 | //canvasQA->cd(1); | |
2845 | graphProtonTPCfinal->SetMarkerStyle(26); | |
2846 | graphProtonTPCfinal->Draw("psame"); | |
2847 | ||
2848 | for(Int_t ip = 0; ip < graphKaonTPCfinal->GetN(); ip ++) { | |
2849 | graphKaonTPCfinal->GetY()[ip] -= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]); | |
2850 | graphKaonTPCfinal->GetY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]); | |
2851 | graphKaonTPCfinal->GetEY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);; | |
2852 | } | |
2853 | canvasQA[1]->cd(); | |
2854 | //canvasQA->cd(2); | |
2855 | graphKaonTPCfinal->SetMarkerStyle(26); | |
2856 | graphKaonTPCfinal->Draw("psame"); | |
2857 | ||
2858 | for(Int_t ip = 0; ip < graphPionTPCfinal->GetN(); ip ++) { | |
2859 | graphPionTPCfinal->GetY()[ip] -= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]); | |
2860 | graphPionTPCfinal->GetY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]); | |
2861 | graphPionTPCfinal->GetEY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);; | |
2862 | } | |
2863 | canvasQA[2]->cd(); | |
2864 | //canvasQA->cd(3); | |
2865 | graphPionTPCfinal->SetMarkerStyle(26); | |
2866 | graphPionTPCfinal->Draw("psame"); | |
2867 | ||
2868 | for(Int_t ip = 0; ip < graphElectronTPCfinal->GetN(); ip ++) { | |
2869 | graphElectronTPCfinal->GetY()[ip] -= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]); | |
2870 | graphElectronTPCfinal->GetY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]); | |
2871 | graphElectronTPCfinal->GetEY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);; | |
2872 | } | |
2873 | canvasQA[3]->cd(); | |
2874 | //canvasQA->cd(4); | |
2875 | graphElectronTPCfinal->SetMarkerStyle(26); | |
2876 | graphElectronTPCfinal->Draw("psame"); | |
2877 | ||
2878 | TFile* fSave = TFile::Open("splines_QA_ResidualPolynomials.root", "RECREATE"); | |
2879 | fSave->cd(); | |
2880 | for (Int_t i = 0; i < 4; i++) | |
2881 | canvasQA[i]->Write(); | |
2882 | ||
2883 | fSave->Close(); | |
2884 | //canvasQA->SaveAs("splines_QA_ResidualPolynomials.root"); | |
2885 | ||
2886 | ||
2887 | delete funcBB; | |
2888 | ||
2889 | return arrResponse; | |
2890 | } | |
2891 | ||
2892 | ||
2893 | //________________________________________________________________________ | |
2894 | TF1* AliTPCcalibResidualPID::FitBB(TObjArray* inputGraphs, Bool_t isMC, Bool_t isPPb, const Bool_t useV0s, | |
2895 | const Double_t * initialParameters, AliTPCcalibResidualPID::FitType fitType) { | |
2896 | // | |
2897 | // Fit Bethe-Bloch parametrisation to data points (inputGraphs) | |
2898 | // | |
2899 | const Int_t nPar = 6; | |
2900 | TGraphErrors * graphAll = (TGraphErrors *) inputGraphs->FindObject("beamDataPoints"); | |
2901 | // | |
2902 | Float_t from = 0.9; //TODO ADJUST -> Very important | |
2903 | Float_t to = graphAll->GetXaxis()->GetXmax() * 2; | |
2904 | ||
2905 | ||
2906 | ||
2907 | TF1* funcBB = 0x0; | |
2908 | ||
2909 | Double_t parametersBBForward[nPar] = { 0, }; | |
2910 | ||
2911 | TVirtualFitter::SetMaxIterations(5e6); | |
2912 | ||
2913 | if (fitType == AliTPCcalibResidualPID::kSaturatedLund) { | |
2914 | printf("Fit function: Saturated Lund\n"); | |
2915 | ||
2916 | funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar); | |
2917 | //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values | |
2918 | //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation | |
2919 | Double_t parametersBB[nPar] = {61.0, 8.7, 1.86, 0.85, 113.4, -38}; // Yields reasonable results for data and MC ~ all periods | |
2920 | ||
2921 | if (isPPb && !isMC) { | |
2922 | parametersBB[0] = 51.6; | |
2923 | parametersBB[1] = 9.7; | |
2924 | parametersBB[2] = 1.62; | |
2925 | parametersBB[3] = 0.99; | |
2926 | parametersBB[4] = 104.4; | |
2927 | parametersBB[5] = -27.0; | |
2928 | } | |
2929 | if (isMC) { | |
2930 | parametersBB[0] = 41.4; | |
2931 | parametersBB[1] = 8.6; | |
2932 | parametersBB[2] = 2.2; | |
2933 | parametersBB[3] = 0.92; | |
2934 | parametersBB[4] = 90.4; | |
2935 | parametersBB[5] = -20.0; | |
2936 | } | |
2937 | funcBB->SetParameters(parametersBB); | |
2938 | ||
2939 | Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10}; | |
2940 | funcBB->SetParErrors(parameterErrorsBB); | |
2941 | ||
2942 | for (Int_t i = 0; i < nPar; i++) | |
2943 | parametersBBForward[i] = parametersBB[i]; | |
2944 | } | |
2945 | else if (fitType == AliTPCcalibResidualPID::kLund) { | |
2946 | printf("Fit function: Lund\n"); | |
2947 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); | |
2948 | ||
2949 | funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar); | |
2950 | //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values | |
2951 | //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation | |
2952 | Double_t parametersBB[nPar] = {35.0, 7., 1.86, 1., 75., 0}; // Yields reasonable results for data and MC ~ all periods | |
2953 | ||
2954 | funcBB->SetParameters(parametersBB); | |
2955 | ||
2956 | Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10}; | |
2957 | funcBB->SetParErrors(parameterErrorsBB); | |
2958 | ||
2959 | funcBB->FixParameter(5, 0.0); // No saturation | |
2960 | ||
2961 | for (Int_t i = 0; i < nPar; i++) | |
2962 | parametersBBForward[i] = parametersBB[i]; | |
2963 | } | |
2964 | else if (fitType == kAlephWithAdditionalParam) { | |
2965 | printf("Fit function: Aleph with additional parameter\n"); | |
2966 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); | |
2967 | ||
2968 | // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability. | |
2969 | // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus, | |
2970 | // fixed to unity | |
2971 | funcBB = new TF1("funcAleph", | |
2972 | "[0]/TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])*([1]-[2]-[5]*TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to); | |
2973 | //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); | |
2974 | ||
2975 | Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06}; | |
2976 | //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1}; | |
2977 | funcBB->SetParameters(parametersBB); | |
2978 | ||
2979 | for (Int_t i = 0; i < nPar; i++) | |
2980 | parametersBBForward[i] = parametersBB[i]; | |
2981 | } | |
2982 | else if (fitType == AliTPCcalibResidualPID::kAleph) { | |
2983 | printf("Fit function: Aleph\n"); | |
2984 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); | |
2985 | ||
2986 | // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability. | |
2987 | // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus, | |
2988 | // fixed to unity | |
2989 | funcBB = new TF1("funcAleph", | |
2990 | "[0]/TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])*([1]-[2]-[5]*TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to); | |
2991 | //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); | |
2992 | //TEST Double_t parametersBB[nPar] = {1.2, 26.0, -30.0, -2.15, 5.6, 1}; | |
2993 | Double_t parametersBB[nPar] = {1.25, 27.5, -29.0, 2.2, 5.2, 1}; | |
2994 | ||
2995 | // For [5] floating: Double_t parametersBB[nPar] = {2.42,15.2,-16,-2.24,2.8, 0.057}; | |
2996 | //OLD with different sign for [3] and [5] not fix: Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06}; | |
2997 | //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1}; | |
2998 | funcBB->SetParameters(parametersBB); | |
2999 | ||
3000 | // Fix parameter 5 to original value of unity | |
3001 | funcBB->FixParameter(5, 1); | |
3002 | ||
3003 | for (Int_t i = 0; i < nPar; i++) | |
3004 | parametersBBForward[i] = parametersBB[i]; | |
3005 | } | |
3006 | else { | |
3007 | printf("Error: fitType not supported!\n"); | |
3008 | return 0x0; | |
3009 | } | |
3010 | ||
3011 | funcBB->SetLineColor(2); | |
3012 | funcBB->SetLineWidth(1); | |
3013 | ||
3014 | // Override initial parameters, if user provides some | |
3015 | if (initialParameters) { | |
3016 | printf("Setting user initial parameters!\n"); | |
3017 | funcBB->SetParameters(initialParameters); | |
3018 | } | |
3019 | ||
3020 | // | |
3021 | // | |
3022 | // | |
3023 | TGraphErrors * graphDelta = new TGraphErrors(*graphAll); | |
3024 | graphDelta->SetName("graphDelta"); | |
3025 | ||
3026 | // In MC case: Remove errors from fit -> Some data points with extremely small errors otherwise completely dominate the fit, | |
3027 | // but these points could still be wrong due to systematics (low p effects, e.g.) | |
3028 | if (isMC) { | |
3029 | for(Int_t ip = 0; ip < graphAll->GetN(); ip++) { | |
3030 | graphAll->SetPointError(ip, 0, 0); | |
3031 | } | |
3032 | } | |
3033 | else if (useV0s) { | |
3034 | // Increase weight in plateau by reducing errors. Neccessary because errors are large there compared to other data points | |
3035 | for(Int_t ip = 0; ip < graphAll->GetN(); ip++) { | |
3036 | if (graphAll->GetX()[ip] >= 2500) { | |
3037 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); | |
3038 | } | |
3039 | // Same for pions in the very rel. rise | |
3040 | if (graphAll->GetX()[ip] >= 25 && graphAll->GetX()[ip] < 1000) { | |
3041 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); | |
3042 | } | |
3043 | // Same for protons in the MIP region | |
3044 | if (graphAll->GetX()[ip] >= 2 && graphAll->GetX()[ip] < 4) { | |
3045 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); | |
3046 | } | |
3047 | } | |
3048 | } | |
3049 | graphAll->Fit(funcBB, "REX0M"); | |
3050 | funcBB->SetRange(from, to); | |
3051 | funcBB->GetParameters(parametersBBForward); | |
3052 | // | |
3053 | // | |
3054 | // | |
3055 | for(Int_t ip = 0; ip < graphDelta->GetN(); ip++) { | |
3056 | graphDelta->GetY()[ip] -= funcBB->Eval(graphDelta->GetX()[ip]); | |
3057 | graphDelta->GetY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]); | |
3058 | graphDelta->GetEY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]); | |
3059 | } | |
3060 | TCanvas * canvDelta_1 = new TCanvas("canvDelta_1","control histogram for Bethe-Bloch fit 1",100,10,1380,800); | |
3061 | TCanvas * canvDelta_2 = new TCanvas("canvDelta_2","control histogram for Bethe-Bloch fit 2",100,10,1380,800); | |
3062 | canvDelta_1->SetGrid(1, 1); | |
3063 | canvDelta_1->SetLogx(); | |
3064 | canvDelta_2->SetGrid(1, 1); | |
3065 | canvDelta_2->SetLogx(); | |
3066 | /* | |
3067 | canvDelta->Divide(2, 1); | |
3068 | canvDelta->GetPad(1)->SetGrid(1, 1); | |
3069 | canvDelta->GetPad(1)->SetLogx(); | |
3070 | canvDelta->GetPad(2)->SetGrid(1, 1); | |
3071 | canvDelta->GetPad(2)->SetLogx(); | |
3072 | */ | |
3073 | ||
3074 | canvDelta_1->cd(); | |
3075 | //canvDelta->cd(1); | |
3076 | TH1F *hBBdummy=new TH1F("hBBdummy","BB fit;#beta#gamma;#LTdE/dx#GT (arb. unit)",100,.8,1e4); | |
3077 | hBBdummy->SetMinimum(45); | |
3078 | hBBdummy->SetMaximum(120); | |
3079 | hBBdummy->GetXaxis()->SetTitleOffset(1.1); | |
3080 | hBBdummy->SetStats(kFALSE); | |
3081 | hBBdummy->Draw(); | |
3082 | ||
3083 | graphAll->SetTitle("BB multi-graph fit"); | |
3084 | graphAll->SetMarkerStyle(22); | |
3085 | graphAll->SetMarkerColor(kMagenta); | |
3086 | graphAll->Draw("p"); | |
3087 | ||
3088 | TLegend *leg=new TLegend(.7,.12,.89,isMC?.4:.6); | |
3089 | leg->SetFillColor(10); | |
3090 | leg->SetBorderSize(1); | |
3091 | ||
3092 | //overlay individual graphs | |
3093 | TGraph *gr=0x0; | |
3094 | if (isMC) { | |
3095 | gr=(TGraph*)inputGraphs->FindObject("Protons_MC_all"); | |
3096 | gr->SetMarkerColor(kRed); | |
3097 | gr->SetMarkerStyle(24); | |
3098 | gr->Draw("p"); | |
3099 | leg->AddEntry(gr,"p (MC)","p"); | |
3100 | gr=(TGraph*)inputGraphs->FindObject("Pions_MC_all"); | |
3101 | gr->SetMarkerColor(kBlue); | |
3102 | gr->SetMarkerStyle(24); | |
3103 | gr->Draw("p"); | |
3104 | leg->AddEntry(gr,"#pi (MC)","p"); | |
3105 | gr=(TGraph*)inputGraphs->FindObject("Kaons_MC_all"); | |
3106 | gr->SetMarkerColor(kGreen); | |
3107 | gr->SetMarkerStyle(24); | |
3108 | gr->Draw("p"); | |
3109 | leg->AddEntry(gr,"K (MC)","p"); | |
3110 | gr=(TGraph*)inputGraphs->FindObject("Electrons_MC_all"); | |
3111 | gr->SetMarkerColor(kBlack); | |
3112 | gr->SetMarkerStyle(24); | |
3113 | gr->Draw("p"); | |
3114 | leg->AddEntry(gr,"e (MC)","p"); | |
3115 | } | |
3116 | else { | |
3117 | gr=(TGraph*)inputGraphs->FindObject("protonTpcGraph"); | |
3118 | gr->SetMarkerColor(kRed); | |
3119 | gr->SetMarkerStyle(26); | |
3120 | gr->Draw("p"); | |
3121 | leg->AddEntry(gr,"p (TPC)","p"); | |
3122 | gr=(TGraph*)inputGraphs->FindObject("protonTofGraph"); | |
3123 | gr->SetMarkerColor(kRed); | |
3124 | gr->SetMarkerStyle(25); | |
3125 | gr->Draw("p"); | |
3126 | leg->AddEntry(gr,"p (TOF)","p"); | |
3127 | gr=(TGraph*)inputGraphs->FindObject("protonV0Graph");//ForBBfit"); | |
3128 | gr->SetMarkerColor(kRed); | |
3129 | gr->SetMarkerStyle(24); | |
3130 | gr->Draw("p"); | |
3131 | leg->AddEntry(gr,"p (V0)","p"); | |
3132 | gr=(TGraph*)inputGraphs->FindObject("protonV0plusTOFGraph");//ForBBfit"); | |
3133 | gr->SetMarkerColor(kRed); | |
3134 | gr->SetMarkerStyle(30); | |
3135 | gr->Draw("p"); | |
3136 | leg->AddEntry(gr,"p (V0+TOF)","p"); | |
3137 | gr=(TGraph*)inputGraphs->FindObject("pionTpcGraph"); | |
3138 | gr->SetMarkerColor(kBlue); | |
3139 | gr->SetMarkerStyle(26); | |
3140 | gr->Draw("p"); | |
3141 | leg->AddEntry(gr,"#pi (TPC)","p"); | |
3142 | gr=(TGraph*)inputGraphs->FindObject("pionTofGraph"); | |
3143 | gr->SetMarkerColor(kBlue); | |
3144 | gr->SetMarkerStyle(25); | |
3145 | gr->Draw("p"); | |
3146 | leg->AddEntry(gr,"#pi (TOF)","p"); | |
3147 | gr=(TGraph*)inputGraphs->FindObject("pionV0Graph");//ForBBfit"); | |
3148 | gr->SetMarkerColor(kBlue); | |
3149 | gr->SetMarkerStyle(24); | |
3150 | gr->Draw("p"); | |
3151 | leg->AddEntry(gr,"#pi (V0)","p"); | |
3152 | gr=(TGraph*)inputGraphs->FindObject("pionV0plusTOFGraph");//ForBBfit"); | |
3153 | gr->SetMarkerColor(kBlue); | |
3154 | gr->SetMarkerStyle(30); | |
3155 | gr->Draw("p"); | |
3156 | leg->AddEntry(gr,"#pi (V0+TOF)","p"); | |
3157 | gr=(TGraph*)inputGraphs->FindObject("kaonTpcGraph"); | |
3158 | gr->SetMarkerColor(kGreen); | |
3159 | gr->SetMarkerStyle(26); | |
3160 | gr->Draw("p"); | |
3161 | leg->AddEntry(gr,"K (TPC)","p"); | |
3162 | gr=(TGraph*)inputGraphs->FindObject("kaonTofGraph"); | |
3163 | gr->SetMarkerColor(kGreen); | |
3164 | gr->SetMarkerStyle(25); | |
3165 | gr->Draw("p"); | |
3166 | leg->AddEntry(gr,"K (TOF)","p"); | |
3167 | gr=(TGraph*)inputGraphs->FindObject("electronGraph"); | |
3168 | gr->SetMarkerColor(kBlack); | |
3169 | gr->SetMarkerStyle(25); | |
3170 | gr->Draw("p"); | |
3171 | leg->AddEntry(gr,"e","p"); | |
3172 | gr=(TGraph*)inputGraphs->FindObject("electronV0Graph");//ForBBfit"); | |
3173 | gr->SetMarkerColor(kBlack); | |
3174 | gr->SetMarkerStyle(24); | |
3175 | gr->Draw("p"); | |
3176 | leg->AddEntry(gr,"e (V0)","p"); | |
3177 | gr=(TGraph*)inputGraphs->FindObject("electronV0plusTOFGraph");//ForBBfit"); | |
3178 | gr->SetMarkerColor(kBlack); | |
3179 | gr->SetMarkerStyle(30); | |
3180 | gr->Draw("p"); | |
3181 | leg->AddEntry(gr,"e (V0+TOF)","p"); | |
3182 | } | |
3183 | ||
3184 | graphAll->GetXaxis()->SetTitle("#beta#gamma"); | |
3185 | graphAll->GetYaxis()->SetTitle("<dE/dx> (arb. unit)"); | |
3186 | graphAll->Draw("p"); | |
3187 | funcBB->GetHistogram()->DrawClone("csame"); | |
3188 | ||
3189 | leg->AddEntry(graphAll,"used","p"); | |
3190 | leg->AddEntry(funcBB,"fit","l"); | |
3191 | leg->Draw("same"); | |
3192 | ||
3193 | canvDelta_2->cd(); | |
3194 | //canvDelta->cd(2); | |
3195 | TH1F *hBBresdummy=new TH1F("hBBresdummy","residuals of BB fit;#beta#gamma;(data-fit)/fit",100,.8,1e4); | |
3196 | hBBresdummy->SetMinimum(-0.04); | |
3197 | hBBresdummy->SetMaximum(0.04); | |
3198 | hBBresdummy->GetXaxis()->SetTitleOffset(1.1); | |
3199 | hBBresdummy->GetYaxis()->SetTitleOffset(0.8); | |
3200 | hBBresdummy->SetStats(kFALSE); | |
3201 | hBBresdummy->Draw(); | |
3202 | ||
3203 | graphDelta->SetTitle("residuals of BB fit to multi-graph"); | |
3204 | graphDelta->GetXaxis()->SetTitle("#beta#gamma"); | |
3205 | graphDelta->GetYaxis()->SetTitle("(data - fit) / fit"); | |
3206 | graphDelta->SetMarkerStyle(22); | |
3207 | graphDelta->SetMarkerColor(4); | |
3208 | graphDelta->Draw("p"); | |
3209 | ||
3210 | TFile* fSave = TFile::Open("splines_QA_BetheBlochFit.root", "RECREATE"); | |
3211 | fSave->cd(); | |
3212 | canvDelta_1->Write(); | |
3213 | canvDelta_2->Write(); | |
3214 | ||
3215 | TString fitResults = "Fit results:\n"; | |
3216 | for (Int_t i = 0; i < nPar; i++) { | |
3217 | fitResults.Append(Form("par%d:\t%f +- %f\n", i, funcBB->GetParameters()[i], funcBB->GetParErrors()[i])); | |
3218 | } | |
3219 | ||
3220 | TNamed* settings = new TNamed(fitResults.Data(), fitResults.Data()); | |
3221 | settings->Write(); | |
3222 | fSave->Close(); | |
3223 | ||
3224 | ||
3225 | printf("\n\n%s", fitResults.Data()); | |
3226 | ||
3227 | return funcBB; | |
3228 | } | |
3229 | ||
3230 | ||
3231 | //________________________________________________________________________ | |
3232 | Double_t AliTPCcalibResidualPID::Lund(Double_t* xx, Double_t* par) | |
3233 | { | |
3234 | // bg is beta*gamma | |
3235 | const Double_t bg = xx[0]; | |
3236 | ||
3237 | const Double_t beta2 = bg*bg / (1.0 + bg*bg); | |
3238 | ||
3239 | const Double_t a = par[0]; | |
3240 | const Double_t b = par[1]; | |
3241 | const Double_t c = par[2]; | |
3242 | const Double_t e = par[3]; | |
3243 | const Double_t f = par[4]; | |
3244 | ||
3245 | const Double_t d = TMath::Exp(c*(a - f) / b); | |
3246 | ||
3247 | const Double_t powbg = TMath::Power(1.0 + bg, c); | |
3248 | ||
3249 | const Double_t value = a / TMath::Power(beta2,e) + | |
3250 | b/c * TMath::Log(powbg / (1.0 + d*powbg)); | |
3251 | ||
3252 | return value; | |
3253 | } | |
3254 | ||
3255 | ||
3256 | //________________________________________________________________________ | |
3257 | Double_t AliTPCcalibResidualPID::SaturatedLund(Double_t* xx, Double_t* par) | |
3258 | { | |
3259 | const Double_t qq = Lund(xx, par); | |
3260 | return qq * TMath::Exp(par[5] / qq); | |
3261 | } | |
3262 | ||
3263 | ||
3264 | //________________________________________________________________________ | |
3265 | void AliTPCcalibResidualPID::FitSlicesY(TH2 *hist, Double_t heightFractionForRange, Int_t cutThreshold, TString fitOption, TObjArray *arr) | |
3266 | { | |
3267 | //heightPercentageForRange | |
3268 | // custom slices fit | |
3269 | // | |
3270 | ||
3271 | if (!hist) return; | |
3272 | if (!arr) return; | |
3273 | ||
3274 | // If old style is to be used | |
3275 | /* | |
3276 | hist->FitSlicesY(0, 0, -1, cutThreshold, fitOption.Data(), &arr); | |
3277 | return; | |
3278 | */ | |
3279 | ||
3280 | ||
3281 | arr->Clear(); | |
3282 | arr->SetOwner(); | |
3283 | arr->Expand(4); | |
3284 | ||
3285 | TAxis *axis=hist->GetXaxis(); | |
3286 | const TArrayD *bins = axis->GetXbins(); | |
3287 | TH1D** hList = new TH1D*[4]; | |
3288 | ||
3289 | for (Int_t i = 0; i < 4; i++) { | |
3290 | delete gDirectory->FindObject(Form("%s_%d", hist->GetName(), i)); | |
3291 | ||
3292 | if (bins->fN == 0) { | |
3293 | hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", | |
3294 | hist->GetNbinsX(), axis->GetXmin(), axis->GetXmax()); | |
3295 | } else { | |
3296 | hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", hist->GetNbinsX(), bins->fArray); | |
3297 | } | |
3298 | ||
3299 | (*arr)[i] = hList[i]; | |
3300 | } | |
3301 | ||
3302 | for (Int_t ibin=axis->GetFirst(); ibin<=axis->GetLast(); ++ibin){ | |
3303 | TH1 *h=hist->ProjectionY("_temp",ibin,ibin); | |
3304 | if (!h) | |
3305 | continue; | |
3306 | ||
3307 | if (h->GetEntries() < cutThreshold) { | |
3308 | delete h; | |
3309 | continue; | |
3310 | } | |
3311 | ||
3312 | // Average around maximum bin -> Might catch some outliers | |
3313 | Int_t maxBin = h->GetMaximumBin(); | |
3314 | Double_t maxVal = h->GetBinContent(maxBin); | |
3315 | ||
3316 | if (maxVal < 2) { // It could happen that all entries are in overflow/underflow; don't fit in this case | |
3317 | delete h; | |
3318 | continue; | |
3319 | } | |
3320 | ||
3321 | UChar_t usedBins = 1; | |
3322 | if (maxBin > 1) { | |
3323 | maxVal += h->GetBinContent(maxBin - 1); | |
3324 | usedBins++; | |
3325 | } | |
3326 | if (maxBin < h->GetNbinsX()) { | |
3327 | maxVal += h->GetBinContent(maxBin + 1); | |
3328 | usedBins++; | |
3329 | } | |
3330 | maxVal /= usedBins; | |
3331 | ||
3332 | Double_t thresholdFraction = heightFractionForRange * maxVal; | |
3333 | Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction)); | |
3334 | Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction)); | |
3335 | ||
3336 | Double_t lowThreshold = h->GetBinCenter(lowThrBin); | |
3337 | Double_t highThreshold = h->GetBinCenter(highThrBin); | |
3338 | ||
3339 | TFitResultPtr res = h->Fit("gaus", Form("%sS", fitOption.Data()), "", lowThreshold, highThreshold); | |
3340 | ||
3341 | if ((Int_t)res == 0) { | |
3342 | Int_t resBin = ibin; | |
3343 | hList[0]->SetBinContent(resBin,res->GetParams()[0]); | |
3344 | hList[0]->SetBinError(resBin,res->GetErrors()[0]); | |
3345 | hList[1]->SetBinContent(resBin,res->GetParams()[1]); | |
3346 | hList[1]->SetBinError(resBin,res->GetErrors()[1]); | |
3347 | hList[2]->SetBinContent(resBin,res->GetParams()[2]); | |
3348 | hList[2]->SetBinError(resBin,res->GetErrors()[2]); | |
3349 | hList[3]->SetBinContent(resBin,res->Ndf()>0?res->Chi2()/res->Ndf():0); | |
3350 | } | |
3351 | ||
3352 | delete h; | |
3353 | } | |
3354 | } | |
3355 | ||
3356 | ||
3357 | //______________________________________________________________________________ | |
3358 | Bool_t AliTPCcalibResidualPID::GetVertexIsOk(AliVEvent* event) const | |
3359 | { | |
3360 | AliAODEvent* aod = 0x0; | |
3361 | AliESDEvent* esd = 0x0; | |
3362 | ||
3363 | aod = dynamic_cast<AliAODEvent*>(event); | |
3364 | if (!aod) { | |
3365 | esd = dynamic_cast<AliESDEvent*>(event); | |
3366 | ||
3367 | if (!esd) { | |
3368 | AliError("Event seems to be neither AOD nor ESD!"); | |
3369 | return kFALSE; | |
3370 | } | |
3371 | } | |
3372 | ||
3373 | if (fIsPbpOrpPb) { | |
3374 | const AliVVertex* trkVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : | |
3375 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex())); | |
3376 | ||
3377 | if (!trkVtx || trkVtx->GetNContributors() <= 0) | |
3378 | return kFALSE; | |
3379 | ||
3380 | TString vtxTtl = trkVtx->GetTitle(); | |
3381 | if (!vtxTtl.Contains("VertexerTracks")) | |
3382 | return kFALSE; | |
3383 | ||
3384 | Float_t zvtx = trkVtx->GetZ(); | |
3385 | const AliVVertex* spdVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) : | |
3386 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD())); | |
3387 | if (spdVtx->GetNContributors() <= 0) | |
3388 | return kFALSE; | |
3389 | ||
3390 | TString vtxTyp = spdVtx->GetTitle(); | |
3391 | Double_t cov[6] = {0}; | |
3392 | spdVtx->GetCovarianceMatrix(cov); | |
3393 | Double_t zRes = TMath::Sqrt(cov[5]); | |
3394 | if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25)) | |
3395 | return kFALSE; | |
3396 | ||
3397 | if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5) | |
3398 | return kFALSE; | |
3399 | ||
3400 | if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm | |
3401 | return kFALSE; | |
3402 | ||
3403 | return kTRUE; | |
3404 | } | |
3405 | ||
3406 | ||
3407 | // pp and PbPb | |
3408 | const AliVVertex* primaryVertex = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : | |
3409 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks())); | |
3410 | ||
3411 | if (!primaryVertex || primaryVertex->GetNContributors() <= 0) | |
3412 | return kFALSE; | |
3413 | ||
3414 | if (TMath::Abs(primaryVertex->GetZ()) >= fZvtxCutEvent) //Default: 10 cm | |
3415 | return kFALSE; | |
3416 | ||
3417 | return kTRUE; | |
3418 | } | |
3419 | ||
3420 | ||
3421 | //______________________________________________________________________________ | |
3422 | void AliTPCcalibResidualPID::FillV0PIDlist(AliESDEvent* event) | |
3423 | { | |
3424 | // | |
3425 | // Fill the PID tag list | |
3426 | // | |
3427 | ||
3428 | // If no event forwarded as parameter (default), cast current input event. | |
3429 | // Dynamic cast to ESD events (DO NOTHING for AOD events) | |
3430 | if (!event) | |
3431 | event = dynamic_cast<AliESDEvent *>(InputEvent()); | |
3432 | ||
3433 | // If this fails, do nothing | |
3434 | if (!event) | |
3435 | return; | |
3436 | ||
3437 | if (!fV0KineCuts) { | |
3438 | AliError("V0KineCuts not available!"); | |
3439 | return; | |
3440 | } | |
3441 | ||
3442 | TString beamType(event->GetBeamType()); | |
3443 | ||
3444 | if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) { | |
3445 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb); | |
3446 | } | |
3447 | else { | |
3448 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP); | |
3449 | } | |
3450 | ||
3451 | // V0 selection | |
3452 | // set event | |
3453 | fV0KineCuts->SetEvent(event); | |
3454 | ||
3455 | const Int_t numTracks = event->GetNumberOfTracks(); | |
3456 | fV0tags = new Char_t[numTracks]; | |
3457 | for (Int_t i = 0; i < numTracks; i++) | |
3458 | fV0tags[i] = 0; | |
3459 | ||
3460 | fV0motherIndex = new Int_t[numTracks]; | |
3461 | for (Int_t i = 0; i < numTracks; i++) | |
3462 | fV0motherIndex[i] = -1; | |
3463 | ||
3464 | fV0motherPDG = new Int_t[numTracks]; | |
3465 | for (Int_t i = 0; i < numTracks; i++) | |
3466 | fV0motherPDG[i] = 0; | |
3467 | ||
3468 | fNumTagsStored = numTracks; | |
3469 | ||
3470 | ||
3471 | ||
3472 | ||
3473 | ||
3474 | // loop over V0 particles | |
3475 | for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) { | |
3476 | AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0); | |
3477 | ||
3478 | if (!v0) | |
3479 | continue; | |
3480 | ||
3481 | // Reject onFly V0's <-> Only take V0's from offline V0 finder | |
3482 | if (v0->GetOnFlyStatus()) | |
3483 | continue; | |
3484 | ||
3485 | // Get the particle selection | |
3486 | Bool_t foundV0 = kFALSE; | |
3487 | Int_t pdgV0 = 0, pdgP = 0, pdgN = 0; | |
3488 | ||
3489 | foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN); | |
3490 | if (!foundV0) | |
3491 | continue; | |
3492 | ||
3493 | Int_t iTrackP = v0->GetPindex(); // positive track | |
3494 | Int_t iTrackN = v0->GetNindex(); // negative track | |
3495 | ||
3496 | /* | |
3497 | AliESDtrack* trackP = event->GetTrack(iTrackP); | |
3498 | AliESDtrack* trackN = event->GetTrack(iTrackN); | |
3499 | ||
3500 | if (!trackP || !trackN) | |
3501 | continue; | |
3502 | ||
3503 | ||
3504 | ||
3505 | Float_t xy = 999, z = 999; | |
3506 | trackP->GetImpactParameters(xy, z); | |
3507 | ||
3508 | Bool_t reject = kFALSE; | |
3509 | if (TMath::Abs(z) > 1) | |
3510 | continue; | |
3511 | ||
3512 | trackN->GetImpactParameters(xy, z); | |
3513 | if (TMath::Abs(z) > 1) | |
3514 | continue; | |
3515 | */ | |
3516 | ||
3517 | ||
3518 | // Fill the Object arrays | |
3519 | // positive particles | |
3520 | if (pdgP == -11) { | |
3521 | fV0tags[iTrackP] = 14; | |
3522 | } | |
3523 | else if (pdgP == 211) { | |
3524 | fV0tags[iTrackP] = 15; | |
3525 | } | |
3526 | else if(pdgP == 2212) { | |
3527 | fV0tags[iTrackP] = 16; | |
3528 | } | |
3529 | ||
3530 | fV0motherIndex[iTrackP] = iv0; | |
3531 | fV0motherPDG[iTrackP] = pdgV0; | |
3532 | ||
3533 | // negative particles | |
3534 | if( pdgN == 11){ | |
3535 | fV0tags[iTrackN] = -14; | |
3536 | } | |
3537 | else if( pdgN == -211){ | |
3538 | fV0tags[iTrackN] = -15; | |
3539 | } | |
3540 | else if( pdgN == -2212){ | |
3541 | fV0tags[iTrackN] = -16; | |
3542 | } | |
3543 | ||
3544 | fV0motherIndex[iTrackN] = iv0; | |
3545 | fV0motherPDG[iTrackN] = pdgV0; | |
3546 | ||
3547 | } | |
3548 | } | |
3549 | ||
3550 | ||
3551 | //______________________________________________________________________________ | |
3552 | void AliTPCcalibResidualPID::ClearV0PIDlist() | |
3553 | { | |
3554 | // | |
3555 | // Clear the PID tag list | |
3556 | // | |
3557 | ||
3558 | delete fV0tags; | |
3559 | fV0tags = 0; | |
3560 | ||
3561 | delete fV0motherIndex; | |
3562 | fV0motherIndex = 0; | |
3563 | ||
3564 | delete fV0motherPDG; | |
3565 | fV0motherPDG = 0; | |
3566 | ||
3567 | fNumTagsStored = 0; | |
3568 | } | |
3569 | ||
3570 | ||
3571 | //______________________________________________________________________________ | |
3572 | Char_t AliTPCcalibResidualPID::GetV0tag(Int_t trackIndex) const | |
3573 | { | |
3574 | // | |
3575 | // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list. | |
3576 | // | |
3577 | ||
3578 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags) | |
3579 | return -99; | |
3580 | ||
3581 | return fV0tags[trackIndex]; | |
3582 | } | |
3583 | ||
3584 | ||
3585 | //______________________________________________________________________________ | |
3586 | Int_t AliTPCcalibResidualPID::GetV0motherIndex(Int_t trackIndex) const | |
3587 | { | |
3588 | // | |
3589 | // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list. | |
3590 | // | |
3591 | ||
3592 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex) | |
3593 | return -99; | |
3594 | ||
3595 | return fV0motherIndex[trackIndex]; | |
3596 | } | |
3597 | ||
3598 | ||
3599 | //______________________________________________________________________________ | |
3600 | Int_t AliTPCcalibResidualPID::GetV0motherPDG(Int_t trackIndex) const | |
3601 | { | |
3602 | // | |
3603 | // Get the PDG code of the V0 mother for the corresponding trackIndex. Returns 0 in case of invalid index/mother index list. | |
3604 | // | |
3605 | ||
3606 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherPDG) | |
3607 | return 0; | |
3608 | ||
3609 | return fV0motherPDG[trackIndex]; | |
3610 | } | |
3611 | ||
3612 | ||
3613 | //______________________________________________________________________________ | |
3614 | Int_t AliTPCcalibResidualPID::MergeGraphErrors(TGraphErrors* mergedGraph, TCollection* li) | |
3615 | { | |
3616 | // Adds all graphs with errors from the collection to this graph. | |
3617 | // Returns the total number of poins in the result or -1 in case of an error. | |
3618 | ||
3619 | // TODO "Normal" merge of latest root trunk will also take into account the error bars, | |
3620 | // so this function can be replaced by the normal root function with the latest root version. | |
3621 | ||
3622 | TIter next(li); | |
3623 | while (TObject* o = next()) { | |
3624 | TGraph *g = dynamic_cast<TGraph*>(o); | |
3625 | if (!g) { | |
3626 | Printf("ERROR: Cannot merge an object which doesn't inherit from TGraph found in the list"); | |
3627 | return -1; | |
3628 | } | |
3629 | int n0 = mergedGraph->GetN(); | |
3630 | int n1 = n0+g->GetN(); | |
3631 | mergedGraph->Set(n1); | |
3632 | Double_t * x = g->GetX(); | |
3633 | Double_t * y = g->GetY(); | |
3634 | Double_t * ex = g->GetEX(); | |
3635 | Double_t * ey = g->GetEY(); | |
3636 | for (Int_t i = 0 ; i < g->GetN(); i++) { | |
3637 | mergedGraph->SetPoint(n0+i, x[i], y[i]); | |
3638 | Double_t exPoint = ex ? ex[i] : 0; | |
3639 | Double_t eyPoint = ey ? ey[i] : 0; | |
3640 | mergedGraph->SetPointError(n0+i, exPoint, eyPoint); | |
3641 | } | |
3642 | } | |
3643 | return mergedGraph->GetN(); | |
3644 | } | |
3645 | ||
3646 | //________________________________________________________________________ | |
3647 | Bool_t AliTPCcalibResidualPID::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer) | |
3648 | { | |
3649 | // | |
3650 | // TPC Cut MIGeo | |
3651 | // | |
3652 | ||
3653 | if (!track || !evt) | |
3654 | return kFALSE; | |
3655 | ||
3656 | const Short_t sign = track->Charge(); | |
3657 | Double_t xyz[50]; | |
3658 | Double_t pxpypz[50]; | |
3659 | Double_t cv[100]; | |
3660 | ||
3661 | track->GetXYZ(xyz); | |
3662 | track->GetPxPyPz(pxpypz); | |
3663 | ||
3664 | AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign); | |
3665 | const AliESDtrack dummy; | |
3666 | ||
3667 | const Double_t magField = evt->GetMagneticField(); | |
3668 | Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0); | |
3669 | Double_t varNcr = track->GetTPCClusterInfo(3, 1); | |
3670 | Double_t varNcls = track->GetTPCsignalN(); | |
3671 | ||
3672 | const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt()); | |
3673 | Bool_t cutGeom = varGeom > fgCutGeo * varEval; | |
3674 | Bool_t cutNcr = varNcr > fgCutNcr * varEval; | |
3675 | Bool_t cutNcls = varNcls > fgCutNcl * varEval; | |
3676 | ||
3677 | Bool_t kout = cutGeom && cutNcr && cutNcls; | |
3678 | ||
3679 | if (streamer) { | |
3680 | Double_t dedx = track->GetTPCsignal(); | |
3681 | ||
3682 | (*streamer)<<"tree"<< | |
3683 | "param.="<< par<< | |
3684 | "varGeom="<<varGeom<< | |
3685 | "varNcr="<<varNcr<< | |
3686 | "varNcls="<<varNcls<< | |
3687 | ||
3688 | "cutGeom="<<cutGeom<< | |
3689 | "cutNcr="<<cutNcr<< | |
3690 | "cutNcls="<<cutNcls<< | |
3691 | ||
3692 | "kout="<<kout<< | |
3693 | "dedx="<<dedx<< | |
3694 | ||
3695 | "\n"; | |
3696 | } | |
3697 | ||
3698 | delete par; | |
3699 | ||
3700 | return kout; | |
3701 | } |