2 // Class to do the sharing correction of FMD ESD data
4 #include "AliFMDSharingFilter.h"
10 #include "AliForwardCorrectionManager.h"
11 // #include "AliFMDAnaParameters.h"
17 ClassImp(AliFMDSharingFilter)
19 ; // This is for Emacs
23 //____________________________________________________________________
24 AliFMDSharingFilter::AliFMDSharingFilter()
28 fCorrectAngles(kFALSE),
33 //____________________________________________________________________
34 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
35 : TNamed("fmdSharingFilter", title),
38 fCorrectAngles(kFALSE),
42 fRingHistos.Add(new RingHistos(1, 'I'));
43 fRingHistos.Add(new RingHistos(2, 'I'));
44 fRingHistos.Add(new RingHistos(2, 'O'));
45 fRingHistos.Add(new RingHistos(3, 'I'));
46 fRingHistos.Add(new RingHistos(3, 'O'));
49 //____________________________________________________________________
50 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
54 fCorrectAngles(o.fCorrectAngles),
58 TIter next(&o.fRingHistos);
60 while ((obj = next())) fRingHistos.Add(obj);
63 //____________________________________________________________________
64 AliFMDSharingFilter::~AliFMDSharingFilter()
69 //____________________________________________________________________
71 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
76 fCorrectAngles = o.fCorrectAngles;
81 TIter next(&o.fRingHistos);
83 while ((obj = next())) fRingHistos.Add(obj);
88 //____________________________________________________________________
89 AliFMDSharingFilter::RingHistos*
90 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
94 case 1: idx = 0; break;
95 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
96 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
98 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
100 return static_cast<RingHistos*>(fRingHistos.At(idx));
103 //____________________________________________________________________
105 AliFMDSharingFilter::Filter(const AliESDFMD& input,
110 TIter next(&fRingHistos);
112 Double_t lowCut = GetLowCut();
113 while ((o = static_cast<RingHistos*>(next())))
116 // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
118 for(UShort_t d = 1; d <= 3; d++) {
119 Int_t nRings = (d == 1 ? 1 : 2);
120 for (UShort_t q = 0; q < nRings; q++) {
121 Char_t r = (q == 0 ? 'I' : 'O');
122 UShort_t nsec = (q == 0 ? 20 : 40);
123 UShort_t nstr = (q == 0 ? 512 : 256);
125 RingHistos* histos = GetRingHistos(d, r);
127 for(UShort_t s = 0; s < nsec; s++) {
128 Bool_t usedThis = kFALSE;
129 Bool_t usedPrev = kFALSE;
131 for(UShort_t t = 0; t < nstr; t++) {
132 output.SetMultiplicity(d,r,s,t,0.);
133 Float_t mult = SignalInStrip(input,d,r,s,t);
135 // Keep dead-channel information.
136 if(mult == AliESDFMD::kInvalidMult)
137 output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
139 // If no signal or dead strip, go on.
140 if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
142 // Get the pseudo-rapidity
143 Double_t eta = input.Eta(d,r,s,t);
145 // Fill the diagnostics histogram
146 histos->fBefore->Fill(mult);
148 // Get next and previous signal - if any
152 prevE = SignalInStrip(input,d,r,s,t-1);
153 if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
156 nextE = SignalInStrip(input,d,r,s,t+1);
157 if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
160 Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
163 if (mergedEnergy > lowCut) histos->Incr();
164 histos->fAfter->Fill(mergedEnergy);
166 output.SetMultiplicity(d,r,s,t,mergedEnergy);
167 output.SetEta(d,r,s,t,eta);
174 while ((o = static_cast<RingHistos*>(next())))
180 //_____________________________________________________________________
182 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
188 Double_t mult = input.Multiplicity(d,r,s,t);
189 if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
191 if (fCorrectAngles && !input.IsAngleCorrected())
192 mult = AngleCorrect(mult, input.Eta(d,r,s,t));
193 else if (!fCorrectAngles && input.IsAngleCorrected())
194 mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
197 //_____________________________________________________________________
199 AliFMDSharingFilter::GetLowCut() const
201 if (fLowCut > 0) return fLowCut;
202 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
203 AliFMDCorrELossFit* fits = fcm.GetELossFit();
204 return fits->GetLowCut();
207 //_____________________________________________________________________
209 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
211 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
214 // Get the high cut. The high cut is defined as the
215 // most-probably-value peak found from the energy distributions, minus
216 // 2 times the width of the corresponding Landau.
217 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
218 Double_t delta = 1024;
221 AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
229 AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
231 Double_t highCut = (delta - fNXi * xi);
235 //_____________________________________________________________________
237 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
249 // IF the multiplicity is very large, or below the cut, reject it, and
250 // flag it as candidate for sharing
251 Double_t lowCut = GetLowCut();
252 if(mult > 12 || mult < lowCut) {
258 // If this was shared with the previous signal, return 0 and mark it as
259 // not a candiate for sharing
266 //analyse and perform sharing on one strip
269 Double_t highCut = GetHighCut(d, r, eta);
271 // If this signal is smaller than the next, and the next signal is smaller
272 // than then the high cut, and it's a low-flux event, then mark this and
273 // the next signal as candidates for sharing
275 // This is the test if the signal is the smaller of two possibly
285 Double_t totalE = mult;
287 // If the previous signal was larger than the low cut, and
288 // the previous was smaller than high cut, and the previous signal
289 // isn't marked as used, then add it's energy to this signal
290 if (prevE > lowCut &&
295 // If the next signal is larger than the low cut, and
296 // the next signal is smaller than the low cut, then add the next signal
297 // to this, and marked the next signal as used
298 if(nextE > lowCut && nextE < highCut ) {
303 // If we're processing on non-angle corrected data, we should do the
304 // angle correction here
306 totalE = AngleCorrect(totalE, eta);
308 // Fill post processing histogram
309 // if(totalE > fLowCut)
310 // GetRingHistos(d,r)->fAfter->Fill(totalE);
313 Double_t mergedEnergy = 0;
316 // If we have a signal, then this is used (sharedPrev=true) and
317 // the signal is set to the result
318 mergedEnergy = totalE;
322 // If we do not have a signal (too low), then this is not shared,
323 // and the next strip is not shared either
330 //____________________________________________________________________
332 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
334 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
335 if (eta < 0) theta -= TMath::Pi();
336 return mult * TMath::Cos(theta);
338 //____________________________________________________________________
340 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
342 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
343 if (eta < 0) theta -= TMath::Pi();
344 return mult / TMath::Cos(theta);
347 //____________________________________________________________________
349 AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
351 if (nEvents <= 0) return;
352 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
355 TIter next(&fRingHistos);
357 while ((o = static_cast<RingHistos*>(next())))
358 o->ScaleHistograms(d, nEvents);
361 //____________________________________________________________________
363 AliFMDSharingFilter::DefineOutput(TList* dir)
365 TList* d = new TList;
366 d->SetName(GetName());
368 TIter next(&fRingHistos);
370 while ((o = static_cast<RingHistos*>(next()))) {
374 //____________________________________________________________________
376 AliFMDSharingFilter::Print(Option_t* /*option*/) const
378 char ind[gROOT->GetDirLevel()+1];
379 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
380 ind[gROOT->GetDirLevel()] = '\0';
381 std::cout << ind << "AliFMDSharingFilter: " << GetName() << '\n'
382 << ind << " Low cut: " << fLowCut << '\n'
383 << ind << " N xi factor: " << fNXi << '\n'
384 << ind << " Use corrected angles: "
385 << (fCorrectAngles ? "yes" : "no") << std::endl;
388 //====================================================================
389 AliFMDSharingFilter::RingHistos::RingHistos()
390 : AliForwardUtil::RingHistos(),
397 //____________________________________________________________________
398 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
399 : AliForwardUtil::RingHistos(d,r),
405 fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
406 Form("Energy loss in %s (reconstruction)", fName.Data()),
408 fAfter = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
409 Form("Energy loss in %s (sharing corrected)",
410 fName.Data()), 1000, 0, 25);
411 fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
412 fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
413 fBefore->SetFillColor(kRed+1);
414 fBefore->SetFillStyle(3001);
416 fBefore->SetDirectory(0);
417 fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
418 fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
419 fAfter->SetFillColor(kBlue+1);
420 fAfter->SetFillStyle(3001);
422 fAfter->SetDirectory(0);
424 fHits = new TH1D(Form("%s_Hits", fName.Data()),
425 Form("Number of hits in %s", fName.Data()),
427 fHits->SetDirectory(0);
428 fHits->SetXTitle("# of hits");
429 fHits->SetFillColor(kGreen+1);
431 //____________________________________________________________________
432 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
433 : AliForwardUtil::RingHistos(o),
440 //____________________________________________________________________
441 AliFMDSharingFilter::RingHistos&
442 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
444 AliForwardUtil::RingHistos::operator=(o);
448 if (fBefore) delete fBefore;
449 if (fAfter) delete fAfter;
450 if (fHits) delete fHits;
452 fBefore = static_cast<TH1D*>(o.fBefore->Clone());
453 fAfter = static_cast<TH1D*>(o.fAfter->Clone());
454 fHits = static_cast<TH1D*>(o.fHits->Clone());
458 //____________________________________________________________________
459 AliFMDSharingFilter::RingHistos::~RingHistos()
461 if (fBefore) delete fBefore;
462 if (fAfter) delete fAfter;
463 if (fHits) delete fHits;
465 //____________________________________________________________________
467 AliFMDSharingFilter::RingHistos::Finish()
472 //____________________________________________________________________
474 AliFMDSharingFilter::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
476 TList* l = GetOutputList(dir);
479 TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
481 TH1D* after = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
483 if (before) before->Scale(1./nEvents);
484 if (after) after->Scale(1./nEvents);
487 //____________________________________________________________________
489 AliFMDSharingFilter::RingHistos::Output(TList* dir)
491 TList* d = DefineOutputList(dir);
498 //____________________________________________________________________