[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDSharingFilter.cxx

//
// Class to do the sharing correction.  That is, a filter that merges 
// adjacent strip signals presumably originating from a single particle 
// that impinges on the detector in such a way that it deposite energy 
// into two or more strips. 
//
// Input: 
//    - AliESDFMD object  - from reconstruction
//
// Output: 
//    - AliESDFMD object  - copy of input, but with signals merged 
//
// Corrections used: 
//    - AliFMDCorrELossFit
//
// Histograms: 
//    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of 
//      signals before and after the filter.  
//    - For each ring (see above), an array of distributions of number of 
//      hit strips for each vertex bin (if enabled - see Init method)
// 
//
//
#include "AliFMDSharingFilter.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TH1.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
#include "AliFMDCorrELossFit.h"
#include <AliLog.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

ClassImp(AliFMDSharingFilter)
#if 0
; // This is for Emacs
#endif 


//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter()
  : TNamed(), 
    fRingHistos(),
    fLowCut(0.),
    fCorrectAngles(kFALSE), 
    fNXi(1),
    fDebug(0)
{
  // 
  // Default Constructor - do not use 
  //
}

//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
  : TNamed("fmdSharingFilter", title), 
    fRingHistos(), 
    fLowCut(0.),
    fCorrectAngles(kFALSE), 
    fNXi(1),
    fDebug(0)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    title Title of object  - not significant 
  //
  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'));
}

//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
  : TNamed(o), 
    fRingHistos(), 
    fLowCut(o.fLowCut),
    fCorrectAngles(o.fCorrectAngles), 
    fNXi(o.fNXi),
    fDebug(o.fDebug)
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

//____________________________________________________________________
AliFMDSharingFilter::~AliFMDSharingFilter()
{
  // 
  // Destructor
  //
  fRingHistos.Delete();
}

//____________________________________________________________________
AliFMDSharingFilter&
AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this 
  //
  TNamed::operator=(o);

  fLowCut        = o.fLowCut;
  fCorrectAngles = o.fCorrectAngles;
  fNXi           = o.fNXi;
  fDebug         = o.fDebug;

  fRingHistos.Delete();
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
  
  return *this;
}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos*
AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
{
  // 
  // Get the ring histogram container 
  // 
  // Parameters:
  //    d Detector
  //    r Ring 
  // 
  // Return:
  //    Ring histogram container 
  //
  Int_t idx = -1;
  switch (d) { 
  case 1: idx = 0; break;
  case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  }
  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
  
  return static_cast<RingHistos*>(fRingHistos.At(idx));
}

//____________________________________________________________________
Bool_t
AliFMDSharingFilter::Filter(const AliESDFMD& input, 
			    Bool_t           lowFlux,
			    AliESDFMD&       output)
{
  // 
  // Filter the input AliESDFMD object
  // 
  // Parameters:
  //    input     Input 
  //    lowFlux   If this is a low-flux event 
  //    output    Output AliESDFMD object 
  // 
  // Return:
  //    True on success, false otherwise 
  //
  output.Clear();
  TIter    next(&fRingHistos);
  RingHistos* o      = 0;
  Double_t    lowCut = GetLowCut();
  while ((o = static_cast<RingHistos*>(next())))
    o->Clear();

  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();

  for(UShort_t d = 1; d <= 3; d++) {
    Int_t nRings = (d == 1 ? 1 : 2);
    for (UShort_t q = 0; q < nRings; q++) {
      Char_t      r    = (q == 0 ? 'I' : 'O');
      UShort_t    nsec = (q == 0 ?  20 :  40);
      UShort_t    nstr = (q == 0 ? 512 : 256);

      RingHistos* histos = GetRingHistos(d, r);
      
      for(UShort_t s = 0; s < nsec;  s++) {
	Bool_t usedThis = kFALSE;
	Bool_t usedPrev = kFALSE;
	
	for(UShort_t t = 0; t < nstr; t++) {
	  output.SetMultiplicity(d,r,s,t,0.);
	  Float_t mult = SignalInStrip(input,d,r,s,t);
	  
	  // Keep dead-channel information. 
	  if(mult == AliESDFMD::kInvalidMult)
	    output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);

	  // If no signal or dead strip, go on. 
	  if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;

	  // Get the pseudo-rapidity 
	  Double_t eta = input.Eta(d,r,s,t);
	  
	  // Fill the diagnostics histogram 
	  histos->fBefore->Fill(mult);

	  // Get next and previous signal - if any 
	  Double_t prevE = 0;
	  Double_t nextE = 0;
	  if (t != 0) {
	    prevE = SignalInStrip(input,d,r,s,t-1);
	    if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
	  }
	  if (t != nstr - 1) {
	    nextE = SignalInStrip(input,d,r,s,t+1);
	    if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
	  }
	  
	  Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
						      lowFlux,d,r,s,t, 
						      usedPrev,usedThis);
	  if (mergedEnergy > lowCut) histos->Incr();
	  histos->fAfter->Fill(mergedEnergy);

	  output.SetMultiplicity(d,r,s,t,mergedEnergy);
	  output.SetEta(d,r,s,t,eta);
	} // for strip
      } // for sector
    } // for ring 
  } // for detector

  next.Reset();
  while ((o = static_cast<RingHistos*>(next())))
    o->Finish();

  return kTRUE;
}

