]>
Commit | Line | Data |
---|---|---|
0b7de667 | 1 | #include "AliFMDMCHitEnergyFitter.h" |
2 | #include "AliMCEvent.h" | |
3 | #include "AliESDEvent.h" | |
4 | #include "AliStack.h" | |
5 | #include "AliFMDHit.h" | |
6 | #include "AliESDFMD.h" | |
7 | #include "AliMCAuxHandler.h" | |
8 | #include "AliGenEventHeader.h" | |
9 | #include "AliLog.h" | |
10 | #include "AliFMDEncodedEdx.h" | |
11 | #include <TH1.h> | |
12 | #include <TH2.h> | |
13 | #include <TNtuple.h> | |
14 | #include <TMath.h> | |
15 | ||
16 | ||
17 | //____________________________________________________________________ | |
18 | AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter() | |
19 | : AliFMDEnergyFitter(), | |
20 | fSumPrimary(0), | |
21 | fSumSecondary(0), | |
22 | fIp(), | |
23 | fNTrack(0), | |
24 | fNPrimary(0), | |
25 | fTuple(0) | |
26 | { | |
27 | } | |
28 | ||
29 | //____________________________________________________________________ | |
30 | AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter(const char* title, | |
31 | Bool_t useTuple) | |
32 | : AliFMDEnergyFitter(title), | |
33 | fSumPrimary(0), | |
34 | fSumSecondary(0), | |
35 | fIp(), | |
36 | fNTrack(0), | |
37 | fNPrimary(0), | |
38 | fTuple(0) | |
39 | { | |
40 | // Some defaults | |
41 | fUseIncreasingBins = true; | |
42 | fDoFits = true; | |
43 | fDoMakeObject = false; | |
44 | fResidualMethod = kNoResiduals; | |
45 | ||
46 | if (!useTuple) return; | |
47 | ||
48 | fTuple = new TNtuple("hits", "Hits", | |
49 | "primary:eta:phi:edep:length:" | |
50 | "detector:ring:sector:strip:ipz"); | |
51 | ||
52 | } | |
53 | //____________________________________________________________________ | |
54 | AliFMDMCHitEnergyFitter::~AliFMDMCHitEnergyFitter() | |
55 | {} | |
56 | ||
57 | //____________________________________________________________________ | |
58 | void | |
59 | AliFMDMCHitEnergyFitter::CreateOutputObjects(TList* dir) | |
60 | { | |
61 | DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir); | |
62 | AliFMDEnergyFitter::CreateOutputObjects(dir); | |
63 | // TList* d = static_cast<TList*>(l->FindObject(GetName())); | |
64 | } | |
65 | //____________________________________________________________________ | |
66 | Bool_t | |
67 | AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput) | |
68 | { | |
69 | DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput); | |
70 | fSumPrimary.Reset(0); | |
71 | fSumSecondary.Reset(0); | |
72 | ||
73 | AliMCEvent& mc = const_cast<AliMCEvent&>(mcInput); | |
74 | AliHeader* header = mc.Header(); | |
75 | AliStack* stack = mc.Stack(); | |
76 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
77 | ||
78 | genHeader->PrimaryVertex(fIp); | |
79 | fNTrack = stack->GetNtrack(); | |
80 | fNPrimary = stack->GetNprimary(); | |
81 | ||
82 | return true; | |
83 | } | |
84 | ||
85 | //____________________________________________________________________ | |
86 | Bool_t | |
87 | AliFMDMCHitEnergyFitter::Event(const AliESDEvent& esdInput, | |
88 | const AliMCEvent& mcInput, | |
89 | AliMCAuxHandler& handler) | |
90 | { | |
91 | DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)", | |
92 | &esdInput, &mcInput, &handler); | |
93 | PreEvent(mcInput); | |
94 | ||
95 | Int_t nEntries = handler.GetNEntry(); | |
96 | DMSG(fDebug,5, "We got a total of %d particles", nEntries); | |
97 | ||
98 | for (Int_t i = 0; i < nEntries; i++) { | |
99 | TClonesArray* hits = handler.GetEntryArray(i); | |
100 | if (!hits) continue; | |
101 | ||
102 | AccumulateHits(mcInput, *hits); | |
103 | } | |
104 | return PostEvent(esdInput); | |
105 | } | |
106 | ||
107 | //____________________________________________________________________ | |
108 | Bool_t | |
109 | AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput, | |
110 | const TClonesArray& hits) | |
111 | { | |
112 | DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)", | |
113 | &mcInput, &hits, hits.GetEntries()); | |
114 | ||
115 | AliStack* stack = const_cast<AliMCEvent&>(mcInput).Stack(); | |
116 | Float_t cache[10]; | |
117 | Int_t nHit = hits.GetEntries(); | |
118 | Int_t oldTr = -1; | |
119 | UShort_t oldD = 0; | |
120 | Char_t oldR = '\0'; | |
121 | for (Int_t j = 0; j < nHit; j++) { | |
122 | ||
123 | // Check the hit | |
124 | AliFMDHit* hit = static_cast<AliFMDHit*>(hits.At(j)); | |
125 | if (!hit) continue; | |
126 | ||
127 | // Zero fill array | |
128 | if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0; | |
129 | ||
130 | // Get the track number | |
131 | Int_t iTr = hit->Track(); | |
132 | if (iTr < 0 || iTr >= fNTrack) | |
133 | AliWarningF("Track # %d seems out of bounds [0,%d]", | |
134 | iTr, fNTrack-1); | |
135 | ||
136 | // Get the track | |
137 | AliMCParticle* particle = | |
138 | static_cast<AliMCParticle*>(mcInput.GetTrack(iTr)); | |
139 | ||
140 | // Particle type 0: unknown, 1: primary, 2: secondary | |
141 | // Int_t flag = 0; - not used at the moment | |
142 | ||
143 | // Check if this charged and a primary | |
144 | if (particle) { | |
145 | // Only charged tracks | |
146 | if (!particle->Charge() != 0) continue; | |
147 | } | |
148 | Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary; | |
149 | Double_t eloss = hit->Edep(); | |
150 | Double_t length = hit->Length(); | |
151 | // Double_t eta = particle->Eta(); | |
152 | Double_t x = hit->X() - fIp[0]; | |
153 | Double_t y = hit->Y() - fIp[1]; | |
154 | Double_t z = hit->Z() - fIp[2]; | |
155 | Double_t r = TMath::Sqrt(x*x + y*y); | |
156 | Double_t phi = TMath::ATan2(y, x); | |
157 | Double_t theta = TMath::ATan2(r, z); | |
158 | Double_t eta = -TMath::Log(TMath::Tan(theta/2)); | |
159 | UShort_t det = hit->Detector(); | |
160 | Char_t rng = hit->Ring(); | |
161 | UShort_t sec = hit->Sector(); | |
162 | UShort_t str = hit->Strip(); | |
163 | ||
164 | Bool_t isDistinct = (iTr != oldTr) || (det != oldD) || (rng != oldR); | |
165 | if (isDistinct) { | |
166 | RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng)); | |
167 | h->fKind->Fill(eta, isPrimary ? .5 : 1.5); | |
168 | } | |
169 | oldTr = iTr; | |
170 | oldD = det; | |
171 | oldR = rng; | |
172 | ||
173 | AliFMDFloatMap& sum = (isPrimary ? fSumPrimary : fSumSecondary); | |
174 | // Float_t old = sum(det, rng, sec, str); | |
175 | sum(det, rng, sec, str) += eloss; | |
176 | #if 0 | |
177 | Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f", | |
178 | det, rng, sec, str, (isPrimary ? "1st" : "2nd"), | |
179 | eloss, old, sum(det, rng, sec, str)); | |
180 | #endif | |
181 | if (!fTuple) continue; | |
182 | ||
183 | // Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz | |
184 | cache[0] = isPrimary ? 1 : 0; | |
185 | cache[1] = eta; | |
186 | cache[2] = phi; | |
187 | cache[3] = eloss; | |
188 | cache[4] = length; | |
189 | cache[5] = hit->Detector(); | |
190 | cache[6] = Int_t(hit->Ring()); | |
191 | cache[7] = hit->Sector(); | |
192 | cache[8] = hit->Strip(); | |
193 | cache[9] = fIp[2]; | |
194 | fTuple->Fill(cache); | |
195 | } | |
196 | return true; | |
197 | } | |
198 | ||
199 | //____________________________________________________________________ | |
200 | Bool_t AliFMDMCHitEnergyFitter::PostEvent(const AliESDEvent& input) | |
201 | { | |
202 | DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input); | |
203 | ||
204 | AliESDFMD* esdFMD = input.GetFMDData(); | |
205 | for (UShort_t d = 1; d <= 3; d++) { | |
206 | UShort_t nQ = (d == 1 ? 1 : 2); | |
207 | for (UShort_t q = 0; q < nQ; q++) { | |
208 | UShort_t nS = (q == 0 ? 20 : 40); | |
209 | UShort_t nT = (q == 0 ? 512 : 256); | |
210 | Char_t r = (q == 0 ? 'I' : 'O'); | |
211 | RingHistos* h = | |
212 | static_cast<AliFMDMCHitEnergyFitter::RingHistos*>(GetRingHistos(d,r)); | |
213 | if (!h) continue; | |
214 | ||
215 | for (UShort_t s = 0; s < nS; s++) { | |
216 | for (UShort_t t = 0; t < nT; t++) { | |
217 | Float_t totPrim = fSumPrimary(d, r, s, t); | |
218 | Float_t totSec = fSumSecondary(d, r, s, t); | |
219 | Float_t totAll = totPrim + totSec; | |
220 | Double_t esdEta = esdFMD->Eta(d, r, s, t); | |
221 | if (totAll > 0) h->FillMC(0, esdEta, totAll); | |
222 | if (totPrim > 0) h->FillMC(1, esdEta, totPrim); | |
223 | if (totSec > 0) h->FillMC(2, esdEta, totSec); | |
224 | } | |
225 | } | |
226 | } | |
227 | } | |
228 | return true; | |
229 | } | |
230 | ||
231 | //____________________________________________________________________ | |
232 | AliFMDEnergyFitter::RingHistos* | |
233 | AliFMDMCHitEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const | |
234 | { | |
235 | // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r); | |
236 | return new AliFMDMCHitEnergyFitter::RingHistos(d,r); | |
237 | } | |
238 | ||
239 | //==================================================================== | |
240 | AliFMDMCHitEnergyFitter::RingHistos::RingHistos() | |
241 | : AliFMDEnergyFitter::RingHistos(), | |
242 | fPrimary(0), | |
243 | fSecondary(0), | |
244 | fKind(0) | |
245 | {} | |
246 | ||
247 | //____________________________________________________________________ | |
248 | AliFMDMCHitEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r) | |
249 | : AliFMDEnergyFitter::RingHistos(d,r), | |
250 | fPrimary(0), | |
251 | fSecondary(0), | |
252 | fKind(0) | |
253 | {} | |
254 | //____________________________________________________________________ | |
255 | TArrayD | |
256 | AliFMDMCHitEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t, | |
257 | Double_t, | |
258 | Double_t) const | |
259 | { | |
260 | // Make an increasing axis for ELoss distributions. | |
261 | // | |
262 | // We use the service function of AliFMDEncodedEdx to do this | |
263 | TArrayD ret; | |
264 | const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis(); | |
265 | s.FillBinArray(ret); | |
266 | ||
267 | return ret; | |
268 | } | |
269 | //____________________________________________________________________ | |
270 | void | |
271 | AliFMDMCHitEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis, | |
272 | const TAxis& /*cAxis*/, | |
273 | Double_t maxDE, | |
274 | Int_t nDEbins, | |
275 | Bool_t useIncrBin) | |
276 | { | |
277 | // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins, | |
278 | // useIncrBin); | |
279 | ||
280 | fHist = Make("eloss", "#sum#Delta_{true} of all", | |
281 | eAxis, maxDE, nDEbins, useIncrBin); | |
282 | fPrimary = Make("primary", "#sum#Delta_{true} of primaries", | |
283 | eAxis, maxDE, nDEbins, useIncrBin); | |
284 | fSecondary = Make("secondary","#sum#Delta_{true} of secondaries", | |
285 | eAxis, maxDE, nDEbins, useIncrBin); | |
286 | fHist->SetXTitle("#Delta_{true}"); | |
287 | fSecondary->SetXTitle("#Delta_{true}"); | |
288 | fPrimary->SetXTitle("#Delta_{true}"); | |
289 | ||
290 | fKind = 0; | |
291 | if (eAxis.GetXbins()->GetArray()) | |
292 | fKind = new TH2I("kind", "Particle types", | |
293 | eAxis.GetNbins(), eAxis.GetXbins()->GetArray(), | |
294 | 2, 0, 2); | |
295 | else | |
296 | fKind = new TH2I("kind", "Particle types", | |
297 | eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), | |
298 | 2, 0, 2); | |
299 | fKind->SetXTitle("#eta"); | |
300 | fKind->GetYaxis()->SetBinLabel(1, "Primary"); | |
301 | fKind->GetYaxis()->SetBinLabel(2, "Secondary"); | |
302 | fKind->SetDirectory(0); | |
303 | fKind->Sumw2(); | |
304 | ||
305 | fList->Add(fHist); | |
306 | fList->Add(fPrimary); | |
307 | fList->Add(fSecondary); | |
308 | fList->Add(fKind); | |
309 | } | |
310 | ||
311 | //____________________________________________________________________ | |
312 | void | |
313 | AliFMDMCHitEnergyFitter::RingHistos::FillMC(UShort_t flag, Double_t eta, | |
314 | Double_t mult) | |
315 | { | |
316 | switch (flag) { | |
317 | // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break; | |
318 | case 0: fHist->Fill(eta, mult); break; | |
319 | case 1: fPrimary->Fill(eta, mult); break; | |
320 | case 2: fSecondary->Fill(eta, mult); break; | |
321 | } | |
322 | } | |
323 | ||
324 | //____________________________________________________________________ | |
325 | TObjArray* | |
326 | AliFMDMCHitEnergyFitter::RingHistos::Fit(TList* dir, | |
327 | Double_t lowCut, | |
328 | UShort_t nParticles, | |
329 | UShort_t minEntries, | |
330 | UShort_t minusBins, | |
331 | Double_t relErrorCut, | |
332 | Double_t chi2nuCut, | |
333 | Double_t minWeight, | |
334 | Double_t regCut, | |
335 | EResidualMethod residuals) const | |
336 | { | |
337 | TObjArray* all = FitSlices(dir, "eloss", lowCut, nParticles, minEntries, | |
338 | minusBins, relErrorCut, chi2nuCut, minWeight, | |
339 | regCut, residuals, false, 0); | |
340 | TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries, | |
341 | minusBins, relErrorCut, chi2nuCut, minWeight, | |
342 | regCut, residuals, false, 0); | |
343 | TObjArray* sec = FitSlices(dir, "secondary", lowCut, nParticles, | |
344 | minEntries, minusBins, relErrorCut, chi2nuCut, | |
345 | minWeight, regCut, residuals, false, 0); | |
346 | if (!all || !prim || !sec) { | |
347 | AliWarningF("Failed to get results for all(%p), primary(%p), and/or " | |
348 | "secondary(%p)", all, prim, sec); | |
349 | return 0; | |
350 | } | |
351 | // Already added to the sub-folders | |
352 | // dir->Add(all); | |
353 | // dir->Add(prim); | |
354 | // dir->Add(sec); | |
355 | ||
356 | Int_t nPrim = prim->GetEntriesFast(); | |
357 | TObjArray* out = new TObjArray; | |
358 | for (Int_t i = 0; i < nPrim; i++) { | |
359 | TH1* h = static_cast<TH1*>(prim->At(i)); | |
360 | if (!h) continue; | |
361 | ||
362 | TAxis* xAxis = h->GetXaxis(); | |
363 | TH2* hh = 0; | |
364 | if (xAxis->GetXbins()->GetArray()) | |
365 | hh = new TH2D(h->GetName(), h->GetTitle(), | |
366 | xAxis->GetNbins(), xAxis->GetXbins()->GetArray(), | |
367 | 3, 0, 3); | |
368 | else | |
369 | hh = new TH2D(h->GetName(), h->GetTitle(), | |
370 | xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(), | |
371 | 3, 0, 3); | |
372 | hh->GetYaxis()->SetBinLabel(1, "Primary"); | |
373 | hh->GetYaxis()->SetBinLabel(2, "Secondary"); | |
374 | hh->GetYaxis()->SetBinLabel(3, "All"); | |
375 | hh->GetXaxis()->SetTitle("#eta"); | |
376 | out->Add(hh); | |
377 | } | |
378 | for (Int_t i = 0; i < nPrim; i++) { | |
379 | TH2* h = static_cast<TH2*>(out->At(i)); | |
380 | if (!h) continue; | |
381 | ||
382 | TH1* hp = static_cast<TH1*>(prim->At(i)); | |
383 | TH1* hs = static_cast<TH1*>(sec->At(i)); | |
384 | TH1* ha = static_cast<TH1*>(all->At(i)); | |
385 | TH1* hh[] = { hp, hs, ha }; | |
386 | for (Int_t j = 0; j < 3; j++) { | |
387 | TH1* ph = hh[j]; | |
388 | if (!ph) continue; | |
389 | ||
390 | for (Int_t k = 1; k <= ph->GetNbinsX(); k++) { | |
391 | Double_t c = ph->GetBinContent(k); | |
392 | Double_t e = ph->GetBinError(k); | |
393 | h->SetBinContent(k, j+1, c); | |
394 | h->SetBinError(k, j+1, e); | |
395 | } | |
396 | } | |
397 | } | |
398 | TList* l = GetOutputList(dir); | |
399 | if (!l) return 0; | |
400 | ||
401 | out->SetName("compartive"); | |
402 | out->SetOwner(); | |
403 | l->Add(out); | |
404 | return out; | |
405 | } |