More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDSharingFilter.cxx
index 26176494dd34d877b201aac3e327dee1b5d74158..ac6fc5630b9c1310d60398fe893f05a08ccc8697 100644 (file)
@@ -7,8 +7,12 @@
 #include <TList.h>
 #include <TH1.h>
 #include <TMath.h>
-#include "AliFMDAnaParameters.h"
+#include "AliForwardCorrectionManager.h"
+// #include "AliFMDAnaParameters.h"
 #include <AliLog.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
 
 ClassImp(AliFMDSharingFilter)
 #if 0
@@ -20,28 +24,26 @@ ClassImp(AliFMDSharingFilter)
 AliFMDSharingFilter::AliFMDSharingFilter()
   : TNamed(), 
     fRingHistos(),
-    fLowCut(0.3),
-    fCorrectAngles(kFALSE),
-    fEtaCorr(0)
+    fLowCut(0.),
+    fCorrectAngles(kFALSE), 
+    fNXi(1),
+    fDebug(0)
 {}
 
 //____________________________________________________________________
 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
   : TNamed("fmdSharingFilter", title), 
     fRingHistos(), 
-    fLowCut(0.3),
-    fCorrectAngles(kFALSE),
-    fEtaCorr(0)
+    fLowCut(0.),
+    fCorrectAngles(kFALSE), 
+    fNXi(1),
+    fDebug(0)
 {
   fRingHistos.Add(new RingHistos(1, 'I'));
   fRingHistos.Add(new RingHistos(2, 'I'));
   fRingHistos.Add(new RingHistos(2, 'O'));
   fRingHistos.Add(new RingHistos(3, 'I'));
   fRingHistos.Add(new RingHistos(3, 'O'));
-  fEtaCorr = new TH2F("etaCorr", "Correction of #eta", 
-                     200, -4, 6,  200, -4, 6);
-  fEtaCorr->SetXTitle("#eta from ESD");
-  fEtaCorr->SetYTitle("#eta from AnaParameters");
 }
 
 //____________________________________________________________________
@@ -49,8 +51,9 @@ AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
   : TNamed(o), 
     fRingHistos(), 
     fLowCut(o.fLowCut),