//_____________________________________________________________________
Double_t 
AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input, 
				   UShort_t         d,
				   Char_t           r,
				   UShort_t         s,
				   UShort_t         t) const
{
  // 
  // Get the signal in a strip 
  // 
  // Parameters:
  //    fmd   ESD object
  //    d     Detector
  //    r     Ring 
  //    s     Sector 
  //    t     Strip
  // 
  // Return:
  //    The energy signal 
  //
  Double_t mult = input.Multiplicity(d,r,s,t);
  if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;

  if (fCorrectAngles && !input.IsAngleCorrected()) 
    mult = AngleCorrect(mult, input.Eta(d,r,s,t));
  else if (!fCorrectAngles && input.IsAngleCorrected()) 
    mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
  return mult;
}
//_____________________________________________________________________
Double_t 
AliFMDSharingFilter::GetLowCut() const
{
  //
  // Get the low cut.  Normally, the low cut is taken to be the lower
  // value of the fit range used when generating the energy loss fits.
  // However, if fLowCut is set (using SetLowCit) to a value greater
  // than 0, then that value is used.
  //
  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
{
  //
  // 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.
  //
  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.
  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;
}
	   
//_____________________________________________________________________
Double_t 
AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
					 Double_t eta,
					 Double_t prevE,
					 Double_t nextE,
					 Bool_t   lowFlux,
					 UShort_t d,
					 Char_t   r,
					 UShort_t /*s*/,
					 UShort_t /*t*/,
					 Bool_t&  usedPrev, 
					 Bool_t&  usedThis) const
{
  // 
  // The actual algorithm 
  // 
  // Parameters:
  //    mult      The unfiltered signal in the strip
  //    eta       Psuedo rapidity 
  //    prevE     Previous strip signal (or 0)
  //    nextE     Next strip signal (or 0) 
  //    lowFlux   Whether this is a low flux event 
  //    d         Detector
  //    r         Ring 
  //    s         Sector 
  //    t         Strip
  //    usedPrev  Whether the previous strip was used in sharing or not
  //    usedThis  Wether this strip was used in sharing or not. 
  // 
  // Return:
  //    The filtered signal in the strip
  //

  // IF the multiplicity is very large, or below the cut, reject it, and 
  // flag it as candidate for sharing 
  Double_t    lowCut = GetLowCut();
  if(mult > 12 || mult < lowCut) {
    usedThis      = kFALSE;
    usedPrev      = kFALSE;
    return 0;
  }

  // If this was shared with the previous signal, return 0 and mark it as 
  // not a candiate for sharing 
  if(usedThis) {
    usedThis      = kFALSE;
    usedPrev      = kTRUE;
    return 0.;
  }

  //analyse and perform sharing on one strip
  
  // 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 
  // than then the high cut, and it's a low-flux event, then mark this and 
  // the next signal as candidates for sharing 
  // 
  // This is the test if the signal is the smaller of two possibly
  // shared signals.   
  if (mult  < nextE   && 
      nextE > highCut && 
      lowFlux) {
    usedThis      = kFALSE;
    usedPrev      = kFALSE;
    return 0;
  }
  
  Double_t totalE  = 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 > lowCut && 
      prevE < highCut && 
      !usedPrev) 
    totalE += prevE;

  // 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 > lowCut && nextE < highCut ) {
    totalE      += nextE;
    usedThis =  kTRUE;
  }

  // If we're processing on non-angle corrected data, we should do the 
  // angle correction here 
  if (!fCorrectAngles) 
    totalE = AngleCorrect(totalE, eta);

  // Fill post processing histogram
  // if(totalE > fLowCut)
  //   GetRingHistos(d,r)->fAfter->Fill(totalE);

  // Return value 
  Double_t mergedEnergy = 0;
  
  if(totalE > 0) {
    // If we have a signal, then this is used (sharedPrev=true) and
    // the signal is set to the result
    mergedEnergy = totalE;
    usedPrev     = kTRUE;
  }
  else {
    // If we do not have a signal (too low), then this is not shared, 
    // and the next strip is not shared either 
    usedThis  = kFALSE;
    usedPrev  = kFALSE;
  }
  
  return mergedEnergy; 
}
//____________________________________________________________________
Double_t
AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
{
  // 
  // Angle correct the signal 
  // 
  // Parameters:
  //    mult Angle Un-corrected Signal 
  //    eta  Pseudo-rapidity 
  // 
  // Return:
  //    Angle corrected signal 
  //
  Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
  if (eta < 0) theta -= TMath::Pi();
  return mult * TMath::Cos(theta);
}
//____________________________________________________________________
Double_t
AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
{
  // 
  // Angle de-correct the signal 
  // 
  // Parameters:
  //    mult Angle corrected Signal 
  //    eta  Pseudo-rapidity 
  // 
  // Return:
  //    Angle un-corrected signal 
  //
  Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
  if (eta < 0) theta -= TMath::Pi();
  return mult / TMath::Cos(theta);
}

