]>
Commit | Line | Data |
---|---|---|
7e4038b5 | 1 | // |
7984e5f7 | 2 | // Class to do the sharing correction. That is, a filter that merges |
3 | // adjacent strip signals presumably originating from a single particle | |
4 | // that impinges on the detector in such a way that it deposite energy | |
5 | // into two or more strips. | |
6 | // | |
7 | // Input: | |
8 | // - AliESDFMD object - from reconstruction | |
9 | // | |
10 | // Output: | |
11 | // - AliESDFMD object - copy of input, but with signals merged | |
12 | // | |
13 | // Corrections used: | |
14 | // - AliFMDCorrELossFit | |
15 | // | |
16 | // Histograms: | |
17 | // - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of | |
18 | // signals before and after the filter. | |
19 | // - For each ring (see above), an array of distributions of number of | |
20 | // hit strips for each vertex bin (if enabled - see Init method) | |
21 | // | |
22 | // | |
7e4038b5 | 23 | // |
24 | #include "AliFMDSharingFilter.h" | |
25 | #include <AliESDFMD.h> | |
26 | #include <TAxis.h> | |
27 | #include <TList.h> | |
28 | #include <TH1.h> | |
29 | #include <TMath.h> | |
0bd4b00f | 30 | #include "AliForwardCorrectionManager.h" |
fb3430ac | 31 | #include "AliFMDCorrELossFit.h" |
7e4038b5 | 32 | #include <AliLog.h> |
0bd4b00f | 33 | #include <TROOT.h> |
5bb5d1f6 | 34 | #include <THStack.h> |
0bd4b00f | 35 | #include <iostream> |
36 | #include <iomanip> | |
7e4038b5 | 37 | |
38 | ClassImp(AliFMDSharingFilter) | |
39 | #if 0 | |
40 | ; // This is for Emacs | |
41 | #endif | |
42 | ||
5bb5d1f6 | 43 | #define DBG(L,M) \ |
44 | do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false) | |
45 | #define DBGL(L,M) \ | |
46 | do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false) | |
47 | ||
48 | ||
7e4038b5 | 49 | |
50 | //____________________________________________________________________ | |
51 | AliFMDSharingFilter::AliFMDSharingFilter() | |
52 | : TNamed(), | |
53 | fRingHistos(), | |
0bd4b00f | 54 | fLowCut(0.), |
ea3e5d95 | 55 | fCorrectAngles(kFALSE), |
0bd4b00f | 56 | fNXi(1), |
5bb5d1f6 | 57 | fIncludeSigma(true), |
58 | fSummed(0), | |
59 | fHighCuts(0), | |
60 | fOper(0), | |
ea3e5d95 | 61 | fDebug(0) |
7984e5f7 | 62 | { |
63 | // | |
64 | // Default Constructor - do not use | |
65 | // | |
66 | } | |
7e4038b5 | 67 | |
68 | //____________________________________________________________________ | |
69 | AliFMDSharingFilter::AliFMDSharingFilter(const char* title) | |
70 | : TNamed("fmdSharingFilter", title), | |
71 | fRingHistos(), | |
0bd4b00f | 72 | fLowCut(0.), |
ea3e5d95 | 73 | fCorrectAngles(kFALSE), |
0bd4b00f | 74 | fNXi(1), |
5bb5d1f6 | 75 | fIncludeSigma(true), |
76 | fSummed(0), | |
77 | fHighCuts(0), | |
78 | fOper(0), | |
ea3e5d95 | 79 | fDebug(0) |
7e4038b5 | 80 | { |
7984e5f7 | 81 | // |
82 | // Constructor | |
83 | // | |
84 | // Parameters: | |
85 | // title Title of object - not significant | |
86 | // | |
7e4038b5 | 87 | fRingHistos.Add(new RingHistos(1, 'I')); |
88 | fRingHistos.Add(new RingHistos(2, 'I')); | |
89 | fRingHistos.Add(new RingHistos(2, 'O')); | |
90 | fRingHistos.Add(new RingHistos(3, 'I')); | |
91 | fRingHistos.Add(new RingHistos(3, 'O')); | |
7e4038b5 | 92 | } |
93 | ||
94 | //____________________________________________________________________ | |
95 | AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o) | |
96 | : TNamed(o), | |
97 | fRingHistos(), | |
98 | fLowCut(o.fLowCut), | |
ea3e5d95 | 99 | fCorrectAngles(o.fCorrectAngles), |
0bd4b00f | 100 | fNXi(o.fNXi), |
5bb5d1f6 | 101 | fIncludeSigma(o.fIncludeSigma), |
102 | fSummed(o.fSummed), | |
103 | fHighCuts(o.fHighCuts), | |
104 | fOper(o.fOper), | |
ea3e5d95 | 105 | fDebug(o.fDebug) |
7e4038b5 | 106 | { |
7984e5f7 | 107 | // |
108 | // Copy constructor | |
109 | // | |
110 | // Parameters: | |
111 | // o Object to copy from | |
112 | // | |
7e4038b5 | 113 | TIter next(&o.fRingHistos); |
114 | TObject* obj = 0; | |
115 | while ((obj = next())) fRingHistos.Add(obj); | |
116 | } | |
117 | ||
118 | //____________________________________________________________________ | |
119 | AliFMDSharingFilter::~AliFMDSharingFilter() | |
120 | { | |
7984e5f7 | 121 | // |
122 | // Destructor | |
123 | // | |
7e4038b5 | 124 | fRingHistos.Delete(); |
125 | } | |
126 | ||
127 | //____________________________________________________________________ | |
128 | AliFMDSharingFilter& | |
129 | AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o) | |
130 | { | |
7984e5f7 | 131 | // |
132 | // Assignment operator | |
133 | // | |
134 | // Parameters: | |
135 | // o Object to assign from | |
136 | // | |
137 | // Return: | |
138 | // Reference to this | |
139 | // | |
ea3e5d95 | 140 | TNamed::operator=(o); |
7e4038b5 | 141 | |
142 | fLowCut = o.fLowCut; | |
143 | fCorrectAngles = o.fCorrectAngles; | |
0bd4b00f | 144 | fNXi = o.fNXi; |
ea3e5d95 | 145 | fDebug = o.fDebug; |
5bb5d1f6 | 146 | fOper = o.fOper; |
147 | fSummed = o.fSummed; | |
148 | fHighCuts = o.fHighCuts; | |
149 | fIncludeSigma = o.fIncludeSigma; | |
7e4038b5 | 150 | |
151 | fRingHistos.Delete(); | |
152 | TIter next(&o.fRingHistos); | |
153 | TObject* obj = 0; | |
154 | while ((obj = next())) fRingHistos.Add(obj); | |
155 | ||
156 | return *this; | |
157 | } | |
158 | ||
159 | //____________________________________________________________________ | |
160 | AliFMDSharingFilter::RingHistos* | |
161 | AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const | |
162 | { | |
7984e5f7 | 163 | // |
164 | // Get the ring histogram container | |
165 | // | |
166 | // Parameters: | |
167 | // d Detector | |
168 | // r Ring | |
169 | // | |
170 | // Return: | |
171 | // Ring histogram container | |
172 | // | |
7e4038b5 | 173 | Int_t idx = -1; |
174 | switch (d) { | |
175 | case 1: idx = 0; break; | |
176 | case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break; | |
177 | case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break; | |
178 | } | |
179 | if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0; | |
180 | ||
181 | return static_cast<RingHistos*>(fRingHistos.At(idx)); | |
182 | } | |
183 | ||
5bb5d1f6 | 184 | //____________________________________________________________________ |
185 | void | |
186 | AliFMDSharingFilter::Init() | |
187 | { | |
188 | // Initialise | |
189 | AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); | |
190 | ||
191 | // Get the high cut. The high cut is defined as the | |
192 | // most-probably-value peak found from the energy distributions, minus | |
193 | // 2 times the width of the corresponding Landau. | |
194 | AliFMDCorrELossFit* fits = fcm.GetELossFit(); | |
195 | const TAxis& eAxis = fits->GetEtaAxis(); | |
196 | ||
197 | UShort_t nEta = eAxis.GetNbins(); | |
198 | fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5); | |
199 | fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i"); | |
200 | fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i"); | |
201 | fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o"); | |
202 | fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i"); | |
203 | fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o"); | |
204 | ||
205 | UShort_t ybin = 0; | |
206 | for (UShort_t d = 1; d <= 3; d++) { | |
207 | UShort_t nr = (d == 1 ? 1 : 2); | |
208 | for (UShort_t q = 0; q < nr; q++) { | |
209 | Char_t r = (q == 0 ? 'I' : 'O'); | |
210 | ybin++; | |
211 | for (UShort_t e = 1; e <= nEta; e++) { | |
212 | Double_t eta = eAxis.GetBinCenter(e); | |
213 | Double_t cut = GetHighCut(d, r, eta, false); | |
214 | if (cut <= 0) continue; | |
215 | fHighCuts->SetBinContent(e, ybin, cut); | |
216 | } | |
217 | } | |
218 | } | |
219 | } | |
220 | ||
7e4038b5 | 221 | //____________________________________________________________________ |
222 | Bool_t | |
223 | AliFMDSharingFilter::Filter(const AliESDFMD& input, | |
224 | Bool_t lowFlux, | |
ea3e5d95 | 225 | AliESDFMD& output) |
7e4038b5 | 226 | { |
7984e5f7 | 227 | // |
228 | // Filter the input AliESDFMD object | |
229 | // | |
230 | // Parameters: | |
231 | // input Input | |
232 | // lowFlux If this is a low-flux event | |
233 | // output Output AliESDFMD object | |
234 | // | |
235 | // Return: | |
236 | // True on success, false otherwise | |
237 | // | |
7e4038b5 | 238 | output.Clear(); |
239 | TIter next(&fRingHistos); | |
0bd4b00f | 240 | RingHistos* o = 0; |
7e4038b5 | 241 | while ((o = static_cast<RingHistos*>(next()))) |
242 | o->Clear(); | |
243 | ||
5bb5d1f6 | 244 | if (fOper) fOper->Reset(0); |
245 | Int_t nNone = 0; | |
246 | Int_t nCandidate = 0; | |
247 | Int_t nMerged = 0; | |
248 | Int_t nSummed = 0; | |
7e4038b5 | 249 | |
5bb5d1f6 | 250 | Status status[512]; |
251 | ||
7e4038b5 | 252 | for(UShort_t d = 1; d <= 3; d++) { |
253 | Int_t nRings = (d == 1 ? 1 : 2); | |
254 | for (UShort_t q = 0; q < nRings; q++) { | |
5bb5d1f6 | 255 | Char_t r = (q == 0 ? 'I' : 'O'); |
256 | UShort_t nsec = (q == 0 ? 20 : 40); | |
257 | UShort_t nstr = (q == 0 ? 512 : 256); | |
7e4038b5 | 258 | RingHistos* histos = GetRingHistos(d, r); |
259 | ||
260 | for(UShort_t s = 0; s < nsec; s++) { | |
5bb5d1f6 | 261 | for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate; |
262 | #ifdef USE_OLDER_MERGING | |
263 | Bool_t usedThis = kFALSE; | |
264 | Bool_t usedPrev = kFALSE; | |
265 | #endif | |
7e4038b5 | 266 | for(UShort_t t = 0; t < nstr; t++) { |
267 | output.SetMultiplicity(d,r,s,t,0.); | |
268 | Float_t mult = SignalInStrip(input,d,r,s,t); | |
5bb5d1f6 | 269 | // Get the pseudo-rapidity |
270 | Double_t eta = input.Eta(d,r,s,t); | |
271 | Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.; | |
272 | if (s == 0) output.SetEta(d,r,s,t,eta); | |
7e4038b5 | 273 | |
274 | // Keep dead-channel information. | |
275 | if(mult == AliESDFMD::kInvalidMult) | |
276 | output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult); | |
277 | ||
278 | // If no signal or dead strip, go on. | |
5bb5d1f6 | 279 | if (mult == AliESDFMD::kInvalidMult || mult == 0) { |
280 | if (mult == 0) histos->fSum->Fill(eta,phi,mult); | |
281 | status[t] = kNone; | |
282 | continue; | |
283 | } | |
7e4038b5 | 284 | |
7e4038b5 | 285 | // Fill the diagnostics histogram |
286 | histos->fBefore->Fill(mult); | |
287 | ||
288 | // Get next and previous signal - if any | |
289 | Double_t prevE = 0; | |
290 | Double_t nextE = 0; | |
5bb5d1f6 | 291 | Status prevStatus = (t == 0 ? kNone : status[t-1]); |
292 | Status thisStatus = status[t]; | |
293 | Status nextStatus = (t == nstr-1 ? kNone : status[t+1]); | |
7e4038b5 | 294 | if (t != 0) { |
295 | prevE = SignalInStrip(input,d,r,s,t-1); | |
296 | if (prevE == AliESDFMD::kInvalidMult) prevE = 0; | |
297 | } | |
298 | if (t != nstr - 1) { | |
299 | nextE = SignalInStrip(input,d,r,s,t+1); | |
300 | if (nextE == AliESDFMD::kInvalidMult) nextE = 0; | |
301 | } | |
5bb5d1f6 | 302 | if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult); |
7e4038b5 | 303 | |
5bb5d1f6 | 304 | #ifdef USE_OLDER_MERGING |
7e4038b5 | 305 | Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE, |
306 | lowFlux,d,r,s,t, | |
307 | usedPrev,usedThis); | |
5bb5d1f6 | 308 | status[t] = (usedPrev ? kMergedWithOther : kNone); |
309 | if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone); | |
310 | #else | |
311 | Double_t mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, | |
312 | eta, lowFlux, | |
313 | d, r, s, t, | |
314 | prevStatus, | |
315 | thisStatus, | |
316 | nextStatus); | |
317 | if (t != 0) status[t-1] = prevStatus; | |
318 | if (t != nstr-1) status[t+1] = nextStatus; | |
319 | status[t] = thisStatus; | |
320 | ||
321 | #endif | |
322 | // If we're processing on non-angle corrected data, we | |
323 | // should do the angle correction here | |
324 | if (!fCorrectAngles) | |
325 | mergedEnergy = AngleCorrect(mergedEnergy, eta); | |
326 | if (mergedEnergy > 0) histos->Incr(); | |
327 | ||
328 | if (t != 0) | |
329 | histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1), | |
330 | mergedEnergy); | |
331 | histos->fBeforeAfter->Fill(mult, mergedEnergy); | |
7e4038b5 | 332 | histos->fAfter->Fill(mergedEnergy); |
5bb5d1f6 | 333 | histos->fSum->Fill(eta,phi,mergedEnergy); |
7e4038b5 | 334 | |
335 | output.SetMultiplicity(d,r,s,t,mergedEnergy); | |
7e4038b5 | 336 | } // for strip |
5bb5d1f6 | 337 | for (UShort_t t = 0; t < nstr; t++) { |
338 | if (fOper) fOper->operator()(d, r, s, t) = status[t]; | |
339 | switch (status[t]) { | |
340 | case kNone: nNone++; break; | |
341 | case kCandidate: nCandidate++; break; | |
342 | case kMergedWithOther: nMerged++; break; | |
343 | case kMergedInto: nSummed++; break; | |
344 | } | |
345 | } | |
7e4038b5 | 346 | } // for sector |
347 | } // for ring | |
348 | } // for detector | |
5bb5d1f6 | 349 | fSummed->Fill(kNone, nNone); |
350 | fSummed->Fill(kCandidate, nCandidate); | |
351 | fSummed->Fill(kMergedWithOther, nMerged); | |
352 | fSummed->Fill(kMergedInto, nSummed); | |
7e4038b5 | 353 | |
5bb5d1f6 | 354 | DBGL(1, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d", |
355 | nNone, nCandidate, nMerged, nSummed)); | |
7e4038b5 | 356 | next.Reset(); |
357 | while ((o = static_cast<RingHistos*>(next()))) | |
358 | o->Finish(); | |
359 | ||
360 | return kTRUE; | |
361 | } | |
362 | ||
363 | //_____________________________________________________________________ | |
364 | Double_t | |
365 | AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input, | |
366 | UShort_t d, | |
367 | Char_t r, | |
368 | UShort_t s, | |
369 | UShort_t t) const | |
370 | { | |
7984e5f7 | 371 | // |
372 | // Get the signal in a strip | |
373 | // | |
374 | // Parameters: | |
375 | // fmd ESD object | |
376 | // d Detector | |
377 | // r Ring | |
378 | // s Sector | |
379 | // t Strip | |
380 | // | |
381 | // Return: | |
382 | // The energy signal | |
383 | // | |
7e4038b5 | 384 | Double_t mult = input.Multiplicity(d,r,s,t); |
5bb5d1f6 | 385 | // In case of |
386 | // - bad value (invalid or 0) | |
387 | // - we want angle corrected and data is | |
388 | // - we don't want angle corrected and data isn't | |
389 | // just return read value | |
390 | if (mult == AliESDFMD::kInvalidMult || | |
391 | mult == 0 || | |
392 | (fCorrectAngles && input.IsAngleCorrected()) || | |
393 | (!fCorrectAngles && !input.IsAngleCorrected())) | |
394 | return mult; | |
395 | ||
396 | // If we want angle corrected data, correct it, | |
397 | // otherwise de-correct it | |
398 | if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t)); | |
399 | else mult = DeAngleCorrect(mult, input.Eta(d,r,s,t)); | |
7e4038b5 | 400 | return mult; |
401 | } | |
0bd4b00f | 402 | //_____________________________________________________________________ |
403 | Double_t | |
5bb5d1f6 | 404 | AliFMDSharingFilter::GetLowCut(UShort_t, Char_t, Double_t) const |
0bd4b00f | 405 | { |
7984e5f7 | 406 | // |
407 | // Get the low cut. Normally, the low cut is taken to be the lower | |
408 | // value of the fit range used when generating the energy loss fits. | |
409 | // However, if fLowCut is set (using SetLowCit) to a value greater | |
410 | // than 0, then that value is used. | |
411 | // | |
0bd4b00f | 412 | if (fLowCut > 0) return fLowCut; |
413 | AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); | |
414 | AliFMDCorrELossFit* fits = fcm.GetELossFit(); | |
415 | return fits->GetLowCut(); | |
416 | } | |
7e4038b5 | 417 | |
418 | //_____________________________________________________________________ | |
419 | Double_t | |
5bb5d1f6 | 420 | AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, |
421 | Double_t eta, Bool_t errors) const | |
7e4038b5 | 422 | { |
7984e5f7 | 423 | // |
424 | // Get the high cut. The high cut is defined as the | |
425 | // most-probably-value peak found from the energy distributions, minus | |
426 | // 2 times the width of the corresponding Landau. | |
427 | // | |
0bd4b00f | 428 | AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); |
429 | ||
7e4038b5 | 430 | |
431 | // Get the high cut. The high cut is defined as the | |
432 | // most-probably-value peak found from the energy distributions, minus | |
433 | // 2 times the width of the corresponding Landau. | |
5bb5d1f6 | 434 | AliFMDCorrELossFit* fits = fcm.GetELossFit(); |
435 | ||
436 | return fits->GetLowerBound(d, r, eta, fNXi, errors, fIncludeSigma); | |
437 | } | |
438 | ||
439 | //_____________________________________________________________________ | |
440 | namespace { | |
441 | const char* status2String(AliFMDSharingFilter::Status s) | |
442 | { | |
443 | switch(s) { | |
444 | case AliFMDSharingFilter::kNone: return "none "; | |
445 | case AliFMDSharingFilter::kCandidate: return "candidate"; | |
446 | case AliFMDSharingFilter::kMergedWithOther: return "merged "; | |
447 | case AliFMDSharingFilter::kMergedInto: return "summed "; | |
448 | } | |
449 | return "unknown "; | |
0bd4b00f | 450 | } |
5bb5d1f6 | 451 | } |
452 | ||
453 | //_____________________________________________________________________ | |
454 | Double_t | |
455 | AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE, | |
456 | Double_t prevE, | |
457 | Double_t nextE, | |
458 | Double_t eta, | |
459 | Bool_t lowFlux, | |
460 | UShort_t d, | |
461 | Char_t r, | |
462 | UShort_t s, | |
463 | UShort_t t, | |
464 | Status& prevStatus, | |
465 | Status& thisStatus, | |
466 | Status& nextStatus) const | |
467 | { | |
468 | Double_t lowCut = GetLowCut(d, r, eta); | |
469 | ||
470 | DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)", | |
471 | d, r, s, t, | |
472 | thisE, status2String(thisStatus), | |
473 | prevE, status2String(prevStatus), | |
474 | nextE, status2String(nextStatus))); | |
475 | ||
476 | // If below cut, then modify zero signal and make sure the next | |
477 | // strip is considered a candidate. | |
478 | if (thisE < lowCut || thisE > 20) { | |
479 | thisStatus = kNone; | |
480 | DBGL(3,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE)); | |
481 | if (prevStatus == kCandidate) prevStatus = kNone; | |
482 | return 0; | |
483 | } | |
484 | // It this strip was merged with the previous strip, then | |
485 | // make the next strip a candidate and zero the value in this strip. | |
486 | if (thisStatus == kMergedWithOther) { | |
487 | DBGL(3,Form(" Merged with other, 0'ed")); | |
488 | return 0; | |
0bd4b00f | 489 | } |
490 | ||
5bb5d1f6 | 491 | // Get the upper cut on the sharing |
492 | Double_t highCut = GetHighCut(d, r, eta ,false); | |
493 | if (highCut <= 0) { | |
494 | thisStatus = kNone; | |
495 | return 0; | |
496 | } | |
497 | ||
498 | // Only for low-flux events | |
499 | if (lowFlux) { | |
500 | // If this is smaller than the next, and the next is larger | |
501 | // then the cut, mark both as candidates for sharing. | |
502 | // They should be be merged on the next iteration | |
503 | if (thisE < nextE && nextE > highCut) { | |
504 | thisStatus = kCandidate; | |
505 | nextStatus = kCandidate; | |
506 | return 0; | |
507 | } | |
508 | } | |
509 | ||
510 | // Variable to sum signal in | |
511 | Double_t totalE = thisE; | |
512 | ||
513 | // If the previous signal was between the two cuts, and is still | |
514 | // marked as a candidate , then merge it in here. | |
515 | if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) { | |
516 | totalE += prevE; | |
517 | prevStatus = kMergedWithOther; | |
518 | DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f", | |
519 | lowCut, prevE, highCut, totalE)); | |
520 | } | |
521 | ||
522 | // If the next signal is between the two cuts, then merge it here | |
523 | if (nextE > lowCut && nextE < highCut) { | |
524 | totalE += nextE; | |
525 | nextStatus = kMergedWithOther; | |
526 | DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE)); | |
527 | } | |
528 | ||
529 | // If the signal is bigger than the high-cut, then return the value | |
530 | if (totalE > highCut) { | |
531 | if (totalE > thisE) | |
532 | thisStatus = kMergedInto; | |
533 | else | |
534 | thisStatus = kNone; | |
535 | DBGL(3, Form(" %9f>%f9 (was %9f) -> %9f %s", | |
536 | totalE, highCut, thisE, totalE,status2String(thisStatus))); | |
537 | return totalE; | |
538 | } | |
539 | else { | |
540 | // This is below cut, so flag next as a candidate | |
541 | DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut)); | |
542 | nextStatus = kCandidate; | |
543 | } | |
544 | ||
545 | // If the total signal is smaller than low cut then zero this and kill this | |
546 | if (totalE < lowCut) { | |
547 | DBGL(3, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE)); | |
548 | thisStatus = kNone; | |
549 | return 0; | |
0bd4b00f | 550 | } |
5bb5d1f6 | 551 | |
552 | // If total signal not above high cut or lower than low cut, | |
553 | // mark this as a candidate for merging into the next, and zero signal | |
554 | DBGL(3, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate", | |
555 | lowCut, totalE, highCut, thisE)); | |
556 | thisStatus = kCandidate; | |
557 | return 0; | |
7e4038b5 | 558 | } |
559 | ||
560 | //_____________________________________________________________________ | |
561 | Double_t | |
562 | AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult, | |
563 | Double_t eta, | |
564 | Double_t prevE, | |
565 | Double_t nextE, | |
566 | Bool_t lowFlux, | |
567 | UShort_t d, | |
568 | Char_t r, | |
569 | UShort_t /*s*/, | |
570 | UShort_t /*t*/, | |
571 | Bool_t& usedPrev, | |
fb3430ac | 572 | Bool_t& usedThis) const |
7e4038b5 | 573 | { |
7984e5f7 | 574 | // |
575 | // The actual algorithm | |
576 | // | |
577 | // Parameters: | |
578 | // mult The unfiltered signal in the strip | |
579 | // eta Psuedo rapidity | |
580 | // prevE Previous strip signal (or 0) | |
581 | // nextE Next strip signal (or 0) | |
582 | // lowFlux Whether this is a low flux event | |
583 | // d Detector | |
584 | // r Ring | |
585 | // s Sector | |
586 | // t Strip | |
587 | // usedPrev Whether the previous strip was used in sharing or not | |
588 | // usedThis Wether this strip was used in sharing or not. | |
589 | // | |
590 | // Return: | |
591 | // The filtered signal in the strip | |
592 | // | |
593 | ||
7e4038b5 | 594 | // IF the multiplicity is very large, or below the cut, reject it, and |
595 | // flag it as candidate for sharing | |
5bb5d1f6 | 596 | Double_t lowCut = GetLowCut(d,r,eta); |
0bd4b00f | 597 | if(mult > 12 || mult < lowCut) { |
7e4038b5 | 598 | usedThis = kFALSE; |
599 | usedPrev = kFALSE; | |
600 | return 0; | |
601 | } | |
602 | ||
603 | // If this was shared with the previous signal, return 0 and mark it as | |
604 | // not a candiate for sharing | |
605 | if(usedThis) { | |
606 | usedThis = kFALSE; | |
607 | usedPrev = kTRUE; | |
608 | return 0.; | |
609 | } | |
610 | ||
611 | //analyse and perform sharing on one strip | |
612 | ||
0bd4b00f | 613 | // Get the high cut |
5bb5d1f6 | 614 | Double_t highCut = GetHighCut(d, r, eta, false); |
615 | if (highCut <= 0) { | |
616 | usedThis = kFALSE; | |
617 | usedPrev = kTRUE; | |
618 | return 0; | |
619 | } | |
7e4038b5 | 620 | |
621 | // If this signal is smaller than the next, and the next signal is smaller | |
622 | // than then the high cut, and it's a low-flux event, then mark this and | |
623 | // the next signal as candidates for sharing | |
624 | // | |
625 | // This is the test if the signal is the smaller of two possibly | |
626 | // shared signals. | |
627 | if (mult < nextE && | |
628 | nextE > highCut && | |
629 | lowFlux) { | |
630 | usedThis = kFALSE; | |
631 | usedPrev = kFALSE; | |
632 | return 0; | |
633 | } | |
634 | ||
635 | Double_t totalE = mult; | |
636 | ||
637 | // If the previous signal was larger than the low cut, and | |
638 | // the previous was smaller than high cut, and the previous signal | |
639 | // isn't marked as used, then add it's energy to this signal | |
0bd4b00f | 640 | if (prevE > lowCut && |
7e4038b5 | 641 | prevE < highCut && |
642 | !usedPrev) | |
643 | totalE += prevE; | |
644 | ||
645 | // If the next signal is larger than the low cut, and | |
5bb5d1f6 | 646 | // the next signal is smaller than the high cut, then add the next signal |
7e4038b5 | 647 | // to this, and marked the next signal as used |
0bd4b00f | 648 | if(nextE > lowCut && nextE < highCut ) { |
7e4038b5 | 649 | totalE += nextE; |
650 | usedThis = kTRUE; | |
651 | } | |
652 | ||
653 | // If we're processing on non-angle corrected data, we should do the | |
654 | // angle correction here | |
5bb5d1f6 | 655 | // if (!fCorrectAngles) |
656 | // totalE = AngleCorrect(totalE, eta); | |
7e4038b5 | 657 | |
658 | // Fill post processing histogram | |
659 | // if(totalE > fLowCut) | |
660 | // GetRingHistos(d,r)->fAfter->Fill(totalE); | |
661 | ||
662 | // Return value | |
663 | Double_t mergedEnergy = 0; | |
664 | ||
665 | if(totalE > 0) { | |
666 | // If we have a signal, then this is used (sharedPrev=true) and | |
667 | // the signal is set to the result | |
668 | mergedEnergy = totalE; | |
669 | usedPrev = kTRUE; | |
670 | } | |
671 | else { | |
672 | // If we do not have a signal (too low), then this is not shared, | |
673 | // and the next strip is not shared either | |
674 | usedThis = kFALSE; | |
675 | usedPrev = kFALSE; | |
676 | } | |
677 | ||
678 | return mergedEnergy; | |
679 | } | |
680 | //____________________________________________________________________ | |
681 | Double_t | |
682 | AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const | |
683 | { | |
7984e5f7 | 684 | // |
685 | // Angle correct the signal | |
686 | // | |
687 | // Parameters: | |
688 | // mult Angle Un-corrected Signal | |
689 | // eta Pseudo-rapidity | |
690 | // | |
691 | // Return: | |
692 | // Angle corrected signal | |
693 | // | |
7e4038b5 | 694 | Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta)); |
695 | if (eta < 0) theta -= TMath::Pi(); | |
696 | return mult * TMath::Cos(theta); | |
697 | } | |
698 | //____________________________________________________________________ | |
699 | Double_t | |
700 | AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const | |
701 | { | |
7984e5f7 | 702 | // |
703 | // Angle de-correct the signal | |
704 | // | |
705 | // Parameters: | |
706 | // mult Angle corrected Signal | |
707 | // eta Pseudo-rapidity | |
708 | // | |
709 | // Return: | |
710 | // Angle un-corrected signal | |
711 | // | |
7e4038b5 | 712 | Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta)); |
713 | if (eta < 0) theta -= TMath::Pi(); | |
714 | return mult / TMath::Cos(theta); | |
715 | } | |
716 | ||
717 | //____________________________________________________________________ | |
718 | void | |
fb3430ac | 719 | AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents) |
7e4038b5 | 720 | { |
7984e5f7 | 721 | // |
722 | // Scale the histograms to the total number of events | |
723 | // | |
724 | // Parameters: | |
725 | // dir Where the output is | |
726 | // nEvents Number of events | |
727 | // | |
7e4038b5 | 728 | if (nEvents <= 0) return; |
9d99b0dd | 729 | TList* d = static_cast<TList*>(dir->FindObject(GetName())); |
730 | if (!d) return; | |
7e4038b5 | 731 | |
732 | TIter next(&fRingHistos); | |
733 | RingHistos* o = 0; | |
5bb5d1f6 | 734 | THStack* sums = new THStack("sums", "Sum of ring signals"); |
735 | while ((o = static_cast<RingHistos*>(next()))) { | |
9d99b0dd | 736 | o->ScaleHistograms(d, nEvents); |
5bb5d1f6 | 737 | TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e"); |
738 | sum->Scale(1., "width"); | |
739 | sum->SetTitle(o->GetName()); | |
740 | sum->SetDirectory(0); | |
741 | sum->SetYTitle("#sum #Delta/#Delta_{mip}"); | |
742 | sums->Add(sum); | |
743 | } | |
744 | d->Add(sums); | |
7e4038b5 | 745 | } |
746 | ||
747 | //____________________________________________________________________ | |
748 | void | |
9d99b0dd | 749 | AliFMDSharingFilter::DefineOutput(TList* dir) |
7e4038b5 | 750 | { |
7984e5f7 | 751 | // |
752 | // Define the output histograms. These are put in a sub list of the | |
753 | // passed list. The histograms are merged before the parent task calls | |
754 | // AliAnalysisTaskSE::Terminate | |
755 | // | |
756 | // Parameters: | |
757 | // dir Directory to add to | |
758 | // | |
7e4038b5 | 759 | TList* d = new TList; |
760 | d->SetName(GetName()); | |
761 | dir->Add(d); | |
5bb5d1f6 | 762 | |
763 | fSummed = new TH2I("operations", "Strip operations", | |
764 | kMergedInto, kNone-.5, kMergedInto+.5, | |
765 | 51201, -.5, 51200.5); | |
766 | fSummed->SetXTitle("Operation"); | |
767 | fSummed->SetYTitle("# of strips"); | |
768 | fSummed->SetZTitle("Events"); | |
769 | fSummed->GetXaxis()->SetBinLabel(kNone, "None"); | |
770 | fSummed->GetXaxis()->SetBinLabel(kCandidate, "Candidate"); | |
771 | fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other"); | |
772 | fSummed->GetXaxis()->SetBinLabel(kMergedInto, "Merged into"); | |
773 | fSummed->SetDirectory(0); | |
774 | d->Add(fSummed); | |
775 | ||
776 | fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1); | |
777 | fHighCuts->SetXTitle("#eta"); | |
778 | fHighCuts->SetDirectory(0); | |
779 | d->Add(fHighCuts); | |
780 | ||
781 | ||
7e4038b5 | 782 | TIter next(&fRingHistos); |
783 | RingHistos* o = 0; | |
784 | while ((o = static_cast<RingHistos*>(next()))) { | |
785 | o->Output(d); | |
786 | } | |
787 | } | |
0bd4b00f | 788 | //____________________________________________________________________ |
789 | void | |
790 | AliFMDSharingFilter::Print(Option_t* /*option*/) const | |
791 | { | |
7984e5f7 | 792 | // |
793 | // Print information | |
794 | // | |
795 | // Parameters: | |
796 | // option Not used | |
797 | // | |
0bd4b00f | 798 | char ind[gROOT->GetDirLevel()+1]; |
799 | for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; | |
800 | ind[gROOT->GetDirLevel()] = '\0'; | |
1f480471 | 801 | std::cout << ind << ClassName() << ": " << GetName() << '\n' |
5bb5d1f6 | 802 | << std::boolalpha |
803 | << ind << " Debug: " << fDebug << "\n" | |
0bd4b00f | 804 | << ind << " Low cut: " << fLowCut << '\n' |
805 | << ind << " N xi factor: " << fNXi << '\n' | |
5bb5d1f6 | 806 | << ind << " Include sigma in cut: " << fIncludeSigma << '\n' |
807 | << ind << " Use corrected angles: " << fCorrectAngles | |
808 | << std::noboolalpha << std::endl; | |
0bd4b00f | 809 | } |
7e4038b5 | 810 | |
811 | //==================================================================== | |
812 | AliFMDSharingFilter::RingHistos::RingHistos() | |
9d99b0dd | 813 | : AliForwardUtil::RingHistos(), |
7e4038b5 | 814 | fBefore(0), |
815 | fAfter(0), | |
5bb5d1f6 | 816 | fBeforeAfter(0), |
817 | fNeighborsBefore(0), | |
818 | fNeighborsAfter(0), | |
819 | fSum(0), | |
7e4038b5 | 820 | fHits(0), |
821 | fNHits(0) | |
7984e5f7 | 822 | { |
823 | // | |
824 | // Default CTOR | |
825 | // | |
826 | ||
827 | } | |
7e4038b5 | 828 | |
829 | //____________________________________________________________________ | |
830 | AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r) | |
9d99b0dd | 831 | : AliForwardUtil::RingHistos(d,r), |
7e4038b5 | 832 | fBefore(0), |
833 | fAfter(0), | |
5bb5d1f6 | 834 | fBeforeAfter(0), |
835 | fNeighborsBefore(0), | |
836 | fNeighborsAfter(0), | |
837 | fSum(0), | |
7e4038b5 | 838 | fHits(0), |
839 | fNHits(0) | |
840 | { | |
7984e5f7 | 841 | // |
842 | // Constructor | |
843 | // | |
844 | // Parameters: | |
845 | // d detector | |
846 | // r ring | |
847 | // | |
5bb5d1f6 | 848 | fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)", |
7e4038b5 | 849 | 1000, 0, 25); |
850 | fBefore->SetXTitle("#Delta E/#Delta E_{mip}"); | |
851 | fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})"); | |
5bb5d1f6 | 852 | fBefore->SetFillColor(Color()); |
7e4038b5 | 853 | fBefore->SetFillStyle(3001); |
7e4038b5 | 854 | fBefore->SetDirectory(0); |
5bb5d1f6 | 855 | |
856 | fAfter = static_cast<TH1D*>(fBefore->Clone("anaEloss")); | |
857 | fAfter->SetTitle("Energy loss in %s (sharing corrected)"); | |
858 | fAfter->SetFillColor(Color()+2); | |
7e4038b5 | 859 | fAfter->SetDirectory(0); |
860 | ||
5bb5d1f6 | 861 | Double_t max = 15; |
862 | Double_t min = -1; | |
863 | Int_t n = int((max-min) / (max / 300)); | |
864 | fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation", | |
865 | n, min, max, n, min, max); | |
866 | fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before"); | |
867 | fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after"); | |
868 | fBeforeAfter->SetZTitle("Correlation"); | |
869 | fBeforeAfter->SetDirectory(0); | |
870 | ||
871 | fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore")); | |
872 | fNeighborsBefore->SetTitle("Correlation of neighbors befire"); | |
873 | fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}"); | |
874 | fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}"); | |
875 | fNeighborsBefore->SetDirectory(0); | |
876 | ||
877 | fNeighborsAfter = | |
878 | static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter")); | |
879 | fNeighborsAfter->SetTitle("Correlation of neighbors after"); | |
880 | fNeighborsAfter->SetDirectory(0); | |
881 | ||
882 | fSum = new TH2D("summed", "Summed signal", 200, -4, 6, | |
883 | (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi()); | |
884 | fSum->SetDirectory(0); | |
885 | fSum->Sumw2(); | |
886 | fSum->SetMarkerColor(Color()); | |
887 | // fSum->SetFillColor(Color()); | |
888 | fSum->SetXTitle("#eta"); | |
889 | fSum->SetYTitle("#varphi [radians]"); | |
890 | fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) "); | |
891 | ||
892 | fHits = new TH1D("hits", "Number of hits", 200, 0, 200000); | |
7e4038b5 | 893 | fHits->SetDirectory(0); |
894 | fHits->SetXTitle("# of hits"); | |
895 | fHits->SetFillColor(kGreen+1); | |
896 | } | |
897 | //____________________________________________________________________ | |
898 | AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o) | |
9d99b0dd | 899 | : AliForwardUtil::RingHistos(o), |
7e4038b5 | 900 | fBefore(o.fBefore), |
901 | fAfter(o.fAfter), | |
5bb5d1f6 | 902 | fBeforeAfter(o.fBeforeAfter), |
903 | fNeighborsBefore(o.fNeighborsBefore), | |
904 | fNeighborsAfter(o.fNeighborsAfter), | |
905 | fSum(o.fSum), | |
7e4038b5 | 906 | fHits(o.fHits), |
907 | fNHits(o.fNHits) | |
7984e5f7 | 908 | { |
909 | // | |
910 | // Copy constructor | |
911 | // | |
912 | // Parameters: | |
913 | // o Object to copy from | |
914 | // | |
915 | } | |
7e4038b5 | 916 | |
917 | //____________________________________________________________________ | |
918 | AliFMDSharingFilter::RingHistos& | |
919 | AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o) | |
920 | { | |
7984e5f7 | 921 | // |
922 | // Assignment operator | |
923 | // | |
924 | // Parameters: | |
925 | // o Object to assign from | |
926 | // | |
927 | // Return: | |
928 | // Reference to this | |
929 | // | |
9d99b0dd | 930 | AliForwardUtil::RingHistos::operator=(o); |
7e4038b5 | 931 | fDet = o.fDet; |
932 | fRing = o.fRing; | |
933 | ||
934 | if (fBefore) delete fBefore; | |
935 | if (fAfter) delete fAfter; | |
936 | if (fHits) delete fHits; | |
937 | ||
5bb5d1f6 | 938 | fBefore = static_cast<TH1D*>(o.fBefore->Clone()); |
939 | fAfter = static_cast<TH1D*>(o.fAfter->Clone()); | |
940 | fBeforeAfter = static_cast<TH2D*>(o.fBeforeAfter->Clone()); | |
941 | fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone()); | |
942 | fNeighborsAfter = static_cast<TH2D*>(o.fNeighborsAfter->Clone()); | |
943 | fHits = static_cast<TH1D*>(o.fHits->Clone()); | |
944 | fSum = static_cast<TH2D*>(o.fSum->Clone()); | |
945 | ||
7e4038b5 | 946 | return *this; |
947 | } | |
948 | //____________________________________________________________________ | |
949 | AliFMDSharingFilter::RingHistos::~RingHistos() | |
950 | { | |
7984e5f7 | 951 | // |
952 | // Destructor | |
953 | // | |
7e4038b5 | 954 | } |
955 | //____________________________________________________________________ | |
956 | void | |
957 | AliFMDSharingFilter::RingHistos::Finish() | |
958 | { | |
7984e5f7 | 959 | // |
960 | // Finish off | |
961 | // | |
962 | // | |
7e4038b5 | 963 | fHits->Fill(fNHits); |
964 | } | |
965 | ||
9d99b0dd | 966 | //____________________________________________________________________ |
967 | void | |
fb3430ac | 968 | AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir, |
969 | Int_t nEvents) | |
9d99b0dd | 970 | { |
7984e5f7 | 971 | // |
972 | // Scale the histograms to the total number of events | |
973 | // | |
974 | // Parameters: | |
975 | // nEvents Number of events | |
976 | // dir Where the output is | |
977 | // | |
9d99b0dd | 978 | TList* l = GetOutputList(dir); |
979 | if (!l) return; | |
980 | ||
5bb5d1f6 | 981 | TH1D* before = static_cast<TH1D*>(l->FindObject("esdELoss")); |
982 | TH1D* after = static_cast<TH1D*>(l->FindObject("anaELoss")); | |
983 | TH2D* summed = static_cast<TH2D*>(l->FindObject("summed")); | |
9d99b0dd | 984 | if (before) before->Scale(1./nEvents); |
985 | if (after) after->Scale(1./nEvents); | |
5bb5d1f6 | 986 | if (summed) summed->Scale(1./nEvents); |
987 | ||
9d99b0dd | 988 | } |
989 | ||
7e4038b5 | 990 | //____________________________________________________________________ |
991 | void | |
992 | AliFMDSharingFilter::RingHistos::Output(TList* dir) | |
993 | { | |
7984e5f7 | 994 | // |
995 | // Make output | |
996 | // | |
997 | // Parameters: | |
998 | // dir where to store | |
999 | // | |
9d99b0dd | 1000 | TList* d = DefineOutputList(dir); |
5bb5d1f6 | 1001 | |
7e4038b5 | 1002 | d->Add(fBefore); |
1003 | d->Add(fAfter); | |
5bb5d1f6 | 1004 | d->Add(fBeforeAfter); |
1005 | d->Add(fNeighborsBefore); | |
1006 | d->Add(fNeighborsAfter); | |
7e4038b5 | 1007 | d->Add(fHits); |
5bb5d1f6 | 1008 | d->Add(fSum); |
1009 | ||
7e4038b5 | 1010 | dir->Add(d); |
1011 | } | |
1012 | ||
1013 | //____________________________________________________________________ | |
1014 | // | |
1015 | // EOF | |
1016 | // |