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.
8 // - AliESDFMD object - from reconstruction
11 // - AliESDFMD object - copy of input, but with signals merged
14 // - AliFMDCorrELossFit
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)
24 #include "AliFMDSharingFilter.h"
25 #include <AliESDFMD.h>
30 #include "AliForwardCorrectionManager.h"
31 #include "AliFMDCorrELossFit.h"
35 #include <TParameter.h>
39 ClassImp(AliFMDSharingFilter)
41 ; // This is for Emacs
45 do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false)
47 do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false)
51 //____________________________________________________________________
52 AliFMDSharingFilter::AliFMDSharingFilter()
55 fCorrectAngles(kFALSE),
61 fZeroSharedHitsBelowThreshold(false),
64 fUseSimpleMerging(false),
65 fThreeStripSharing(true),
66 fRecalculateEta(false)
69 // Default Constructor - do not use
71 DGUARD(fDebug,0, "Default CTOR for AliFMDSharingFilter");
74 //____________________________________________________________________
75 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
76 : TNamed("fmdSharingFilter", title),
78 fCorrectAngles(kFALSE),
84 fZeroSharedHitsBelowThreshold(false),
87 fUseSimpleMerging(false),
88 fThreeStripSharing(true),
89 fRecalculateEta(false)
95 // title Title of object - not significant
97 DGUARD(fDebug,0, "Named CTOR for AliFMDSharingFilter: %s", title);
98 fRingHistos.SetName(GetName());
99 fRingHistos.SetOwner();
101 fRingHistos.Add(new RingHistos(1, 'I'));
102 fRingHistos.Add(new RingHistos(2, 'I'));
103 fRingHistos.Add(new RingHistos(2, 'O'));
104 fRingHistos.Add(new RingHistos(3, 'I'));
105 fRingHistos.Add(new RingHistos(3, 'O'));
108 fHCuts.SetIncludeSigma(1);
109 fLCuts.SetMultCuts(.15);
112 //____________________________________________________________________
113 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
116 fCorrectAngles(o.fCorrectAngles),
118 fHighCuts(o.fHighCuts),
119 fLowCuts(o.fLowCuts),
122 fZeroSharedHitsBelowThreshold(o.fZeroSharedHitsBelowThreshold),
125 fUseSimpleMerging(o.fUseSimpleMerging),
126 fThreeStripSharing(o.fThreeStripSharing),
127 fRecalculateEta(o.fRecalculateEta)
133 // o Object to copy from
135 DGUARD(fDebug,0, "Copy CTOR for AliFMDSharingFilter");
136 TIter next(&o.fRingHistos);
138 while ((obj = next())) fRingHistos.Add(obj);
141 //____________________________________________________________________
142 AliFMDSharingFilter::~AliFMDSharingFilter()
147 DGUARD(fDebug,3, "DTOR for AliFMDSharingFilter");
148 fRingHistos.Delete();
151 //____________________________________________________________________
153 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
156 // Assignment operator
159 // o Object to assign from
164 DGUARD(fDebug,3, "Assigment for AliFMDSharingFilter");
165 if (&o == this) return *this;
166 TNamed::operator=(o);
168 fCorrectAngles = o.fCorrectAngles;
172 fHighCuts = o.fHighCuts;
173 fLowCuts = o.fLowCuts;
174 fZeroSharedHitsBelowThreshold = o.fZeroSharedHitsBelowThreshold;
177 fUseSimpleMerging = o.fUseSimpleMerging;
178 fThreeStripSharing = o.fThreeStripSharing;
179 fRecalculateEta = o.fRecalculateEta;
181 fRingHistos.Delete();
182 TIter next(&o.fRingHistos);
184 while ((obj = next())) fRingHistos.Add(obj);
189 //____________________________________________________________________
190 AliFMDSharingFilter::RingHistos*
191 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
194 // Get the ring histogram container
201 // Ring histogram container
205 case 1: idx = 0; break;
206 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
207 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
209 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
211 return static_cast<RingHistos*>(fRingHistos.At(idx));
214 //____________________________________________________________________
216 AliFMDSharingFilter::Init(const TAxis& axis)
219 DGUARD(fDebug,1, "Initialize for AliFMDSharingFilter");
220 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
221 AliFMDCorrELossFit* fits = fcm.GetELossFit();
223 // Get the high cut. The high cut is defined as the
224 // most-probably-value peak found from the energy distributions, minus
225 // 2 times the width of the corresponding Landau.
227 TAxis eAxis(axis.GetNbins(),
231 eAxis.Set(fits->GetEtaAxis().GetNbins(),
232 fits->GetEtaAxis().GetXmin(),
233 fits->GetEtaAxis().GetXmax());
235 UShort_t nEta = eAxis.GetNbins();
237 fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
238 fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
239 fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
240 fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
241 fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
242 fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
244 fLowCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
245 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
246 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
247 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
248 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
249 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
252 for (UShort_t d = 1; d <= 3; d++) {
253 UShort_t nr = (d == 1 ? 1 : 2);
254 for (UShort_t q = 0; q < nr; q++) {
255 Char_t r = (q == 0 ? 'I' : 'O');
257 for (UShort_t e = 1; e <= nEta; e++) {
258 Double_t eta = eAxis.GetBinCenter(e);
260 if (fDebug > 3) fHCuts.Print();
262 Double_t hcut = GetHighCut(d, r, eta, false);
263 Double_t lcut = GetLowCut(d, r, eta);
265 if (hcut > 0) fHighCuts->SetBinContent(e, ybin, hcut);
266 if (lcut > 0) fLowCuts ->SetBinContent(e, ybin, lcut);
272 //____________________________________________________________________
274 AliFMDSharingFilter::Filter(const AliESDFMD& input,
280 // Filter the input AliESDFMD object
284 // lowFlux If this is a low-flux event
285 // output Output AliESDFMD object
288 // True on success, false otherwise
290 DGUARD(fDebug,1, "Filter event in AliFMDSharingFilter");
292 TIter next(&fRingHistos);
294 while ((o = static_cast<RingHistos*>(next())))
297 if (fOper) fOper->Reset(0);
299 Int_t nCandidate = 0;
305 for(UShort_t d = 1; d <= 3; d++) {
306 Int_t nRings = (d == 1 ? 1 : 2);
307 for (UShort_t q = 0; q < nRings; q++) {
308 Char_t r = (q == 0 ? 'I' : 'O');
309 UShort_t nsec = (q == 0 ? 20 : 40);
310 UShort_t nstr = (q == 0 ? 512 : 256);
311 RingHistos* histos = GetRingHistos(d, r);
313 for(UShort_t s = 0; s < nsec; s++) {
315 for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate;
317 #ifdef USE_OLDER_MERGING
318 Bool_t usedThis = kFALSE;
319 Bool_t usedPrev = kFALSE;
322 Bool_t used = kFALSE;
323 Double_t eTotal = -1;
324 Int_t nDistanceBefore = -1;
325 Int_t nDistanceAfter = -1;
326 Bool_t twoLow = kFALSE;
327 for(UShort_t t = 0; t < nstr; t++) {
331 output.SetMultiplicity(d,r,s,t,0.);
332 Float_t mult = SignalInStrip(input,d,r,s,t);
334 Float_t multNext = 0;
335 Float_t multNextNext = 0;
337 if(t<nstr-1) multNext = SignalInStrip(input,d,r,s,t+1);
338 if(t<nstr-2) multNextNext = SignalInStrip(input,d,r,s,t+2);
339 if(multNext == AliESDFMD::kInvalidMult) multNext = 0;
340 if(multNextNext == AliESDFMD::kInvalidMult) multNextNext = 0;
341 if(!fThreeStripSharing) multNextNext = 0;
342 // Get the pseudo-rapidity
343 Double_t eta = input.Eta(d,r,s,t);
344 Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
345 if (s == 0) output.SetEta(d,r,s,t,eta);
347 Double_t etaOld = eta;
348 Double_t etaCalc = 0;
349 if(fRecalculateEta) {
350 etaCalc = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
353 if(fRecalculateEta && mult > 0 && mult != AliESDFMD::kInvalidMult ) {
354 Double_t cosOld = TMath::Cos(2*TMath::ATan(TMath::Exp(-1*TMath::Abs(etaOld))));
355 Double_t cosNew = TMath::Cos(2*TMath::ATan(TMath::Exp(-1*TMath::Abs(etaCalc))));
356 if(mult > 0) mult = (mult/cosOld)*cosNew;
357 if(multNext > 0) multNext = (multNext/cosOld)*cosNew;
358 if(multNextNext > 0) multNextNext = (multNextNext/cosOld)*cosNew;
360 //No corrections beyond this point
362 // mult = 0; multNext = 0; multNextNext = 0;
366 // Keep dead-channel information.
367 if(mult == AliESDFMD::kInvalidMult)
368 output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
370 // If no signal or dead strip, go on.
371 if (mult == AliESDFMD::kInvalidMult || mult == 0) {
372 if (mult == 0) histos->fSum->Fill(eta,phi,mult);
377 // Fill the diagnostics histogram
378 histos->fBefore->Fill(mult);
380 Double_t mergedEnergy = 0;
382 if(fUseSimpleMerging) {
384 if (t < nstr-1) histos->fNeighborsBefore->Fill(mult,multNext);
385 if(mult > GetHighCut(d, r, eta ,false)) {
386 histos->fDistanceBefore->Fill(nDistanceBefore);
387 nDistanceBefore = -1;
391 if(fThreeStripSharing && multNext > GetLowCut(d, r, eta) &&
392 (multNext < GetHighCut(d, r, eta ,false) || twoLow)) {
393 eTotal = eTotal + multNext;
395 histos->fTriple->Fill(eTotal);
400 histos->fDouble->Fill(eTotal);
406 if(used) {used = kFALSE; continue; }
407 if(mult > GetLowCut(d, r, eta)) etot = mult;
409 if(mult > GetLowCut(d, r, eta) &&
410 multNext > GetLowCut(d, r, eta) &&
411 (mult < GetHighCut(d, r, eta ,false) ||
412 multNext < GetHighCut(d, r, eta ,false))) {
414 if(mult < GetHighCut(d, r, eta ,false) &&
415 multNext < GetHighCut(d, r, eta ,false) )
418 if(mult>multNext && multNextNext < GetLowCut(d, r, eta))
420 etot = mult + multNext;
422 histos->fDouble->Fill(etot);
426 eTotal = mult + multNext;
431 histos->fSingle->Fill(etot);
432 histos->fSinglePerStrip->Fill(etot,t);
438 if(mergedEnergy > GetHighCut(d, r, eta ,false) ) {
439 histos->fDistanceAfter->Fill(nDistanceAfter);
442 //if(mult>0 && multNext >0)
443 // std::cout<<mult<<" "<<multNext<<" "<<mergedEnergy<<std::endl;
446 // Get next and previous signal - if any
449 Status prevStatus = (t == 0 ? kNone : status[t-1]);
450 Status thisStatus = status[t];
451 Status nextStatus = (t == nstr-1 ? kNone : status[t+1]);
453 prevE = SignalInStrip(input,d,r,s,t-1);
454 if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
457 nextE = SignalInStrip(input,d,r,s,t+1);
458 if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
460 if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult);
462 #ifdef USE_OLDER_MERGING
463 /*Double_t*/ mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
466 status[t] = (usedPrev ? kMergedWithOther : kNone);
467 if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone);
469 /*Double_t*/ mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE,
475 if (t != 0) status[t-1] = prevStatus;
476 if (t != nstr-1) status[t+1] = nextStatus;
477 status[t] = thisStatus;
480 // If we're processing on non-angle corrected data, we
481 // should do the angle correction here
484 mergedEnergy = AngleCorrect(mergedEnergy, eta);
485 if (mergedEnergy > 0) histos->Incr();
488 histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1),
490 histos->fBeforeAfter->Fill(mult, mergedEnergy);
492 histos->fAfter->Fill(mergedEnergy);
493 histos->fSum->Fill(eta,phi,mergedEnergy);
495 output.SetMultiplicity(d,r,s,t,mergedEnergy);
497 for (UShort_t t = 0; t < nstr; t++) {
498 if (fOper) fOper->operator()(d, r, s, t) = status[t];
500 case kNone: nNone++; break;
501 case kCandidate: nCandidate++; break;
502 case kMergedWithOther: nMerged++; break;
503 case kMergedInto: nSummed++; break;
509 fSummed->Fill(kNone, nNone);
510 fSummed->Fill(kCandidate, nCandidate);
511 fSummed->Fill(kMergedWithOther, nMerged);
512 fSummed->Fill(kMergedInto, nSummed);
514 DBGL(5, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d",
515 nNone, nCandidate, nMerged, nSummed));
517 while ((o = static_cast<RingHistos*>(next())))
523 //_____________________________________________________________________
525 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
532 // Get the signal in a strip
544 Double_t mult = input.Multiplicity(d,r,s,t);
546 // - bad value (invalid or 0)
547 // - we want angle corrected and data is
548 // - we don't want angle corrected and data isn't
549 // just return read value
550 if (mult == AliESDFMD::kInvalidMult ||
552 (fCorrectAngles && input.IsAngleCorrected()) ||
553 (!fCorrectAngles && !input.IsAngleCorrected()))
556 // If we want angle corrected data, correct it,
557 // otherwise de-correct it
558 if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
559 else mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
562 //_____________________________________________________________________
564 AliFMDSharingFilter::GetLowCut(UShort_t d, Char_t r, Double_t eta) const
567 // Get the low cut. Normally, the low cut is taken to be the lower
568 // value of the fit range used when generating the energy loss fits.
569 // However, if fLowCut is set (using SetLowCit) to a value greater
570 // than 0, then that value is used.
572 return fLCuts.GetMultCut(d,r,eta,false);
574 if (!fCutAtFractionOfMPV && fLowCut > 0) return fLowCut;
576 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
577 AliFMDCorrELossFit* fits = fcm.GetELossFit();
579 if (fCutAtFractionOfMPV) {
580 AliFMDCorrELossFit::ELossFit* func = fits->GetFit(d,r,eta);
581 return fFractionOfMPV*func->GetDelta() ;
583 return fits->GetLowCut();
587 //_____________________________________________________________________
589 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r,
590 Double_t eta, Bool_t errors) const
593 // Get the high cut. The high cut is defined as the
594 // most-probably-value peak found from the energy distributions, minus
595 // 2 times the width of the corresponding Landau.
597 return fHCuts.GetMultCut(d,r,eta,errors);
600 //_____________________________________________________________________
602 const char* status2String(AliFMDSharingFilter::Status s)
605 case AliFMDSharingFilter::kNone: return "none ";
606 case AliFMDSharingFilter::kCandidate: return "candidate";
607 case AliFMDSharingFilter::kMergedWithOther: return "merged ";
608 case AliFMDSharingFilter::kMergedInto: return "summed ";
614 //_____________________________________________________________________
616 AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE,
627 Status& nextStatus) const
629 Double_t lowCut = GetLowCut(d, r, eta);
631 DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)",
633 thisE, status2String(thisStatus),
634 prevE, status2String(prevStatus),
635 nextE, status2String(nextStatus)));
637 // If below cut, then modify zero signal and make sure the next
638 // strip is considered a candidate.
639 if (thisE < lowCut || thisE > 20) {
641 DBGL(5,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE));
642 if (prevStatus == kCandidate) prevStatus = kNone;
645 // It this strip was merged with the previous strip, then
646 // make the next strip a candidate and zero the value in this strip.
647 if (thisStatus == kMergedWithOther) {
648 DBGL(5,Form(" Merged with other, 0'ed"));
652 // Get the upper cut on the sharing
653 Double_t highCut = GetHighCut(d, r, eta ,false);
659 // Only for low-flux events
661 // If this is smaller than the next, and the next is larger
662 // then the cut, mark both as candidates for sharing.
663 // They should be be merged on the next iteration
664 if (thisE < nextE && nextE > highCut) {
665 thisStatus = kCandidate;
666 nextStatus = kCandidate;
671 // Variable to sum signal in
672 Double_t totalE = thisE;
674 // If the previous signal was between the two cuts, and is still
675 // marked as a candidate , then merge it in here.
676 if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) {
678 prevStatus = kMergedWithOther;
679 DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f",
680 lowCut, prevE, highCut, totalE));
683 // If the next signal is between the two cuts, then merge it here
684 if (nextE > lowCut && nextE < highCut) {
686 nextStatus = kMergedWithOther;
687 DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE));
690 // If the signal is bigger than the high-cut, then return the value
691 if (totalE > highCut) {
693 thisStatus = kMergedInto;
696 DBGL(5, Form(" %9f>%f9 (was %9f) -> %9f %s",
697 totalE, highCut, thisE, totalE,status2String(thisStatus)));
701 // This is below cut, so flag next as a candidate
702 DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut));
703 nextStatus = kCandidate;
706 // If the total signal is smaller than low cut then zero this and kill this
707 if (totalE < lowCut) {
708 DBGL(5, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE));
713 // If total signal not above high cut or lower than low cut,
714 // mark this as a candidate for merging into the next, and zero signal
715 DBGL(5, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate",
716 lowCut, totalE, highCut, thisE));
717 thisStatus = kCandidate;
718 return (fZeroSharedHitsBelowThreshold ? 0 : totalE);
721 //_____________________________________________________________________
723 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
733 Bool_t& usedThis) const
736 // The actual algorithm
739 // mult The unfiltered signal in the strip
740 // eta Psuedo rapidity
741 // prevE Previous strip signal (or 0)
742 // nextE Next strip signal (or 0)
743 // lowFlux Whether this is a low flux event
748 // usedPrev Whether the previous strip was used in sharing or not
749 // usedThis Wether this strip was used in sharing or not.
752 // The filtered signal in the strip
755 // IF the multiplicity is very large, or below the cut, reject it, and
756 // flag it as candidate for sharing
757 Double_t lowCut = GetLowCut(d,r,eta);
758 if(mult > 12 || mult < lowCut) {
764 // If this was shared with the previous signal, return 0 and mark it as
765 // not a candiate for sharing
772 //analyse and perform sharing on one strip
775 Double_t highCut = GetHighCut(d, r, eta, false);
782 // If this signal is smaller than the next, and the next signal is smaller
783 // than then the high cut, and it's a low-flux event, then mark this and
784 // the next signal as candidates for sharing
786 // This is the test if the signal is the smaller of two possibly
796 Double_t totalE = mult;
798 // If the previous signal was larger than the low cut, and
799 // the previous was smaller than high cut, and the previous signal
800 // isn't marked as used, then add it's energy to this signal
801 if (prevE > lowCut &&
806 // If the next signal is larger than the low cut, and
807 // the next signal is smaller than the high cut, then add the next signal
808 // to this, and marked the next signal as used
809 if(nextE > lowCut && nextE < highCut ) {
814 // If we're processing on non-angle corrected data, we should do the
815 // angle correction here
816 // if (!fCorrectAngles)
817 // totalE = AngleCorrect(totalE, eta);
819 // Fill post processing histogram
820 // if(totalE > fLowCut)
821 // GetRingHistos(d,r)->fAfter->Fill(totalE);
824 Double_t mergedEnergy = 0;
827 // If we have a signal, then this is used (sharedPrev=true) and
828 // the signal is set to the result
829 mergedEnergy = totalE;
833 // If we do not have a signal (too low), then this is not shared,
834 // and the next strip is not shared either
841 //____________________________________________________________________
843 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
846 // Angle correct the signal
849 // mult Angle Un-corrected Signal
850 // eta Pseudo-rapidity
853 // Angle corrected signal
855 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
856 if (eta < 0) theta -= TMath::Pi();
857 return mult * TMath::Cos(theta);
859 //____________________________________________________________________
861 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
864 // Angle de-correct the signal
867 // mult Angle corrected Signal
868 // eta Pseudo-rapidity
871 // Angle un-corrected signal
873 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
874 if (eta < 0) theta -= TMath::Pi();
875 return mult / TMath::Cos(theta);
878 //____________________________________________________________________
880 AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
883 // Scale the histograms to the total number of events
886 // dir Where the output is
887 // nEvents Number of events
889 DGUARD(fDebug,1, "Scale histograms in AliFMDSharingFilter");
890 if (nEvents <= 0) return;
891 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
894 TIter next(&fRingHistos);
896 THStack* sums = new THStack("sums", "Sum of ring signals");
897 while ((o = static_cast<RingHistos*>(next()))) {
898 o->ScaleHistograms(d, nEvents);
899 TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
900 sum->Scale(1., "width");
901 sum->SetTitle(o->GetName());
902 sum->SetDirectory(0);
903 sum->SetYTitle("#sum #Delta/#Delta_{mip}");
909 //____________________________________________________________________
911 AliFMDSharingFilter::DefineOutput(TList* dir)
914 // Define the output histograms. These are put in a sub list of the
915 // passed list. The histograms are merged before the parent task calls
916 // AliAnalysisTaskSE::Terminate
919 // dir Directory to add to
921 DGUARD(fDebug,1, "Define output in AliFMDSharingFilter");
922 TList* d = new TList;
923 d->SetName(GetName());
926 fSummed = new TH2I("operations", "Strip operations",
927 kMergedInto, kNone-.5, kMergedInto+.5,
928 51201, -.5, 51200.5);
929 fSummed->SetXTitle("Operation");
930 fSummed->SetYTitle("# of strips");
931 fSummed->SetZTitle("Events");
932 fSummed->GetXaxis()->SetBinLabel(kNone, "None");
933 fSummed->GetXaxis()->SetBinLabel(kCandidate, "Candidate");
934 fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
935 fSummed->GetXaxis()->SetBinLabel(kMergedInto, "Merged into");
936 fSummed->SetDirectory(0);
939 fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
940 fHighCuts->SetXTitle("#eta");
941 fHighCuts->SetDirectory(0);
944 fLowCuts = new TH2D("lowCuts", "Low cuts used", 1,0,1, 1,0,1);
945 fLowCuts->SetXTitle("#eta");
946 fLowCuts->SetDirectory(0);
949 TNamed* angle = new TNamed("angle", fCorrectAngles ?
950 "corrected" : "uncorrected");
951 angle->SetUniqueID(fCorrectAngles);
952 TNamed* low = new TNamed("lowSignal", fZeroSharedHitsBelowThreshold ?
954 low->SetUniqueID(fZeroSharedHitsBelowThreshold);
955 TNamed* simple = new TNamed("simple", fUseSimpleMerging ? "yes" : "no");
956 simple->SetUniqueID(fUseSimpleMerging);
964 fLCuts.Output(d,"lCuts");
965 fHCuts.Output(d,"hCuts");
967 TIter next(&fRingHistos);
969 while ((o = static_cast<RingHistos*>(next()))) {
973 //____________________________________________________________________
975 AliFMDSharingFilter::Print(Option_t* /*option*/) const
983 char ind[gROOT->GetDirLevel()+1];
984 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
985 ind[gROOT->GetDirLevel()] = '\0';
986 std::cout << ind << ClassName() << ": " << GetName() << '\n'
988 << ind << " Debug: " << fDebug << "\n"
989 << ind << " Use corrected angles: " << fCorrectAngles << '\n'
990 << ind << " Zero below threshold: "
991 << fZeroSharedHitsBelowThreshold << '\n'
992 << ind << " Use simple sharing: " << fUseSimpleMerging << '\n'
993 << std::noboolalpha << std::endl;
994 std::cout << ind << " Low cuts: " << std::endl;
996 std::cout << ind << " High cuts: " << std::endl;
1000 //====================================================================
1001 AliFMDSharingFilter::RingHistos::RingHistos()
1002 : AliForwardUtil::RingHistos(),
1012 fNeighborsBefore(0),
1024 //____________________________________________________________________
1025 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
1026 : AliForwardUtil::RingHistos(d,r),
1036 fNeighborsBefore(0),
1049 fBefore = new TH1D("esdEloss", Form("Energy loss in %s (reconstruction)",
1050 GetName()), 600, 0, 15);
1051 fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
1052 fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
1053 fBefore->SetFillColor(Color());
1054 fBefore->SetFillStyle(3001);
1055 fBefore->SetLineColor(kBlack);
1056 fBefore->SetLineStyle(2);
1057 fBefore->SetDirectory(0);
1059 fAfter = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
1060 fAfter->SetTitle(Form("Energy loss in %s (sharing corrected)", GetName()));
1061 fAfter->SetFillColor(Color()+2);
1062 fAfter->SetLineStyle(1);
1063 fAfter->SetDirectory(0);
1065 fSingle = new TH1D("singleEloss", "Energy loss (single strips)",
1067 fSingle->SetXTitle("#Delta/#Delta_{mip}");
1068 fSingle->SetYTitle("P(#Delta/#Delta_{mip})");
1069 fSingle->SetFillColor(Color());
1070 fSingle->SetFillStyle(3001);
1071 fSingle->SetLineColor(kBlack);
1072 fSingle->SetLineStyle(2);
1073 fSingle->SetDirectory(0);
1075 fDouble = static_cast<TH1D*>(fSingle->Clone("doubleEloss"));
1076 fDouble->SetTitle("Energy loss (two strips)");
1077 fDouble->SetFillColor(Color()+1);
1078 fDouble->SetDirectory(0);
1080 fTriple = static_cast<TH1D*>(fSingle->Clone("tripleEloss"));
1081 fTriple->SetTitle("Energy loss (three strips)");
1082 fTriple->SetFillColor(Color()+2);
1083 fTriple->SetDirectory(0);
1085 //Int_t nBinsForInner = (r == 'I' ? 32 : 16);
1086 Int_t nBinsForInner = (r == 'I' ? 512 : 256);
1087 Int_t nStrips = (r == 'I' ? 512 : 256);
1089 fSinglePerStrip = new TH2D("singlePerStrip", "SinglePerStrip",
1090 600,0,15, nBinsForInner,0,nStrips);
1091 fSinglePerStrip->SetXTitle("#Delta/#Delta_{mip}");
1092 fSinglePerStrip->SetYTitle("Strip #");
1093 fSinglePerStrip->SetZTitle("Counts");
1094 fSinglePerStrip->SetDirectory(0);
1096 fDistanceBefore = new TH1D("distanceBefore", "Distance before sharing",
1097 nStrips , 0,nStrips );
1098 fDistanceBefore->SetXTitle("Distance");
1099 fDistanceBefore->SetYTitle("Counts");
1100 fDistanceBefore->SetFillColor(kGreen+2);
1101 fDistanceBefore->SetFillStyle(3001);
1102 fDistanceBefore->SetLineColor(kBlack);
1103 fDistanceBefore->SetLineStyle(2);
1104 fDistanceBefore->SetDirectory(0);
1106 fDistanceAfter = static_cast<TH1D*>(fDistanceBefore->Clone("distanceAfter"));
1107 fDistanceAfter->SetTitle("Distance after sharing");
1108 fDistanceAfter->SetFillColor(kGreen+1);
1109 fDistanceAfter->SetDirectory(0);
1114 Int_t n = int((max-min) / (max / 300));
1115 fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation",
1116 n, min, max, n, min, max);
1117 fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
1118 fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
1119 fBeforeAfter->SetZTitle("Correlation");
1120 fBeforeAfter->SetDirectory(0);
1122 fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
1123 fNeighborsBefore->SetTitle("Correlation of neighbors before");
1124 fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
1125 fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
1126 fNeighborsBefore->SetDirectory(0);
1129 static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
1130 fNeighborsAfter->SetTitle("Correlation of neighbors after");
1131 fNeighborsAfter->SetDirectory(0);
1133 fSum = new TH2D("summed", "Summed signal", 200, -4, 6,
1134 (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1135 fSum->SetDirectory(0);
1137 fSum->SetMarkerColor(Color());
1138 // fSum->SetFillColor(Color());
1139 fSum->SetXTitle("#eta");
1140 fSum->SetYTitle("#varphi [radians]");
1141 fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) ");
1143 fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
1144 fHits->SetDirectory(0);
1145 fHits->SetXTitle("# of hits");
1146 fHits->SetFillColor(kGreen+1);
1148 //____________________________________________________________________
1149 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
1150 : AliForwardUtil::RingHistos(o),
1156 fSinglePerStrip(o.fSinglePerStrip),
1157 fDistanceBefore(o.fDistanceBefore),
1158 fDistanceAfter(o.fDistanceAfter),
1159 fBeforeAfter(o.fBeforeAfter),
1160 fNeighborsBefore(o.fNeighborsBefore),
1161 fNeighborsAfter(o.fNeighborsAfter),
1170 // o Object to copy from
1174 //____________________________________________________________________
1175 AliFMDSharingFilter::RingHistos&
1176 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
1179 // Assignment operator
1182 // o Object to assign from
1185 // Reference to this
1187 if (&o == this) return *this;
1188 AliForwardUtil::RingHistos::operator=(o);
1192 if (fBefore) delete fBefore;
1193 if (fAfter) delete fAfter;
1194 if (fSingle) delete fSingle;
1195 if (fDouble) delete fDouble;
1196 if (fTriple) delete fTriple;
1197 if (fSinglePerStrip) delete fSinglePerStrip;
1198 if (fDistanceBefore) delete fDistanceBefore;
1199 if (fDistanceAfter) delete fDistanceAfter;
1200 if (fHits) delete fHits;
1203 fBefore = static_cast<TH1D*>(o.fBefore->Clone());
1204 fAfter = static_cast<TH1D*>(o.fAfter->Clone());
1205 fSingle = static_cast<TH1D*>(o.fSingle->Clone());
1206 fDouble = static_cast<TH1D*>(o.fDouble->Clone());
1207 fTriple = static_cast<TH1D*>(o.fTriple->Clone());
1208 fSinglePerStrip = static_cast<TH2D*>(o.fSinglePerStrip->Clone());
1209 fDistanceBefore = static_cast<TH1D*>(o.fDistanceBefore->Clone());
1210 fDistanceAfter = static_cast<TH1D*>(o.fDistanceAfter->Clone());
1211 fBeforeAfter = static_cast<TH2D*>(o.fBeforeAfter->Clone());
1212 fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
1213 fNeighborsAfter = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
1214 fHits = static_cast<TH1D*>(o.fHits->Clone());
1215 fSum = static_cast<TH2D*>(o.fSum->Clone());
1219 //____________________________________________________________________
1220 AliFMDSharingFilter::RingHistos::~RingHistos()
1226 //____________________________________________________________________
1228 AliFMDSharingFilter::RingHistos::Finish()
1234 fHits->Fill(fNHits);
1237 //____________________________________________________________________
1239 AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir,
1243 // Scale the histograms to the total number of events
1246 // nEvents Number of events
1247 // dir Where the output is
1249 TList* l = GetOutputList(dir);
1253 TH1D* before = static_cast<TH1D*>(l->FindObject("esdEloss"));
1254 TH1D* after = static_cast<TH1D*>(l->FindObject("anaEloss"));
1255 if (before) before->Scale(1./nEvents);
1256 if (after) after->Scale(1./nEvents);
1259 TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
1260 if (summed) summed->Scale(1./nEvents);
1264 //____________________________________________________________________
1266 AliFMDSharingFilter::RingHistos::Output(TList* dir)
1272 // dir where to store
1274 TList* d = DefineOutputList(dir);
1281 d->Add(fSinglePerStrip);
1282 d->Add(fDistanceBefore);
1283 d->Add(fDistanceAfter);
1284 d->Add(fBeforeAfter);
1285 d->Add(fNeighborsBefore);
1286 d->Add(fNeighborsAfter);
1290 // Removed to avoid doubly adding the list which destroys
1295 //____________________________________________________________________