//____________________________________________________________________
void
AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
{
  // 
  // Scale the histograms to the total number of events 
  // 
  // Parameters:
  //    dir     Where the output is 
  //    nEvents Number of events 
  //
  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->ScaleHistograms(d, nEvents);
}

//____________________________________________________________________
void
AliFMDSharingFilter::DefineOutput(TList* dir)
{
  // 
  // Define the output histograms.  These are put in a sub list of the
  // passed list.   The histograms are merged before the parent task calls 
  // AliAnalysisTaskSE::Terminate 
  // 
  // Parameters:
  //    dir Directory to add to 
  //
  TList* d = new TList;
  d->SetName(GetName());
  dir->Add(d);
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->Output(d);
  }
}
//____________________________________________________________________
void
AliFMDSharingFilter::Print(Option_t* /*option*/) const
{
  // 
  // Print information
  // 
  // Parameters:
  //    option Not used 
  //
  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()
  : AliForwardUtil::RingHistos(), 
    fBefore(0), 
    fAfter(0), 
    fHits(0),
    fNHits(0)
{
  // 
  // Default CTOR
  //

}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r), 
    fBefore(0), 
    fAfter(0),
    fHits(0),
    fNHits(0)
{
  // 
  // Constructor
  // 
  // Parameters:
  //    d detector
  //    r ring 
  //
  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);
  fBefore->SetFillStyle(3001);
  // fBefore->Sumw2();
  fBefore->SetDirectory(0);
  fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
  fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
  fAfter->SetFillColor(kBlue+1);
  fAfter->SetFillStyle(3001);
  // fAfter->Sumw2();
  fAfter->SetDirectory(0);

  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");
  fHits->SetFillColor(kGreen+1);
}
//____________________________________________________________________
AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
  : AliForwardUtil::RingHistos(o), 
    fBefore(o.fBefore), 
    fAfter(o.fAfter),
    fHits(o.fHits),
    fNHits(o.fNHits)
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos&
AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this 
  //
  AliForwardUtil::RingHistos::operator=(o);
  fDet = o.fDet;
  fRing = o.fRing;
  
  if (fBefore) delete  fBefore;
  if (fAfter)  delete  fAfter;
  if (fHits)   delete fHits;
  
  fBefore = static_cast<TH1D*>(o.fBefore->Clone());
  fAfter  = static_cast<TH1D*>(o.fAfter->Clone());
  fHits  = static_cast<TH1D*>(o.fHits->Clone());
  
  return *this;
}
//____________________________________________________________________
AliFMDSharingFilter::RingHistos::~RingHistos()
{
  // 
  // Destructor 
  //
  if (fBefore) delete fBefore;
  if (fAfter)  delete fAfter;
  if (fHits)   delete fHits;
}
//____________________________________________________________________
void
AliFMDSharingFilter::RingHistos::Finish()
{
  // 
  // Finish off 
  // 
  //
  fHits->Fill(fNHits);
}

