AliPhysics  63e47e1 (63e47e1)
AliFMDSharingFilter.cxx
Go to the documentation of this file.
1 //
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 //
23 //
24 #include "AliFMDSharingFilter.h"
25 #include "AliFMDStripIndex.h"
26 #include <AliESDFMD.h>
27 #include <TAxis.h>
28 #include <TList.h>
29 #include <TH1.h>
30 #include <TMath.h>
32 #include "AliFMDCorrELossFit.h"
33 #include <AliLog.h>
34 #include <TROOT.h>
35 #include <THStack.h>
36 #include <TParameter.h>
37 #include <iostream>
38 #include <iomanip>
39 
40 ClassImp(AliFMDSharingFilter)
41 #if 0
42 ; // This is for Emacs
43 #endif
44 
45 #define DBG(L,M) \
46  do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false)
47 #define DBGL(L,M) \
48  do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false)
49 
50 
51 
52 //____________________________________________________________________
54  : TNamed(),
55  fRingHistos(),
56  fCorrectAngles(kFALSE),
57  // fSummed(0),
58  fHighCuts(0),
59  fLowCuts(0),
60  // fOper(0),
61  fDebug(0),
62  fZeroSharedHitsBelowThreshold(false),
63  fLCuts(),
64  fHCuts(),
65  fUseSimpleMerging(false),
66  fThreeStripSharing(true),
67  fMergingDisabled(false),
68  fIgnoreESDForAngleCorrection(false)
69 {
70  //
71  // Default Constructor - do not use
72  //
73  DGUARD(fDebug,1, "Default CTOR for AliFMDSharingFilter");
74 }
75 
76 //____________________________________________________________________
78  : TNamed("fmdSharingFilter", title),
79  fRingHistos(),
80  fCorrectAngles(kFALSE),
81  // fSummed(0),
82  fHighCuts(0),
83  fLowCuts(0),
84  // fOper(0),
85  fDebug(0),
87  fLCuts(),
88  fHCuts(),
89  fUseSimpleMerging(false),
90  fThreeStripSharing(true),
91  fMergingDisabled(false),
93 {
94  //
95  // Constructor
96  //
97  // Parameters:
98  // title Title of object - not significant
99  //
100  DGUARD(fDebug,1, "Named CTOR for AliFMDSharingFilter: %s", title);
101  fRingHistos.SetName(GetName());
102  fRingHistos.SetOwner();
103 
106 
107  // fExtraDead.Reset(-1);
108 }
109 
110 //____________________________________________________________________
112 {
113  //
114  // Destructor
115  //
116  DGUARD(fDebug,3, "DTOR for AliFMDSharingFilter");
117  // fRingHistos.Delete();
118 }
119 
120 //____________________________________________________________________
123 {
124  //
125  // Get the ring histogram container
126  //
127  // Parameters:
128  // d Detector
129  // r Ring
130  //
131  // Return:
132  // Ring histogram container
133  //
134  Int_t idx = -1;
135  switch (d) {
136  case 1: idx = 0; break;
137  case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
138  case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
139  }
140  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
141 
142  return static_cast<RingHistos*>(fRingHistos.At(idx));
143 }
144 
145 //____________________________________________________________________
146 void
148 {
149  // Initialise - called on first event
150  DGUARD(fDebug,1, "Initialize for AliFMDSharingFilter");
152  const AliFMDCorrELossFit* fits = fcm.GetELossFit();
153 
154  // Get the high cut. The high cut is defined as the
155  // most-probably-value peak found from the energy distributions, minus
156  // 2 times the width of the corresponding Landau.
157 
158  TAxis eAxis(axis.GetNbins(),
159  axis.GetXmin(),
160  axis.GetXmax());
161  if(fits)
162  eAxis.Set(fits->GetEtaAxis().GetNbins(),
163  fits->GetEtaAxis().GetXmin(),
164  fits->GetEtaAxis().GetXmax());
165 
166  UShort_t nEta = eAxis.GetNbins();
167 
168  fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
169  fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
170  fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
171  fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
172  fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
173  fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
174 
175  fLowCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
176  fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
177  fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
178  fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
179  fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
180  fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
181 
182  // Cache our cuts in histograms
185 }
186 
187 //____________________________________________________________________
188 #define ETA2COS(ETA) \
189  TMath::Cos(2*TMath::ATan(TMath::Exp(-TMath::Abs(ETA))))
190 
191 Bool_t
193  Bool_t /*lowFlux*/,
194  AliESDFMD& output,
195  Double_t /*zvtx*/)
196 {
197  //
198  // Filter the input AliESDFMD object
199  //
200  // Parameters:
201  // input Input
202  // lowFlux If this is a low-flux event
203  // output Output AliESDFMD object
204  //
205  // Return:
206  // True on success, false otherwise
207  //
208  DGUARD(fDebug,1, "Filter event in AliFMDSharingFilter");
209  output.Clear();
210  TIter next(&fRingHistos);
211  RingHistos* o = 0;
212  while ((o = static_cast<RingHistos*>(next()))) o->Clear();
213 
214  Int_t nSingle = 0;
215  Int_t nDouble = 0;
216  Int_t nTriple = 0;
217 
218  for(UShort_t d = 1; d <= 3; d++) {
219  Int_t nRings = (d == 1 ? 1 : 2);
220  for (UShort_t q = 0; q < nRings; q++) {
221  Char_t r = (q == 0 ? 'I' : 'O');
222  UShort_t nsec = (q == 0 ? 20 : 40);
223  UShort_t nstr = (q == 0 ? 512 : 256);
224  RingHistos* histos = GetRingHistos(d, r);
225 
226  for(UShort_t s = 0; s < nsec; s++) {
227  // `used' flags if the _current_ strip was used by _previous_
228  // iteration.
229  Bool_t used = kFALSE;
230  // `eTotal' contains the current sum of merged signals so far
231  Double_t eTotal = -1;
232  // Int_t nDistanceBefore = -1;
233  // Int_t nDistanceAfter = -1;
234  // `twoLow' flags if we saw two consequtive strips with a
235  // signal between the two cuts.
236  Bool_t twoLow = kFALSE;
237  Int_t nStripsAboveCut = 0;
238 
239  for(UShort_t t = 0; t < nstr; t++) {
240  // nDistanceBefore++;
241  // nDistanceAfter++;
242 
243  output.SetMultiplicity(d,r,s,t,0.);
244  Float_t mult = SignalInStrip(input,d,r,s,t);
245  Float_t multNext = (t<nstr-1) ? SignalInStrip(input,d,r,s,t+1) :0;
246  Float_t multNextNext = (t<nstr-2) ? SignalInStrip(input,d,r,s,t+2) :0;
247  if (multNext == AliESDFMD::kInvalidMult) multNext = 0;
248  if (multNextNext == AliESDFMD::kInvalidMult) multNextNext = 0;
249  if(!fThreeStripSharing) multNextNext = 0;
250 
251  // Get the pseudo-rapidity
252  Double_t eta = input.Eta(d,r,s,t);
253  Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
254  if (s == 0) output.SetEta(d,r,s,t,eta);
255 
256  // Keep dead-channel information - either from the ESD (but
257  // see above for older data) or from the settings in the
258  // ForwardAODConfig.C file.
259  if (mult == AliESDFMD::kInvalidMult) {
260  output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
261  histos->fBefore->Fill(-1);
262  mult = AliESDFMD::kInvalidMult;
263  }
264 
265  Double_t lowCut = GetLowCut(d, r, eta);
266  Double_t highCut = GetHighCut(d, r, eta, false);
267  if (mult != AliESDFMD::kInvalidMult && mult > lowCut) {
268  // Always fill the ESD sum histogram
269  histos->fSumESD->Fill(eta, phi, mult);
270  }
271 
272  // If no signal or dead strip, go on.
273  if (mult == AliESDFMD::kInvalidMult || mult == 0) {
274  if (mult == 0) histos->fSum->Fill(eta,phi,mult);
275  // Flush a possible signal
276  if (eTotal > 0 && t > 0)
277  output.SetMultiplicity(d,r,s,t-1,eTotal);
278  // Reset states so we do not try to merge over a dead strip.
279  eTotal = -1;
280  used = false;
281  twoLow = false;
282  if (t > 0)
283  histos->fNConsecutive->Fill(nStripsAboveCut);
284  if (mult == AliESDFMD::kInvalidMult)
285  // Why not fill immidiately here?
286  nStripsAboveCut = -1;
287  else
288  // Why not fill immidiately here?
289  nStripsAboveCut = 0;
290  continue;
291  }
292 
293  // Fill the diagnostics histogram
294  histos->fBefore->Fill(mult);
295 
296  Double_t mergedEnergy = mult;
297  // it seems to me that this logic could be condensed a bit
298  if(mult > lowCut) {
299  if(nStripsAboveCut < 1) {
300  if(t > 0)
301  histos->fNConsecutive->Fill(nStripsAboveCut);
302  nStripsAboveCut=0;
303  }
304  nStripsAboveCut++;
305  }
306  else {
307  if (t > 0)
308  histos->fNConsecutive->Fill(nStripsAboveCut);
309  nStripsAboveCut=0;
310  }
311 
312  if (!fMergingDisabled) {
313  mergedEnergy = 0;
314 
315  // The current sum
316  Float_t etot = 0;
317 
318  // Fill in neighbor information
319  if (t < nstr-1) histos->fNeighborsBefore->Fill(mult,multNext);
320 
321  Bool_t thisValid = mult > lowCut;
322  Bool_t nextValid = multNext > lowCut;
323  Bool_t thisSmall = mult < highCut;
324  Bool_t nextSmall = multNext < highCut;
325 
326  // If this strips signal is above the high cut, reset distance
327  // if (!thisSmall) {
328  // histos->fDistanceBefore->Fill(nDistanceBefore);
329  // nDistanceBefore = -1;
330  // }
331 
332  // If the total signal in the past 1 or 2 strips are non-zero
333  // we need to check
334  if (eTotal > 0) {
335  // Here, we have already flagged one strip as a candidate
336 
337  // If 3-strip merging is enabled, then check the next
338  // strip to see that it falls within cut, or if we have
339  // two low signals
340  if (fThreeStripSharing && nextValid && (nextSmall || twoLow)) {
341  eTotal = eTotal + multNext;
342  used = kTRUE;
343  histos->fTriple->Fill(eTotal);
344  nTriple++;
345  twoLow = kFALSE;
346  }
347  // Otherwise, we got a double hit before, and that
348  // should be stored.
349  else {
350  used = kFALSE;
351  histos->fDouble->Fill(eTotal);
352  nDouble++;
353  }
354  // Store energy loss and reset sum
355  etot = eTotal;
356  eTotal = -1;
357  } // if (eTotal>0)
358  else {
359  // If we have no current sum
360 
361  // Check if this is marked as used, and if so, continue
362  if (used) {used = kFALSE; continue; }
363 
364  // If the signal is abvoe the cut, set current
365  if (thisValid) etot = mult;
366 
367  // If the signal is abiove the cut, and so is the next
368  // signal and either of them are below the high cut,
369  if (thisValid && nextValid && (thisSmall || nextSmall)) {
370 
371  // If this is below the high cut, and the next is too, then
372  // we have two low signals
373  if (thisSmall && nextSmall) twoLow = kTRUE;
374 
375  // If this signal is bigger than the next, and the
376  // one after that is below the low-cut, then update
377  // the sum
378  if (mult>multNext && multNextNext < lowCut) {
379  etot = mult + multNext;
380  used = kTRUE;
381  histos->fDouble->Fill(etot);
382  nDouble++;
383  }
384  // Otherwise, we may need to merge with a third strip
385  else {
386  etot = 0;
387  eTotal = mult + multNext;
388  }
389  }
390  // This is a signle hit
391  else if(etot > 0) {
392  histos->fSingle->Fill(etot);
393  histos->fSinglePerStrip->Fill(etot,t);
394  nSingle++;
395  }
396  } // else if (etotal >= 0)
397 
398  mergedEnergy = etot;
399  // if (mergedEnergy > GetHighCut(d, r, eta ,false)) {
400  // histos->fDistanceAfter->Fill(nDistanceAfter);
401  // nDistanceAfter = -1;
402  // }
403  //if(mult>0 && multNext >0)
404  // std::cout<<mult<<" "<<multNext<<" "<<mergedEnergy<<std::endl;
405  } // if (!fMergingDisabled)
406 
407  if (!fCorrectAngles)
408  mergedEnergy = AngleCorrect(mergedEnergy, eta);
409  // if (mergedEnergy > 0) histos->Incr();
410 
411  if (t != 0)
412  histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1),
413  mergedEnergy);
414  histos->fBeforeAfter->Fill(mult, mergedEnergy);
415  if(mergedEnergy > 0)
416  histos->fAfter->Fill(mergedEnergy);
417  histos->fSum->Fill(eta,phi,mergedEnergy);
418 
419  output.SetMultiplicity(d,r,s,t,mergedEnergy);
420  } // for strip
421  histos->fNConsecutive->Fill(nStripsAboveCut); // fill the last sector
422  } // for sector
423  } // for ring
424  } // for detector
425  DMSG(fDebug, 3,"single=%9d, double=%9d, triple=%9d",
426  nSingle, nDouble, nTriple);
427  next.Reset();
428  // while ((o = static_cast<RingHistos*>(next()))) o->Finish();
429 
430  return kTRUE;
431 }
432 
433 //_____________________________________________________________________
434 Double_t
436  UShort_t d,
437  Char_t r,
438  UShort_t s,
439  UShort_t t) const
440 {
441  //
442  // Get the signal in a strip
443  //
444  // Parameters:
445  // fmd ESD object
446  // d Detector
447  // r Ring
448  // s Sector
449  // t Strip
450  //
451  // Return:
452  // The energy signal
453  //
454  Double_t mult = input.Multiplicity(d,r,s,t);
455  // In case of
456  // - bad value (invalid or 0)
457  // - we want angle corrected and data is
458  // - we don't want angle corrected and data isn't
459  // just return read value
460  if (mult == AliESDFMD::kInvalidMult ||
461  mult == 0 ||
462  (fCorrectAngles && (fIgnoreESDForAngleCorrection || input.IsAngleCorrected())) ||
463  (!fCorrectAngles && !fIgnoreESDForAngleCorrection && !input.IsAngleCorrected()))
464  return mult;
465 
466  // If we want angle corrected data, correct it,
467  // otherwise de-correct it
468  if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
469  else mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
470  return mult;
471 }
472 
473 namespace {
474  Double_t Rng2Cut(UShort_t d, Char_t r, Double_t eta, TH2* h) {
475  Double_t ret = 1024;
476  Int_t ybin = 0;
477  switch(d) {
478  case 1: ybin = 1; break;
479  case 2: ybin = (r=='i' || r=='I') ? 2 : 3; break;
480  case 3: ybin = (r=='i' || r=='I') ? 4 : 5; break;
481  default: return ret;
482  }
483  Int_t xbin = h->GetXaxis()->FindBin(eta);
484  if (xbin < 1 && xbin > h->GetXaxis()->GetNbins()) return ret;
485  ret = h->GetBinContent(xbin,ybin);
486  return ret;
487  }
488 }
489 
490  //_____________________________________________________________________
491 Double_t
493 {
494  //
495  // Get the low cut. Normally, the low cut is taken to be the lower
496  // value of the fit range used when generating the energy loss fits.
497  // However, if fLowCut is set (using SetLowCit) to a value greater
498  // than 0, then that value is used.
499  //
500  return Rng2Cut(d, r, eta, fLowCuts);
501  // return fLCuts.GetMultCut(d,r,eta,false);
502 }
503 
504 //_____________________________________________________________________
505 Double_t
507  Double_t eta, Bool_t /*errors*/) const
508 {
509  //
510  // Get the high cut. The high cut is defined as the
511  // most-probably-value peak found from the energy distributions, minus
512  // 2 times the width of the corresponding Landau.
513  //
514  return Rng2Cut(d, r, eta, fHighCuts);
515  // return fHCuts.GetMultCut(d,r,eta,errors);
516 }
517 
518 //____________________________________________________________________
519 Double_t
521 {
522  //
523  // Angle correct the signal
524  //
525  // Parameters:
526  // mult Angle Un-corrected Signal
527  // eta Pseudo-rapidity
528  //
529  // Return:
530  // Angle corrected signal
531  //
532  Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
533  if (eta < 0) theta -= TMath::Pi();
534  return mult * TMath::Cos(theta);
535 }
536 //____________________________________________________________________
537 Double_t
539 {
540  //
541  // Angle de-correct the signal
542  //
543  // Parameters:
544  // mult Angle corrected Signal
545  // eta Pseudo-rapidity
546  //
547  // Return:
548  // Angle un-corrected signal
549  //
550  Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
551  if (eta < 0) theta -= TMath::Pi();
552  return mult / TMath::Cos(theta);
553 }
554 
555 //____________________________________________________________________
556 void
558 {
559  //
560  // Scale the histograms to the total number of events
561  //
562  // Parameters:
563  // dir Where the output is
564  // nEvents Number of events
565  //
566  DGUARD(fDebug,1, "Scale histograms in AliFMDSharingFilter");
567  if (nEvents <= 0) return;
568  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
569  if (!d) return;
570 
571  TList* out = new TList;
572  out->SetName(d->GetName());
573  out->SetOwner();
574 
575  TParameter<int>* nFiles =
576  static_cast<TParameter<int>*>(d->FindObject("nFiles"));
577 
578  TH2* lowCuts = static_cast<TH2*>(d->FindObject("lowCuts"));
579  TH2* highCuts = static_cast<TH2*>(d->FindObject("highCuts"));
580  if (lowCuts && nFiles) {
581  TH1* oh = static_cast<TH1*>(lowCuts->Clone());
582  oh->Scale(1. / nFiles->GetVal());
583  oh->SetBit(BIT(20));
584  out->Add(oh);
585  }
586  else
587  AliWarning("low cuts histogram not found in input list");
588  if (highCuts && nFiles) {
589  TH1* oh = static_cast<TH1*>(highCuts->Clone());
590  oh->Scale(1. / nFiles->GetVal());
591  oh->SetBit(BIT(20));
592  out->Add(oh);
593  }
594  else
595  AliWarning("high cuts histogram not found in input list");
596 
597  fRingHistos.Add(new RingHistos(1, 'I'));
598  fRingHistos.Add(new RingHistos(2, 'I'));
599  fRingHistos.Add(new RingHistos(2, 'O'));
600  fRingHistos.Add(new RingHistos(3, 'I'));
601  fRingHistos.Add(new RingHistos(3, 'O'));
602 
603  TIter next(&fRingHistos);
604  RingHistos* o = 0;
605  THStack* sums = new THStack("sums", "Sum of ring signals");
606  THStack* sumsESD = new THStack("sumsESD", "Sum of ring ESD signals");
607  while ((o = static_cast<RingHistos*>(next()))) {
608  o->Terminate(d, nEvents);
609  if (!o->fSum) {
610  Warning("Terminate", "No sum histogram found for ring %s", o->GetName());
611  continue;
612  }
613  TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
614  sum->Scale(1., "width");
615  sum->SetTitle(o->GetName());
616  sum->SetDirectory(0);
617  sum->SetYTitle("#sum_{c} #Delta/#Delta_{mip}");
618  sums->Add(sum);
619 
620 
621  if (o->fSumESD) {
622  sum = o->fSumESD->ProjectionX(o->GetName(), 1, o->fSumESD->GetNbinsY(),"e");
623  sum->Scale(1., "width");
624  sum->SetTitle(o->GetName());
625  sum->SetDirectory(0);
626  sum->SetYTitle("#sum_{s} #Delta/#Delta_{mip}");
627  sumsESD->Add(sum);
628  }
629  }
630  out->Add(sums);
631  out->Add(sumsESD);
632  output->Add(out);
633 }
634 
635 //____________________________________________________________________
636 void
638 {
639  //
640  // Define the output histograms. These are put in a sub list of the
641  // passed list. The histograms are merged before the parent task calls
642  // AliAnalysisTaskSE::Terminate
643  //
644  // Parameters:
645  // dir Directory to add to
646  //
647  DGUARD(fDebug,1, "Define output in AliFMDSharingFilter");
648  TList* d = new TList;
649  d->SetOwner();
650  d->SetName(GetName());
651  dir->Add(d);
652 
653 #if 0
654  fSummed = new TH2I("operations", "Strip operations",
655  kMergedInto, kNone-.5, kMergedInto+.5,
656  51201, -.5, 51200.5);
657  fSummed->SetXTitle("Operation");
658  fSummed->SetYTitle("# of strips");
659  fSummed->SetZTitle("Events");
660  fSummed->GetXaxis()->SetBinLabel(kNone, "None");
661  fSummed->GetXaxis()->SetBinLabel(kCandidate, "Candidate");
662  fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
663  fSummed->GetXaxis()->SetBinLabel(kMergedInto, "Merged into");
664  fSummed->SetDirectory(0);
665  d->Add(fSummed);
666 #endif
667 
668  fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
669  fHighCuts->SetXTitle("#eta");
670  fHighCuts->SetDirectory(0);
671  d->Add(fHighCuts);
672 
673  fLowCuts = new TH2D("lowCuts", "Low cuts used", 1,0,1, 1,0,1);
674  fLowCuts->SetXTitle("#eta");
675  fLowCuts->SetDirectory(0);
676  d->Add(fLowCuts);
677 
678  // d->Add(lowCut);
679  // d->Add(nXi);
680  // d->Add(sigma);
682  d->Add(AliForwardUtil::MakeParameter("lowSignal",
686  d->Add(AliForwardUtil::MakeParameter("disabled", fMergingDisabled));
687  TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
688  nFiles->SetMergeMode('+');
689  d->Add(nFiles);
690 
691  fLCuts.Output(d,"lCuts");
692  fHCuts.Output(d,"hCuts");
693 
694  fRingHistos.Add(new RingHistos(1, 'I'));
695  fRingHistos.Add(new RingHistos(2, 'I'));
696  fRingHistos.Add(new RingHistos(2, 'O'));
697  fRingHistos.Add(new RingHistos(3, 'I'));
698  fRingHistos.Add(new RingHistos(3, 'O'));
699 
700  TIter next(&fRingHistos);
701  RingHistos* o = 0;
702  while ((o = static_cast<RingHistos*>(next()))) {
703  o->CreateOutputObjects(d);
704  }
705 }
706 #define PF(N,V,...) \
707  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
708 #define PFB(N,FLAG) \
709  do { \
710  AliForwardUtil::PrintName(N); \
711  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
712  } while(false)
713 #define PFV(N,VALUE) \
714  do { \
715  AliForwardUtil::PrintName(N); \
716  std::cout << (VALUE) << std::endl; } while(false)
717 
718 //____________________________________________________________________
719 void
721 {
722  //
723  // Print information
724  //
725  // Parameters:
726  // option Not used
727  //
729  gROOT->IncreaseDirLevel();
730 
731  PFB("Use corrected angles", fCorrectAngles);
732  PFB("Zero below threshold", fZeroSharedHitsBelowThreshold);
733  PFB("Use simple sharing", fUseSimpleMerging);
734  PFB("Allow 3 strip merging", fThreeStripSharing);
735  PF("Low cuts", "");
736  fLCuts.Print();
737  PF("High cuts", "");
738  fHCuts.Print();
739  gROOT->DecreaseDirLevel();
740 }
741 
742 //====================================================================
745  fBefore(0),
746  fAfter(0),
747  fSingle(0),
748  fDouble(0),
749  fTriple(0),
750  fSinglePerStrip(0),
751  // fDistanceBefore(0),
752  // fDistanceAfter(0),
753  fBeforeAfter(0),
754  fNeighborsBefore(0),
755  fNeighborsAfter(0),
756  fSumESD(0),
757  fSum(0),
758  fNConsecutive(0)
759  // ,
760  // fHits(0),
761  // fNHits(0)
762 {
763  //
764  // Default CTOR
765  //
766 
767 }
768 
769 //____________________________________________________________________
771  : AliForwardUtil::RingHistos(d,r),
772  fBefore(0),
773  fAfter(0),
774  fSingle(0),
775  fDouble(0),
776  fTriple(0),
777  fSinglePerStrip(0),
778  // fDistanceBefore(0),
779  // fDistanceAfter(0),
780  fBeforeAfter(0),
781  fNeighborsBefore(0),
782  fNeighborsAfter(0),
783  fSumESD(0),
784  fSum(0),
785  fNConsecutive(0)
786  //,
787  // fHits(0),
788  // fNHits(0)
789 {
790  //
791  // Constructor
792  //
793  // Parameters:
794  // d detector
795  // r ring
796  //
797  fBefore = new TH1D("esdEloss", Form("Energy loss in %s (reconstruction)",
798  GetName()), 640, -1, 15);
799  fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
800  fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
801  fBefore->SetFillColor(Color());
802  fBefore->SetFillStyle(3001);
803  fBefore->SetLineColor(kBlack);
804  fBefore->SetLineStyle(2);
805  fBefore->SetDirectory(0);
806 
807  fAfter = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
808  fAfter->SetTitle(Form("Energy loss in %s (sharing corrected)", GetName()));
809  fAfter->SetFillColor(Color()+2);
810  fAfter->SetLineStyle(1);
811  fAfter->SetDirectory(0);
812 
813  fSingle = new TH1D("singleEloss", "Energy loss (single strips)",
814  600, 0, 15);
815  fSingle->SetXTitle("#Delta/#Delta_{mip}");
816  fSingle->SetYTitle("P(#Delta/#Delta_{mip})");
817  fSingle->SetFillColor(Color());
818  fSingle->SetFillStyle(3001);
819  fSingle->SetLineColor(kBlack);
820  fSingle->SetLineStyle(2);
821  fSingle->SetDirectory(0);
822 
823  fDouble = static_cast<TH1D*>(fSingle->Clone("doubleEloss"));
824  fDouble->SetTitle("Energy loss (two strips)");
825  fDouble->SetFillColor(Color()+1);
826  fDouble->SetDirectory(0);
827 
828  fTriple = static_cast<TH1D*>(fSingle->Clone("tripleEloss"));
829  fTriple->SetTitle("Energy loss (three strips)");
830  fTriple->SetFillColor(Color()+2);
831  fTriple->SetDirectory(0);
832 
833  //Int_t nBinsForInner = (r == 'I' ? 32 : 16);
834  Int_t nBinsForInner = (r == 'I' ? 512 : 256);
835  Int_t nStrips = (r == 'I' ? 512 : 256);
836 
837  fSinglePerStrip = new TH2D("singlePerStrip", "SinglePerStrip",
838  600,0,15, nBinsForInner,0,nStrips);
839  fSinglePerStrip->SetXTitle("#Delta/#Delta_{mip}");
840  fSinglePerStrip->SetYTitle("Strip #");
841  fSinglePerStrip->SetZTitle("Counts");
842  fSinglePerStrip->SetDirectory(0);
843 
844 #if 0
845  fDistanceBefore = new TH1D("distanceBefore", "Distance before sharing",
846  nStrips , 0,nStrips );
847  fDistanceBefore->SetXTitle("Distance");
848  fDistanceBefore->SetYTitle("Counts");
849  fDistanceBefore->SetFillColor(kGreen+2);
850  fDistanceBefore->SetFillStyle(3001);
851  fDistanceBefore->SetLineColor(kBlack);
852  fDistanceBefore->SetLineStyle(2);
853  fDistanceBefore->SetDirectory(0);
854 
855  fDistanceAfter = static_cast<TH1D*>(fDistanceBefore->Clone("distanceAfter"));
856  fDistanceAfter->SetTitle("Distance after sharing");
857  fDistanceAfter->SetFillColor(kGreen+1);
858  fDistanceAfter->SetDirectory(0);
859 #endif
860 
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 before");
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  fSumESD = new TH2D("summedESD", "Summed ESD signal", 200, -4, 6,
883  NSector(), 0, 2*TMath::Pi());
884  fSumESD->SetDirectory(0);
885  fSumESD->Sumw2();
886  fSumESD->SetMarkerColor(Color());
887  // fSum->SetFillColor(Color());
888  fSumESD->SetXTitle("#eta");
889  fSumESD->SetYTitle("#varphi [radians]");
890  fSumESD->SetZTitle("#sum_{strip} #Delta/#Delta_{mip}(#eta,#varphi) ");
891 
892  fSum = static_cast<TH2D*>(fSumESD->Clone("summed"));
893  fSum->SetTitle("Summed cluster signal");
894  fSum->SetZTitle("#sum_{cluster} #Delta/#Delta_{mip}(#eta,#varphi) ");
895  fSum->SetDirectory(0);
896 
897  // Perhaps we need to ensure that this histogram has enough range to
898  // accommondate all possible ranges - that is, from -1 to the number
899  // of strips in this ring(-type) - i.e., NStrips(). Perhaps the
900  // axis should be defined with increasin bin size - e.g.,
901  //
902  // -1.5,-.5,.5,1.5,...,100.5,128.5,192.5,...,NStrips()
903  //
904  fNConsecutive = new TH1D("nConsecutive","# consecutive strips above low cut",
905  201,-1.5,199.5);
906  fNConsecutive->SetXTitle("N_{strips}");
907  fNConsecutive->SetYTitle("N_{entries}");
908  fNConsecutive->SetFillColor(kYellow+2);
909  fNConsecutive->SetFillStyle(3001);
910  fNConsecutive->SetDirectory(0);
911 
912 
913 #if 0
914  fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
915  fHits->SetDirectory(0);
916  fHits->SetXTitle("# of hits");
917  fHits->SetFillColor(kGreen+1);
918 #endif
919 }
920 //____________________________________________________________________
922  : AliForwardUtil::RingHistos(o),
923  fBefore(o.fBefore),
924  fAfter(o.fAfter),
925  fSingle(o.fSingle),
926  fDouble(o.fDouble),
927  fTriple(o.fTriple),
929  // fDistanceBefore(o.fDistanceBefore),
930  // fDistanceAfter(o.fDistanceAfter),
934  fSumESD(o.fSumESD), //,
935  fSum(o.fSum),
937  //,
938  // fHits(o.fHits),
939  // fNHits(o.fNHits)
940 {
941  //
942  // Copy constructor
943  //
944  // Parameters:
945  // o Object to copy from
946  //
947 }
948 
949 //____________________________________________________________________
952 {
953  //
954  // Assignment operator
955  //
956  // Parameters:
957  // o Object to assign from
958  //
959  // Return:
960  // Reference to this
961  //
962  if (&o == this) return *this;
964  fDet = o.fDet;
965  fRing = o.fRing;
966 
967  if (fBefore) delete fBefore;
968  if (fAfter) delete fAfter;
969  if (fSingle) delete fSingle;
970  if (fDouble) delete fDouble;
971  if (fTriple) delete fTriple;
972  if (fSinglePerStrip) delete fSinglePerStrip;
973  if (fNConsecutive) delete fNConsecutive;
974  // if (fDistanceBefore) delete fDistanceBefore;
975  // if (fDistanceAfter) delete fDistanceAfter;
976  // if (fHits) delete fHits;
977 
978 
979  fBefore = static_cast<TH1D*>(o.fBefore->Clone());
980  fAfter = static_cast<TH1D*>(o.fAfter->Clone());
981  fSingle = static_cast<TH1D*>(o.fSingle->Clone());
982  fDouble = static_cast<TH1D*>(o.fDouble->Clone());
983  fTriple = static_cast<TH1D*>(o.fTriple->Clone());
984  fSinglePerStrip = static_cast<TH2D*>(o.fSinglePerStrip->Clone());
985  // fDistanceBefore = static_cast<TH1D*>(o.fDistanceBefore->Clone());
986  // fDistanceAfter = static_cast<TH1D*>(o.fDistanceAfter->Clone());
987  fBeforeAfter = static_cast<TH2D*>(o.fBeforeAfter->Clone());
988  fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
989  fNeighborsAfter = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
990  // fHits = static_cast<TH1D*>(o.fHits->Clone());
991  fSumESD = static_cast<TH2D*>(o.fSumESD->Clone());
992  fSum = static_cast<TH2D*>(o.fSum->Clone());
993  fNConsecutive = static_cast<TH1D*>(o.fNConsecutive->Clone());
994 
995  return *this;
996 }
997 //____________________________________________________________________
999 {
1000  //
1001  // Destructor
1002  //
1003 }
1004 #if 0
1005 //____________________________________________________________________
1006 void
1007 AliFMDSharingFilter::RingHistos::Finish()
1008 {
1009  //
1010  // Finish off
1011  //
1012  //
1013  // fHits->Fill(fNHits);
1014 }
1015 #endif
1016 //____________________________________________________________________
1017 void
1019 {
1020  //
1021  // Scale the histograms to the total number of events
1022  //
1023  // Parameters:
1024  // nEvents Number of events
1025  // dir Where the output is
1026  //
1027  TList* l = GetOutputList(dir);
1028  if (!l) return;
1029 
1030 #if 0
1031  TH1D* before = static_cast<TH1D*>(l->FindObject("esdEloss"));
1032  TH1D* after = static_cast<TH1D*>(l->FindObject("anaEloss"));
1033  if (before) before->Scale(1./nEvents);
1034  if (after) after->Scale(1./nEvents);
1035 #endif
1036 
1037  TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
1038  if (summed) summed->Scale(1./nEvents);
1039  fSum = summed;
1040 
1041  TH2D* summedESD = static_cast<TH2D*>(l->FindObject("summedESD"));
1042  if (summedESD) summedESD->Scale(1./nEvents);
1043  fSumESD = summedESD;
1044 
1045  TH1D* consecutive = static_cast<TH1D*>(l->FindObject("nConsecutive"));
1046  if (consecutive) consecutive->Scale(1./nEvents);
1047  fNConsecutive= consecutive;
1048 }
1049 
1050 //____________________________________________________________________
1051 void
1053 {
1054  //
1055  // Make output
1056  //
1057  // Parameters:
1058  // dir where to store
1059  //
1060  TList* d = DefineOutputList(dir);
1061 
1062  d->Add(fBefore);
1063  d->Add(fAfter);
1064  d->Add(fSingle);
1065  d->Add(fDouble);
1066  d->Add(fTriple);
1067  d->Add(fSinglePerStrip);
1068  // d->Add(fDistanceBefore);
1069  // d->Add(fDistanceAfter);
1070  d->Add(fBeforeAfter);
1071  d->Add(fNeighborsBefore);
1072  d->Add(fNeighborsAfter);
1073  // d->Add(fHits);
1074  d->Add(fSumESD);
1075  d->Add(fSum);
1076  d->Add(fNConsecutive);
1077 
1078  // Removed to avoid doubly adding the list which destroys
1079  // the merging
1080  //dir->Add(d);
1081 }
1082 
1083 //____________________________________________________________________
1084 //
1085 // EOF
1086 //
RingHistos & operator=(const RingHistos &o)
void Print(Option_t *option="") const
void Terminate(const TList *dir, Int_t nEvents)
double Double_t
Definition: External.C:58
Double_t SignalInStrip(const AliESDFMD &fmd, UShort_t d, Char_t r, UShort_t s, UShort_t t) const
Double_t AngleCorrect(Double_t mult, Double_t eta) const
const char * title
Definition: MakeQAPdf.C:27
const TAxis & GetEtaAxis() const
char Char_t
Definition: External.C:18
virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta, Bool_t errors=true) const
void SetupForData(const TAxis &axis)
#define DMSG(L, N, F,...)
const char * GetName() const
#define PF(N, V,...)
virtual void Terminate(const TList *dir, TList *output, Int_t nEvents)
int Int_t
Definition: External.C:63
void Output(TList *l, const char *name=0) const
float Float_t
Definition: External.C:68
Bool_t Filter(const AliESDFMD &input, Bool_t lowFlux, AliESDFMD &output, Double_t zvtx)
Definition: External.C:228
Definition: External.C:212
void Set(EMethod method, Double_t fmd1i, Double_t fmd2i=-1, Double_t fmd2o=-1, Double_t fmd3i=-1, Double_t fmd3o=-1)
Double_t nEvents
plot quality messages
const AliFMDCorrELossFit * GetELossFit() const
#define DGUARD(L, N, F,...)
static void PrintTask(const TObject &o)
virtual void CreateOutputObjects(TList *dir)
static TObject * MakeParameter(const char *name, UShort_t value)
const UShort_t & NSector() const
Double_t DeAngleCorrect(Double_t mult, Double_t eta) const
Definition: External.C:220
virtual void Print(Option_t *option="") const
void FillHistogram(TH2 *h) const
virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const
RingHistos * GetRingHistos(UShort_t d, Char_t r) const
unsigned short UShort_t
Definition: External.C:28
TList * DefineOutputList(TList *d) const
const char Option_t
Definition: External.C:48
RingHistos & operator=(const RingHistos &o)
#define PFB(N, FLAG)
bool Bool_t
Definition: External.C:53
Definition: External.C:196
static AliForwardCorrectionManager & Instance()
TList * GetOutputList(const TList *d) const
TDirectoryFile * dir