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