//____________________________________________________________________
void
AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir, 
						 Int_t nEvents)
{
  // 
  // Scale the histograms to the total number of events 
  // 
  // Parameters:
  //    nEvents Number of events 
  //    dir     Where the output is 
  //
  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)
{
  // 
  // Make output 
  // 
  // Parameters:
  //    dir where to store 
  //
  TList* d = DefineOutputList(dir);
  d->Add(fBefore);
  d->Add(fAfter);
  d->Add(fHits);
  dir->Add(d);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
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>
	34	#include <iostream>
	35	#include <iomanip>
7e4038b5	36
	37	ClassImp(AliFMDSharingFilter)
	38	#if 0
	39	; // This is for Emacs
	40	#endif
	41
	42
	43	//____________________________________________________________________
	44	AliFMDSharingFilter::AliFMDSharingFilter()
	45	: TNamed(),
	46	fRingHistos(),
0bd4b00f	47	fLowCut(0.),
ea3e5d95	48	fCorrectAngles(kFALSE),
0bd4b00f	49	fNXi(1),
ea3e5d95	50	fDebug(0)
7984e5f7	51	{
	52	//
	53	// Default Constructor - do not use
	54	//
	55	}
7e4038b5	56
	57	//____________________________________________________________________
	58	AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
	59	: TNamed("fmdSharingFilter", title),
	60	fRingHistos(),
0bd4b00f	61	fLowCut(0.),
ea3e5d95	62	fCorrectAngles(kFALSE),
0bd4b00f	63	fNXi(1),
ea3e5d95	64	fDebug(0)
7e4038b5	65	{
7984e5f7	66	//
	67	// Constructor
	68	//
	69	// Parameters:
	70	// title Title of object - not significant
	71	//
7e4038b5	72	fRingHistos.Add(new RingHistos(1, 'I'));
	73	fRingHistos.Add(new RingHistos(2, 'I'));
	74	fRingHistos.Add(new RingHistos(2, 'O'));
	75	fRingHistos.Add(new RingHistos(3, 'I'));
	76	fRingHistos.Add(new RingHistos(3, 'O'));
7e4038b5	77	}
	78
	79	//____________________________________________________________________
	80	AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
	81	: TNamed(o),
	82	fRingHistos(),
	83	fLowCut(o.fLowCut),
ea3e5d95	84	fCorrectAngles(o.fCorrectAngles),
0bd4b00f	85	fNXi(o.fNXi),
ea3e5d95	86	fDebug(o.fDebug)
7e4038b5	87	{
7984e5f7	88	//
	89	// Copy constructor
	90	//
	91	// Parameters:
	92	// o Object to copy from
	93	//
7e4038b5	94	TIter next(&o.fRingHistos);
	95	TObject* obj = 0;
	96	while ((obj = next())) fRingHistos.Add(obj);
	97	}
	98
	99	//____________________________________________________________________
	100	AliFMDSharingFilter::~AliFMDSharingFilter()
	101	{
7984e5f7	102	//
	103	// Destructor
	104	//
7e4038b5	105	fRingHistos.Delete();
	106	}
	107
	108	//____________________________________________________________________
	109	AliFMDSharingFilter&
	110	AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
	111	{
7984e5f7	112	//
	113	// Assignment operator
	114	//
	115	// Parameters:
	116	// o Object to assign from
	117	//
	118	// Return:
	119	// Reference to this
	120	//
ea3e5d95	121	TNamed::operator=(o);
7e4038b5	122
	123	fLowCut = o.fLowCut;
	124	fCorrectAngles = o.fCorrectAngles;
0bd4b00f	125	fNXi = o.fNXi;
ea3e5d95	126	fDebug = o.fDebug;
7e4038b5	127
	128	fRingHistos.Delete();
	129	TIter next(&o.fRingHistos);
	130	TObject* obj = 0;
	131	while ((obj = next())) fRingHistos.Add(obj);
	132
	133	return *this;
	134	}
	135
	136	//____________________________________________________________________
	137	AliFMDSharingFilter::RingHistos*
	138	AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
	139	{
7984e5f7	140	//
	141	// Get the ring histogram container
	142	//
	143	// Parameters:
	144	// d Detector
	145	// r Ring
	146	//
	147	// Return:
	148	// Ring histogram container
	149	//
7e4038b5	150	Int_t idx = -1;
	151	switch (d) {
	152	case 1: idx = 0; break;
	153	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	154	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	155	}
	156	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
	157
	158	return static_cast<RingHistos*>(fRingHistos.At(idx));
	159	}
	160
	161	//____________________________________________________________________
	162	Bool_t
	163	AliFMDSharingFilter::Filter(const AliESDFMD& input,
	164	Bool_t lowFlux,
ea3e5d95	165	AliESDFMD& output)
7e4038b5	166	{
7984e5f7	167	//
	168	// Filter the input AliESDFMD object
	169	//
	170	// Parameters:
	171	// input Input
	172	// lowFlux If this is a low-flux event
	173	// output Output AliESDFMD object
	174	//
	175	// Return:
	176	// True on success, false otherwise
	177	//
7e4038b5	178	output.Clear();
7e4038b5	179	TIter next(&fRingHistos);
0bd4b00f	180	RingHistos* o = 0;
0bd4b00f	181	Double_t lowCut = GetLowCut();
7e4038b5	182	while ((o = static_cast<RingHistos*>(next())))
	183	o->Clear();
	184
ea3e5d95	185	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
7e4038b5	186
	187	for(UShort_t d = 1; d <= 3; d++) {
	188	Int_t nRings = (d == 1 ? 1 : 2);
	189	for (UShort_t q = 0; q < nRings; q++) {
	190	Char_t r = (q == 0 ? 'I' : 'O');
	191	UShort_t nsec = (q == 0 ? 20 : 40);
	192	UShort_t nstr = (q == 0 ? 512 : 256);
	193
	194	RingHistos* histos = GetRingHistos(d, r);
	195
	196	for(UShort_t s = 0; s < nsec; s++) {
	197	Bool_t usedThis = kFALSE;
	198	Bool_t usedPrev = kFALSE;
	199
	200	for(UShort_t t = 0; t < nstr; t++) {
	201	output.SetMultiplicity(d,r,s,t,0.);
	202	Float_t mult = SignalInStrip(input,d,r,s,t);
	203
	204	// Keep dead-channel information.
	205	if(mult == AliESDFMD::kInvalidMult)
	206	output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
	207
	208	// If no signal or dead strip, go on.
	209	if (mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
	210
	211	// Get the pseudo-rapidity
ea3e5d95	212	Double_t eta = input.Eta(d,r,s,t);
7e4038b5	213
	214	// Fill the diagnostics histogram
	215	histos->fBefore->Fill(mult);
	216
	217	// Get next and previous signal - if any
	218	Double_t prevE = 0;
	219	Double_t nextE = 0;
	220	if (t != 0) {
	221	prevE = SignalInStrip(input,d,r,s,t-1);
	222	if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
	223	}
	224	if (t != nstr - 1) {
	225	nextE = SignalInStrip(input,d,r,s,t+1);
	226	if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
	227	}
	228
	229	Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
	230	lowFlux,d,r,s,t,
	231	usedPrev,usedThis);
0bd4b00f	232	if (mergedEnergy > lowCut) histos->Incr();
7e4038b5	233	histos->fAfter->Fill(mergedEnergy);
	234
	235	output.SetMultiplicity(d,r,s,t,mergedEnergy);
	236	output.SetEta(d,r,s,t,eta);
	237	} // for strip
	238	} // for sector
	239	} // for ring
	240	} // for detector
	241
	242	next.Reset();
	243	while ((o = static_cast<RingHistos*>(next())))
	244	o->Finish();
	245
	246	return kTRUE;
	247	}
	248
	249	//_____________________________________________________________________
	250	Double_t
	251	AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
	252	UShort_t d,
	253	Char_t r,
	254	UShort_t s,
	255	UShort_t t) const
	256	{
7984e5f7	257	//
	258	// Get the signal in a strip
	259	//
	260	// Parameters:
	261	// fmd ESD object
	262	// d Detector
	263	// r Ring
	264	// s Sector
	265	// t Strip
	266	//
	267	// Return:
	268	// The energy signal
	269	//
7e4038b5	270	Double_t mult = input.Multiplicity(d,r,s,t);
	271	if (mult == AliESDFMD::kInvalidMult \|\| mult == 0) return mult;
	272
	273	if (fCorrectAngles && !input.IsAngleCorrected())
	274	mult = AngleCorrect(mult, input.Eta(d,r,s,t));
	275	else if (!fCorrectAngles && input.IsAngleCorrected())
	276	mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
	277	return mult;
	278	}
