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