-    fCorrectAngles(o.fCorrectAngles),
-    fEtaCorr(o.fEtaCorr)
+    fCorrectAngles(o.fCorrectAngles), 
+    fNXi(o.fNXi),
+    fDebug(o.fDebug)
 {
   TIter    next(&o.fRingHistos);
   TObject* obj = 0;
@@ -67,11 +70,12 @@ AliFMDSharingFilter::~AliFMDSharingFilter()
 AliFMDSharingFilter&
 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
 {
-  SetName(o.GetName());
-  SetTitle(o.GetTitle());
+  TNamed::operator=(o);
 
   fLowCut        = o.fLowCut;
   fCorrectAngles = o.fCorrectAngles;
+  fNXi           = o.fNXi;
+  fDebug         = o.fDebug;
 
   fRingHistos.Delete();
   TIter    next(&o.fRingHistos);
@@ -100,16 +104,16 @@ AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
 Bool_t
 AliFMDSharingFilter::Filter(const AliESDFMD& input, 
                            Bool_t           lowFlux,
-                           AliESDFMD&       output,
-                           Double_t         vz)
+                           AliESDFMD&       output)
 {
   output.Clear();
   TIter    next(&fRingHistos);
-  RingHistos* o = 0;
+  RingHistos* o      = 0;
+  Double_t    lowCut = GetLowCut();
   while ((o = static_cast<RingHistos*>(next())))
     o->Clear();
 
-  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 
   for(UShort_t d = 1; d <= 3; d++) {
     Int_t nRings = (d == 1 ? 1 : 2);
@@ -136,10 +140,7 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input,
          if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 
          // Get the pseudo-rapidity 
-         Double_t eta1 = input.Eta(d,r,s,t);
-         Double_t eta2 = pars->GetEtaFromStrip(d,r,s,t,vz);
-         fEtaCorr->Fill(eta1, eta2);
-         Double_t eta = eta2;
+         Double_t eta = input.Eta(d,r,s,t);
          
          // Fill the diagnostics histogram 
          histos->fBefore->Fill(mult);
@@ -159,7 +160,7 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input,
          Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
                                                      lowFlux,d,r,s,t, 
                                                      usedPrev,usedThis);
-         if (mergedEnergy > fLowCut) histos->Incr();
+         if (mergedEnergy > lowCut) histos->Incr();
          histos->fAfter->Fill(mergedEnergy);
 
          output.SetMultiplicity(d,r,s,t,mergedEnergy);
@@ -193,21 +194,41 @@ AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
     mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
   return mult;
 }
+//_____________________________________________________________________
+Double_t 
+AliFMDSharingFilter::GetLowCut() const
+{
+  if (fLowCut > 0) return fLowCut;
+  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
+  AliFMDCorrELossFit* fits = fcm.GetELossFit();
+  return fits->GetLowCut();
+}
                        
 //_____________________________________________________________________
 Double_t 
 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
 {
-  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
+
  
   // Get the high cut.  The high cut is defined as the 
   // most-probably-value peak found from the energy distributions, minus 
   // 2 times the width of the corresponding Landau.
-  Double_t mpv     = pars->GetMPV(d,r,eta);
-  Double_t w       = pars->GetSigma(d,r,eta);
-  if (mpv > 100) 
-    AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
-  Double_t highCut = (mpv - 2 * w);
+  AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
+  Double_t delta = 1024;
+  Double_t xi    = 1024;
+  if (!fit) {
+    AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
+  }
+  else {
+    delta = fit->fDelta;
+    xi    = fit->fXi;
+  }
+
+  if (delta > 100) {
+    AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
+  }
+  Double_t highCut = (delta - fNXi * xi);
   return highCut;
 }
           
@@ -227,7 +248,8 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 {
   // IF the multiplicity is very large, or below the cut, reject it, and 
   // flag it as candidate for sharing 
-  if(mult > 12 || mult < fLowCut) {
+  Double_t    lowCut = GetLowCut();
+  if(mult > 12 || mult < lowCut) {
     usedThis      = kFALSE;
     usedPrev      = kFALSE;
     return 0;
@@ -243,7 +265,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 
   //analyse and perform sharing on one strip
   
-  // Calculate the total energy loss 
+  // Get the high cut 
   Double_t highCut = GetHighCut(d, r, eta);
 
   // If this signal is smaller than the next, and the next signal is smaller 
@@ -265,7 +287,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
   // If the previous signal was larger than the low cut, and 
   // the previous was smaller than high cut, and the previous signal 
   // isn't marked as used, then add it's energy to this signal 
-  if (prevE > fLowCut && 
+  if (prevE > lowCut && 
       prevE < highCut && 
       !usedPrev) 
     totalE += prevE;
@@ -273,7 +295,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
   // If the next signal is larger than the low cut, and 
   // the next signal is smaller than the low cut, then add the next signal 
   // to this, and marked the next signal as used 
-  if(nextE > fLowCut && nextE < highCut ) {
+  if(nextE > lowCut && nextE < highCut ) {
     totalE      += nextE;
     usedThis =  kTRUE;
   }
@@ -324,37 +346,48 @@ AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
 
 //____________________________________________________________________
 void
-AliFMDSharingFilter::ScaleHistograms(Int_t nEvents)
+AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
 {
   if (nEvents <= 0) return;
+  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
+  if (!d) return;
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
-  while ((o = static_cast<RingHistos*>(next()))) {
-    o->fBefore->Scale(1. / nEvents);
-    o->fAfter->Scale(1. / nEvents);
-  }
+  while ((o = static_cast<RingHistos*>(next())))
+    o->ScaleHistograms(d, nEvents);
 }
 
 //____________________________________________________________________
 void
-AliFMDSharingFilter::Output(TList* dir)
+AliFMDSharingFilter::DefineOutput(TList* dir)
 {
   TList* d = new TList;
   d->SetName(GetName());
   dir->Add(d);
-  d->Add(fEtaCorr);
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next()))) {
     o->Output(d);
   }
 }
+//____________________________________________________________________
+void
+AliFMDSharingFilter::Print(Option_t* /*option*/) const
+{
+  char ind[gROOT->GetDirLevel()+1];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  std::cout << ind << "AliFMDSharingFilter: " << GetName() << '\n'
+           << ind << " Low cut:                " << fLowCut << '\n'
+           << ind << " N xi factor:            " << fNXi    << '\n'
+           << ind << " Use corrected angles:   " 
+           << (fCorrectAngles ? "yes" : "no") << std::endl;
+}
   
 //====================================================================
 AliFMDSharingFilter::RingHistos::RingHistos()
-  : fDet(0),
-    fRing('\0'),
+  : AliForwardUtil::RingHistos(), 
     fBefore(0), 
     fAfter(0), 
     fHits(0),
@@ -363,19 +396,18 @@ AliFMDSharingFilter::RingHistos::RingHistos()
 
 //____________________________________________________________________
 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
-  : fDet(d), 
-    fRing(r),
+  : AliForwardUtil::RingHistos(d,r), 
     fBefore(0), 
     fAfter(0),
     fHits(0),
     fNHits(0)
 {
-  fBefore = new TH1D(Form("FMD%d%c_ESD_Eloss", d, r), 
-                    Form("Energy loss in FMD%d%c (reconstruction)", d, r), 
-                    1000, 0, 25);
-  fAfter  = new TH1D(Form("FMD%d%c_Ana_Eloss", d, r), 
-                    Form("Energy loss in FMD%d%c (sharing corrected)", d, r), 
+  fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()), 
+                    Form("Energy loss in %s (reconstruction)", fName.Data()), 
                     1000, 0, 25);
+  fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()), 
+                    Form("Energy loss in %s (sharing corrected)",
+                         fName.Data()), 1000, 0, 25);
   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
   fBefore->SetFillColor(kRed+1);
@@ -389,8 +421,8 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
   // fAfter->Sumw2();
   fAfter->SetDirectory(0);
 
-  fHits = new TH1D(Form("FMD%d%c_Hits", d, r), 
-                  Form("Number of hits in FMD%d%c", d, r), 
+  fHits = new TH1D(Form("%s_Hits", fName.Data()), 
+                  Form("Number of hits in %s", fName.Data()), 
                   200, 0, 200000);
   fHits->SetDirectory(0);
   fHits->SetXTitle("# of hits");
@@ -398,9 +430,7 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
 }
 //____________________________________________________________________
 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
-  : TObject(o), 
-    fDet(o.fDet), 
-    fRing(o.fRing), 
+  : AliForwardUtil::RingHistos(o), 
     fBefore(o.fBefore), 
     fAfter(o.fAfter),
     fHits(o.fHits),
@@ -411,6 +441,7 @@ AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
 AliFMDSharingFilter::RingHistos&
 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
 {
+  AliForwardUtil::RingHistos::operator=(o);
   fDet = o.fDet;
   fRing = o.fRing;
   
@@ -438,12 +469,26 @@ AliFMDSharingFilter::RingHistos::Finish()
   fHits->Fill(fNHits);
 }
 
+//____________________________________________________________________
+void
+AliFMDSharingFilter::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
+{
+  TList* l = GetOutputList(dir);
+  if (!l) return; 
+
+  TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
+                                                      fName.Data())));
+  TH1D* after  = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
+                                                      fName.Data())));
+  if (before) before->Scale(1./nEvents);
+  if (after)  after->Scale(1./nEvents);
+}
+
 //____________________________________________________________________
 void
 AliFMDSharingFilter::RingHistos::Output(TList* dir)
 {
-  TList* d = new TList;
-  d->SetName(Form("FMD%d%c", fDet, fRing)); 
+  TList* d = DefineOutputList(dir);
   d->Add(fBefore);
   d->Add(fAfter);
   d->Add(fHits);