0bd4b00f	279	//_____________________________________________________________________
	280	Double_t
	281	AliFMDSharingFilter::GetLowCut() const
	282	{
7984e5f7	283	//
	284	// Get the low cut. Normally, the low cut is taken to be the lower
	285	// value of the fit range used when generating the energy loss fits.
	286	// However, if fLowCut is set (using SetLowCit) to a value greater
	287	// than 0, then that value is used.
	288	//
0bd4b00f	289	if (fLowCut > 0) return fLowCut;
	290	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	291	AliFMDCorrELossFit* fits = fcm.GetELossFit();
	292	return fits->GetLowCut();
	293	}
7e4038b5	294
	295	//_____________________________________________________________________
	296	Double_t
	297	AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
	298	{
7984e5f7	299	//
	300	// Get the high cut. The high cut is defined as the
	301	// most-probably-value peak found from the energy distributions, minus
	302	// 2 times the width of the corresponding Landau.
	303	//
0bd4b00f	304	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
0bd4b00f	305
7e4038b5	306
	307	// Get the high cut. The high cut is defined as the
	308	// most-probably-value peak found from the energy distributions, minus
	309	// 2 times the width of the corresponding Landau.
0bd4b00f	310	AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
	311	Double_t delta = 1024;
	312	Double_t xi = 1024;
	313	if (!fit) {
	314	AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
	315	}
	316	else {
	317	delta = fit->fDelta;
	318	xi = fit->fXi;
	319	}
	320
	321	if (delta > 100) {
	322	AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
	323	}
	324	Double_t highCut = (delta - fNXi * xi);
7e4038b5	325	return highCut;
	326	}
	327
	328	//_____________________________________________________________________
	329	Double_t
	330	AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
	331	Double_t eta,
	332	Double_t prevE,
	333	Double_t nextE,
	334	Bool_t lowFlux,
	335	UShort_t d,
	336	Char_t r,
	337	UShort_t /s/,
	338	UShort_t /t/,
	339	Bool_t& usedPrev,
fb3430ac	340	Bool_t& usedThis) const
7e4038b5	341	{
7984e5f7	342	//
	343	// The actual algorithm
	344	//
	345	// Parameters:
	346	// mult The unfiltered signal in the strip
	347	// eta Psuedo rapidity
	348	// prevE Previous strip signal (or 0)
	349	// nextE Next strip signal (or 0)
	350	// lowFlux Whether this is a low flux event
	351	// d Detector
	352	// r Ring
	353	// s Sector
	354	// t Strip
	355	// usedPrev Whether the previous strip was used in sharing or not
	356	// usedThis Wether this strip was used in sharing or not.
	357	//
	358	// Return:
	359	// The filtered signal in the strip
	360	//
	361
7e4038b5	362	// IF the multiplicity is very large, or below the cut, reject it, and
7e4038b5	363	// flag it as candidate for sharing
0bd4b00f	364	Double_t lowCut = GetLowCut();
0bd4b00f	365	if(mult > 12 \|\| mult < lowCut) {
7e4038b5	366	usedThis = kFALSE;
	367	usedPrev = kFALSE;
	368	return 0;
	369	}
	370
	371	// If this was shared with the previous signal, return 0 and mark it as
	372	// not a candiate for sharing
	373	if(usedThis) {
	374	usedThis = kFALSE;
	375	usedPrev = kTRUE;
	376	return 0.;
	377	}
	378
	379	//analyse and perform sharing on one strip
	380
0bd4b00f	381	// Get the high cut
7e4038b5	382	Double_t highCut = GetHighCut(d, r, eta);
	383
	384	// If this signal is smaller than the next, and the next signal is smaller
	385	// than then the high cut, and it's a low-flux event, then mark this and
	386	// the next signal as candidates for sharing
	387	//
	388	// This is the test if the signal is the smaller of two possibly
	389	// shared signals.
	390	if (mult < nextE &&
	391	nextE > highCut &&
	392	lowFlux) {
	393	usedThis = kFALSE;
	394	usedPrev = kFALSE;
	395	return 0;
	396	}
	397
	398	Double_t totalE = mult;
	399
	400	// If the previous signal was larger than the low cut, and
	401	// the previous was smaller than high cut, and the previous signal
	402	// isn't marked as used, then add it's energy to this signal
0bd4b00f	403	if (prevE > lowCut &&
7e4038b5	404	prevE < highCut &&
	405	!usedPrev)
	406	totalE += prevE;
	407
	408	// If the next signal is larger than the low cut, and
	409	// the next signal is smaller than the low cut, then add the next signal
	410	// to this, and marked the next signal as used
0bd4b00f	411	if(nextE > lowCut && nextE < highCut ) {
7e4038b5	412	totalE += nextE;
	413	usedThis = kTRUE;
	414	}
	415
	416	// If we're processing on non-angle corrected data, we should do the
	417	// angle correction here
	418	if (!fCorrectAngles)
	419	totalE = AngleCorrect(totalE, eta);
	420
	421	// Fill post processing histogram
	422	// if(totalE > fLowCut)
	423	// GetRingHistos(d,r)->fAfter->Fill(totalE);
	424
	425	// Return value
	426	Double_t mergedEnergy = 0;
	427
	428	if(totalE > 0) {
	429	// If we have a signal, then this is used (sharedPrev=true) and
	430	// the signal is set to the result
	431	mergedEnergy = totalE;
	432	usedPrev = kTRUE;
	433	}
	434	else {
	435	// If we do not have a signal (too low), then this is not shared,
	436	// and the next strip is not shared either
	437	usedThis = kFALSE;
	438	usedPrev = kFALSE;
	439	}
	440
	441	return mergedEnergy;
	442	}
	443	//____________________________________________________________________
	444	Double_t
	445	AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
	446	{
7984e5f7	447	//
	448	// Angle correct the signal
	449	//
	450	// Parameters:
	451	// mult Angle Un-corrected Signal
	452	// eta Pseudo-rapidity
	453	//
	454	// Return:
	455	// Angle corrected signal
	456	//
7e4038b5	457	Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
	458	if (eta < 0) theta -= TMath::Pi();
	459	return mult * TMath::Cos(theta);
	460	}
	461	//____________________________________________________________________
	462	Double_t
	463	AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
	464	{
7984e5f7	465	//
	466	// Angle de-correct the signal
	467	//
	468	// Parameters:
	469	// mult Angle corrected Signal
	470	// eta Pseudo-rapidity
	471	//
	472	// Return:
	473	// Angle un-corrected signal
	474	//
7e4038b5	475	Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
	476	if (eta < 0) theta -= TMath::Pi();
	477	return mult / TMath::Cos(theta);
	478	}
	479
	480	//____________________________________________________________________
	481	void
fb3430ac	482	AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
7e4038b5	483	{
7984e5f7	484	//
	485	// Scale the histograms to the total number of events
	486	//
	487	// Parameters:
	488	// dir Where the output is
	489	// nEvents Number of events
	490	//
7e4038b5	491	if (nEvents <= 0) return;
9d99b0dd	492	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
9d99b0dd	493	if (!d) return;
7e4038b5	494
	495	TIter next(&fRingHistos);
	496	RingHistos* o = 0;
9d99b0dd	497	while ((o = static_cast<RingHistos*>(next())))
9d99b0dd	498	o->ScaleHistograms(d, nEvents);
7e4038b5	499	}
	500
	501	//____________________________________________________________________
	502	void
9d99b0dd	503	AliFMDSharingFilter::DefineOutput(TList* dir)
7e4038b5	504	{
7984e5f7	505	//
	506	// Define the output histograms. These are put in a sub list of the
	507	// passed list. The histograms are merged before the parent task calls
	508	// AliAnalysisTaskSE::Terminate
	509	//
	510	// Parameters:
	511	// dir Directory to add to
	512	//
7e4038b5	513	TList* d = new TList;
	514	d->SetName(GetName());
	515	dir->Add(d);
7e4038b5	516	TIter next(&fRingHistos);
	517	RingHistos* o = 0;
	518	while ((o = static_cast<RingHistos*>(next()))) {
	519	o->Output(d);
	520	}
	521	}
0bd4b00f	522	//____________________________________________________________________
	523	void
	524	AliFMDSharingFilter::Print(Option_t* /option/) const
	525	{
7984e5f7	526	//
	527	// Print information
	528	//
	529	// Parameters:
	530	// option Not used
	531	//
0bd4b00f	532	char ind[gROOT->GetDirLevel()+1];
	533	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	534	ind[gROOT->GetDirLevel()] = '\0';
	535	std::cout << ind << "AliFMDSharingFilter: " << GetName() << '\n'
	536	<< ind << " Low cut: " << fLowCut << '\n'
	537	<< ind << " N xi factor: " << fNXi << '\n'
	538	<< ind << " Use corrected angles: "
	539	<< (fCorrectAngles ? "yes" : "no") << std::endl;
	540	}
7e4038b5	541
	542	//====================================================================
	543	AliFMDSharingFilter::RingHistos::RingHistos()
9d99b0dd	544	: AliForwardUtil::RingHistos(),
7e4038b5	545	fBefore(0),
	546	fAfter(0),
	547	fHits(0),
	548	fNHits(0)
7984e5f7	549	{
	550	//
	551	// Default CTOR
	552	//
	553
	554	}
7e4038b5	555
	556	//____________________________________________________________________
	557	AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd	558	: AliForwardUtil::RingHistos(d,r),
7e4038b5	559	fBefore(0),
	560	fAfter(0),
	561	fHits(0),
	562	fNHits(0)
	563	{
7984e5f7	564	//
	565	// Constructor
	566	//
	567	// Parameters:
	568	// d detector
	569	// r ring
	570	//
9d99b0dd	571	fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
9d99b0dd	572	Form("Energy loss in %s (reconstruction)", fName.Data()),
7e4038b5	573	1000, 0, 25);
9d99b0dd	574	fAfter = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
	575	Form("Energy loss in %s (sharing corrected)",
	576	fName.Data()), 1000, 0, 25);
7e4038b5	577	fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
	578	fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
	579	fBefore->SetFillColor(kRed+1);
	580	fBefore->SetFillStyle(3001);
	581	// fBefore->Sumw2();
	582	fBefore->SetDirectory(0);
	583	fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
	584	fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
	585	fAfter->SetFillColor(kBlue+1);
	586	fAfter->SetFillStyle(3001);
	587	// fAfter->Sumw2();
	588	fAfter->SetDirectory(0);
	589
9d99b0dd	590	fHits = new TH1D(Form("%s_Hits", fName.Data()),
9d99b0dd	591	Form("Number of hits in %s", fName.Data()),
7e4038b5	592	200, 0, 200000);
	593	fHits->SetDirectory(0);
	594	fHits->SetXTitle("# of hits");
	595	fHits->SetFillColor(kGreen+1);
	596	}
	597	//____________________________________________________________________
	598	AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd	599	: AliForwardUtil::RingHistos(o),
7e4038b5	600	fBefore(o.fBefore),
	601	fAfter(o.fAfter),
	602	fHits(o.fHits),
	603	fNHits(o.fNHits)
7984e5f7	604	{
	605	//
	606	// Copy constructor
	607	//
	608	// Parameters:
	609	// o Object to copy from
	610	//
	611	}
7e4038b5	612
	613	//____________________________________________________________________
	614	AliFMDSharingFilter::RingHistos&
	615	AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
	616	{
7984e5f7	617	//
	618	// Assignment operator
	619	//
	620	// Parameters:
	621	// o Object to assign from
	622	//
	623	// Return:
	624	// Reference to this
	625	//
9d99b0dd	626	AliForwardUtil::RingHistos::operator=(o);
7e4038b5	627	fDet = o.fDet;
	628	fRing = o.fRing;
	629
	630	if (fBefore) delete fBefore;
	631	if (fAfter) delete fAfter;
	632	if (fHits) delete fHits;
	633
	634	fBefore = static_cast<TH1D*>(o.fBefore->Clone());
	635	fAfter = static_cast<TH1D*>(o.fAfter->Clone());
	636	fHits = static_cast<TH1D*>(o.fHits->Clone());
	637
	638	return *this;
	639	}
	640	//____________________________________________________________________
	641	AliFMDSharingFilter::RingHistos::~RingHistos()
	642	{
7984e5f7	643	//
	644	// Destructor
	645	//
7e4038b5	646	if (fBefore) delete fBefore;
	647	if (fAfter) delete fAfter;
	648	if (fHits) delete fHits;
	649	}
	650	//____________________________________________________________________
	651	void
	652	AliFMDSharingFilter::RingHistos::Finish()
	653	{
7984e5f7	654	//
	655	// Finish off
	656	//
	657	//
7e4038b5	658	fHits->Fill(fNHits);
	659	}
	660
9d99b0dd	661	//____________________________________________________________________
9d99b0dd	662	void
fb3430ac	663	AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir,
fb3430ac	664	Int_t nEvents)
9d99b0dd	665	{
7984e5f7	666	//
	667	// Scale the histograms to the total number of events
	668	//
	669	// Parameters:
	670	// nEvents Number of events
	671	// dir Where the output is
	672	//
9d99b0dd	673	TList* l = GetOutputList(dir);
	674	if (!l) return;
	675
	676	TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
	677	fName.Data())));
	678	TH1D* after = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
	679	fName.Data())));
	680	if (before) before->Scale(1./nEvents);
	681	if (after) after->Scale(1./nEvents);
	682	}
	683
7e4038b5	684	//____________________________________________________________________
	685	void
	686	AliFMDSharingFilter::RingHistos::Output(TList* dir)
	687	{
7984e5f7	688	//
	689	// Make output
	690	//
	691	// Parameters:
	692	// dir where to store
	693	//
9d99b0dd	694	TList* d = DefineOutputList(dir);
7e4038b5	695	d->Add(fBefore);
	696	d->Add(fAfter);
	697	d->Add(fHits);
	698	dir->Add(d);
	699	}
	700
	701	//____________________________________________________________________
	702	//
	703	// EOF
	704	//