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

//
// Class to do the sharing correction of FMD ESD data
//
#include "AliFMDSharingFilter.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TH1.h>
#include <TMath.h>
#include "AliFMDAnaParameters.h"
#include <AliLog.h>

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


//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter()
  : TNamed(), 
    fRingHistos(),
    fLowCut(0.3),
    fCorrectAngles(kFALSE),
    fEtaCorr(0)
{}

//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
  : TNamed("fmdSharingFilter", title), 
    fRingHistos(), 
    fLowCut(0.3),
    fCorrectAngles(kFALSE),
    fEtaCorr(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");
}

//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
  : TNamed(o), 
    fRingHistos(), 
    fLowCut(o.fLowCut),
    fCorrectAngles(o.fCorrectAngles),
    fEtaCorr(o.fEtaCorr)
{
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

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

//____________________________________________________________________
AliFMDSharingFilter&
AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
{
  SetName(o.GetName());
  SetTitle(o.GetTitle());

  fLowCut        = o.fLowCut;
  fCorrectAngles = o.fCorrectAngles;

  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
{
  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,
			    Double_t         vz)
{
  output.Clear();
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  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 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;
	  
	  // 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 > fLowCut) 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
{
  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::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::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);
  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) 
{
  // IF the multiplicity is very large, or below the cut, reject it, and 
  // flag it as candidate for sharing 
  if(mult > 12 || mult < fLowCut) {
    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
  
  // Calculate the total energy loss 
  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 > fLowCut && 
      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 > fLowCut && 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
{
  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
{
  Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
  if (eta < 0) theta -= TMath::Pi();
  return mult / TMath::Cos(theta);
}

//____________________________________________________________________
void
AliFMDSharingFilter::ScaleHistograms(Int_t nEvents)
{
  if (nEvents <= 0) return;

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->fBefore->Scale(1. / nEvents);
    o->fAfter->Scale(1. / nEvents);
  }
}

//____________________________________________________________________
void
AliFMDSharingFilter::Output(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);
  }
}
  
//====================================================================
AliFMDSharingFilter::RingHistos::RingHistos()
  : fDet(0),
    fRing('\0'),
    fBefore(0), 
    fAfter(0), 
    fHits(0),
    fNHits(0)
{}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
  : fDet(d), 
    fRing(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), 
		     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("FMD%d%c_Hits", d, r), 
		   Form("Number of hits in FMD%d%c", d, r), 
		   200, 0, 200000);
  fHits->SetDirectory(0);
  fHits->SetXTitle("# of hits");
  fHits->SetFillColor(kGreen+1);
}
//____________________________________________________________________
AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
  : TObject(o), 
    fDet(o.fDet), 
    fRing(o.fRing), 
    fBefore(o.fBefore), 
    fAfter(o.fAfter),
    fHits(o.fHits),
    fNHits(o.fNHits)
{}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos&
AliFMDSharingFilter::RingHistos::operator=(const RingHistos& 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()
{
  if (fBefore) delete fBefore;
  if (fAfter)  delete fAfter;
  if (fHits)   delete fHits;
}
//____________________________________________________________________
void
AliFMDSharingFilter::RingHistos::Finish()
{
  fHits->Fill(fNHits);
}

//____________________________________________________________________
void
AliFMDSharingFilter::RingHistos::Output(TList* dir)
{
  TList* d = new TList;
  d->SetName(Form("FMD%d%c", fDet, fRing)); 
  d->Add(fBefore);
  d->Add(fAfter);
  d->Add(fHits);
  dir->Add(d);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
	1	//
	2	// Class to do the sharing correction of FMD ESD data
	3	//
	4	#include "AliFMDSharingFilter.h"
	5	#include <AliESDFMD.h>
	6	#include <TAxis.h>
	7	#include <TList.h>
	8	#include <TH1.h>
	9	#include <TMath.h>
	10	#include "AliFMDAnaParameters.h"
	11	#include <AliLog.h>
	12
	13	ClassImp(AliFMDSharingFilter)
	14	#if 0
	15	; // This is for Emacs
	16	#endif
	17
	18
	19	//____________________________________________________________________
	20	AliFMDSharingFilter::AliFMDSharingFilter()
	21	: TNamed(),
	22	fRingHistos(),
	23	fLowCut(0.3),
	24	fCorrectAngles(kFALSE),
	25	fEtaCorr(0)
	26	{}
	27
	28	//____________________________________________________________________
	29	AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
	30	: TNamed("fmdSharingFilter", title),
	31	fRingHistos(),
	32	fLowCut(0.3),
	33	fCorrectAngles(kFALSE),
	34	fEtaCorr(0)
	35	{
	36	fRingHistos.Add(new RingHistos(1, 'I'));
	37	fRingHistos.Add(new RingHistos(2, 'I'));
	38	fRingHistos.Add(new RingHistos(2, 'O'));
	39	fRingHistos.Add(new RingHistos(3, 'I'));
	40	fRingHistos.Add(new RingHistos(3, 'O'));
	41	fEtaCorr = new TH2F("etaCorr", "Correction of #eta",
	42	200, -4, 6, 200, -4, 6);
	43	fEtaCorr->SetXTitle("#eta from ESD");
	44	fEtaCorr->SetYTitle("#eta from AnaParameters");
	45	}
	46
	47	//____________________________________________________________________
	48	AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
	49	: TNamed(o),
	50	fRingHistos(),
	51	fLowCut(o.fLowCut),
	52	fCorrectAngles(o.fCorrectAngles),
	53	fEtaCorr(o.fEtaCorr)
	54	{
	55	TIter next(&o.fRingHistos);
	56	TObject* obj = 0;
	57	while ((obj = next())) fRingHistos.Add(obj);
	58	}
	59
	60	//____________________________________________________________________
	61	AliFMDSharingFilter::~AliFMDSharingFilter()
	62	{
	63	fRingHistos.Delete();
	64	}
	65
	66	//____________________________________________________________________
	67	AliFMDSharingFilter&
	68	AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
	69	{
	70	SetName(o.GetName());
	71	SetTitle(o.GetTitle());
	72
	73	fLowCut = o.fLowCut;
	74	fCorrectAngles = o.fCorrectAngles;
	75
	76	fRingHistos.Delete();
	77	TIter next(&o.fRingHistos);
	78	TObject* obj = 0;
	79	while ((obj = next())) fRingHistos.Add(obj);
	80
	81	return *this;
	82	}
	83
	84	//____________________________________________________________________
	85	AliFMDSharingFilter::RingHistos*
	86	AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
	87	{
	88	Int_t idx = -1;
	89	switch (d) {
	90	case 1: idx = 0; break;
	91	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	92	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	93	}
	94	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
	95
	96	return static_cast<RingHistos*>(fRingHistos.At(idx));
	97	}
	98
	99	//____________________________________________________________________
	100	Bool_t
	101	AliFMDSharingFilter::Filter(const AliESDFMD& input,
	102	Bool_t lowFlux,
	103	AliESDFMD& output,
	104	Double_t vz)
	105	{
	106	output.Clear();
	107	TIter next(&fRingHistos);
	108	RingHistos* o = 0;
	109	while ((o = static_cast<RingHistos*>(next())))
	110	o->Clear();
	111
	112	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	113
	114	for(UShort_t d = 1; d <= 3; d++) {
	115	Int_t nRings = (d == 1 ? 1 : 2);
	116	for (UShort_t q = 0; q < nRings; q++) {
	117	Char_t r = (q == 0 ? 'I' : 'O');
	118	UShort_t nsec = (q == 0 ? 20 : 40);
	119	UShort_t nstr = (q == 0 ? 512 : 256);
	120
	121	RingHistos* histos = GetRingHistos(d, r);
	122
	123	for(UShort_t s = 0; s < nsec; s++) {
	124	Bool_t usedThis = kFALSE;
	125	Bool_t usedPrev = kFALSE;
	126
	127	for(UShort_t t = 0; t < nstr; t++) {
	128	output.SetMultiplicity(d,r,s,t,0.);
	129	Float_t mult = SignalInStrip(input,d,r,s,t);
	130
	131	// Keep dead-channel information.
	132	if(mult == AliESDFMD::kInvalidMult)
	133	output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
	134
	135	// If no signal or dead strip, go on.
	136	if (mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
	137
	138	// Get the pseudo-rapidity
	139	Double_t eta1 = input.Eta(d,r,s,t);
	140	Double_t eta2 = pars->GetEtaFromStrip(d,r,s,t,vz);
	141	fEtaCorr->Fill(eta1, eta2);
	142	Double_t eta = eta2;
	143
	144	// Fill the diagnostics histogram
	145	histos->fBefore->Fill(mult);
	146
	147	// Get next and previous signal - if any
	148	Double_t prevE = 0;
	149	Double_t nextE = 0;
	150	if (t != 0) {
	151	prevE = SignalInStrip(input,d,r,s,t-1);
	152	if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
	153	}
	154	if (t != nstr - 1) {
	155	nextE = SignalInStrip(input,d,r,s,t+1);
	156	if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
	157	}
	158
	159	Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
	160	lowFlux,d,r,s,t,
	161	usedPrev,usedThis);
	162	if (mergedEnergy > fLowCut) histos->Incr();
	163	histos->fAfter->Fill(mergedEnergy);
	164
	165	output.SetMultiplicity(d,r,s,t,mergedEnergy);
	166	output.SetEta(d,r,s,t,eta);
	167	} // for strip
	168	} // for sector
	169	} // for ring
	170	} // for detector
	171
	172	next.Reset();
	173	while ((o = static_cast<RingHistos*>(next())))
	174	o->Finish();
	175
	176	return kTRUE;
	177	}
	178
	179	//_____________________________________________________________________
	180	Double_t
	181	AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
	182	UShort_t d,
	183	Char_t r,
	184	UShort_t s,
	185	UShort_t t) const
	186	{
	187	Double_t mult = input.Multiplicity(d,r,s,t);
	188	if (mult == AliESDFMD::kInvalidMult \|\| mult == 0) return mult;
	189
	190	if (fCorrectAngles && !input.IsAngleCorrected())
	191	mult = AngleCorrect(mult, input.Eta(d,r,s,t));
	192	else if (!fCorrectAngles && input.IsAngleCorrected())
	193	mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
	194	return mult;
	195	}
	196
	197	//_____________________________________________________________________
	198	Double_t
	199	AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
	200	{
	201	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	202
	203	// Get the high cut. The high cut is defined as the
	204	// most-probably-value peak found from the energy distributions, minus
	205	// 2 times the width of the corresponding Landau.
	206	Double_t mpv = pars->GetMPV(d,r,eta);
	207	Double_t w = pars->GetSigma(d,r,eta);
	208	if (mpv > 100)
	209	AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
	210	Double_t highCut = (mpv - 2 * w);
	211	return highCut;
	212	}
	213
	214	//_____________________________________________________________________
	215	Double_t
	216	AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
	217	Double_t eta,
	218	Double_t prevE,
	219	Double_t nextE,
	220	Bool_t lowFlux,
	221	UShort_t d,
	222	Char_t r,
	223	UShort_t /s/,
	224	UShort_t /t/,
	225	Bool_t& usedPrev,
	226	Bool_t& usedThis)
	227	{
	228	// IF the multiplicity is very large, or below the cut, reject it, and
	229	// flag it as candidate for sharing
	230	if(mult > 12 \|\| mult < fLowCut) {
	231	usedThis = kFALSE;
	232	usedPrev = kFALSE;
	233	return 0;
	234	}
	235
	236	// If this was shared with the previous signal, return 0 and mark it as
	237	// not a candiate for sharing
	238	if(usedThis) {
	239	usedThis = kFALSE;
	240	usedPrev = kTRUE;
	241	return 0.;
	242	}
	243
	244	//analyse and perform sharing on one strip
	245
	246	// Calculate the total energy loss
	247	Double_t highCut = GetHighCut(d, r, eta);
	248
	249	// If this signal is smaller than the next, and the next signal is smaller
	250	// than then the high cut, and it's a low-flux event, then mark this and
	251	// the next signal as candidates for sharing
	252	//
	253	// This is the test if the signal is the smaller of two possibly
	254	// shared signals.
	255	if (mult < nextE &&
	256	nextE > highCut &&
	257	lowFlux) {
	258	usedThis = kFALSE;
	259	usedPrev = kFALSE;
	260	return 0;
	261	}
	262
	263	Double_t totalE = mult;
	264
	265	// If the previous signal was larger than the low cut, and
	266	// the previous was smaller than high cut, and the previous signal
	267	// isn't marked as used, then add it's energy to this signal
	268	if (prevE > fLowCut &&
	269	prevE < highCut &&
	270	!usedPrev)
	271	totalE += prevE;
	272
	273	// If the next signal is larger than the low cut, and
	274	// the next signal is smaller than the low cut, then add the next signal
	275	// to this, and marked the next signal as used
	276	if(nextE > fLowCut && nextE < highCut ) {
	277	totalE += nextE;
	278	usedThis = kTRUE;
	279	}
	280
	281	// If we're processing on non-angle corrected data, we should do the
	282	// angle correction here
	283	if (!fCorrectAngles)
	284	totalE = AngleCorrect(totalE, eta);
	285
	286	// Fill post processing histogram
	287	// if(totalE > fLowCut)
	288	// GetRingHistos(d,r)->fAfter->Fill(totalE);
	289
	290	// Return value
	291	Double_t mergedEnergy = 0;
	292
	293	if(totalE > 0) {
	294	// If we have a signal, then this is used (sharedPrev=true) and
	295	// the signal is set to the result
	296	mergedEnergy = totalE;
	297	usedPrev = kTRUE;
	298	}
	299	else {
	300	// If we do not have a signal (too low), then this is not shared,
	301	// and the next strip is not shared either
	302	usedThis = kFALSE;
	303	usedPrev = kFALSE;
	304	}
	305
	306	return mergedEnergy;
	307	}
	308	//____________________________________________________________________
	309	Double_t
	310	AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
	311	{
	312	Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
	313	if (eta < 0) theta -= TMath::Pi();
	314	return mult * TMath::Cos(theta);
	315	}
	316	//____________________________________________________________________
	317	Double_t
	318	AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
	319	{
	320	Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
	321	if (eta < 0) theta -= TMath::Pi();
	322	return mult / TMath::Cos(theta);
	323	}
	324
	325	//____________________________________________________________________
	326	void
	327	AliFMDSharingFilter::ScaleHistograms(Int_t nEvents)
	328	{
	329	if (nEvents <= 0) return;
	330
	331	TIter next(&fRingHistos);
	332	RingHistos* o = 0;
	333	while ((o = static_cast<RingHistos*>(next()))) {
	334	o->fBefore->Scale(1. / nEvents);
	335	o->fAfter->Scale(1. / nEvents);
	336	}
	337	}
	338
	339	//____________________________________________________________________
	340	void
	341	AliFMDSharingFilter::Output(TList* dir)
	342	{
	343	TList* d = new TList;
	344	d->SetName(GetName());
	345	dir->Add(d);
	346	d->Add(fEtaCorr);
	347	TIter next(&fRingHistos);
	348	RingHistos* o = 0;
	349	while ((o = static_cast<RingHistos*>(next()))) {
	350	o->Output(d);
	351	}
	352	}
	353
	354	//====================================================================
	355	AliFMDSharingFilter::RingHistos::RingHistos()
	356	: fDet(0),
	357	fRing('\0'),
	358	fBefore(0),
	359	fAfter(0),
	360	fHits(0),
	361	fNHits(0)
	362	{}
	363
	364	//____________________________________________________________________
	365	AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
	366	: fDet(d),
	367	fRing(r),
	368	fBefore(0),
	369	fAfter(0),
	370	fHits(0),
	371	fNHits(0)
	372	{
	373	fBefore = new TH1D(Form("FMD%d%c_ESD_Eloss", d, r),
	374	Form("Energy loss in FMD%d%c (reconstruction)", d, r),
	375	1000, 0, 25);
	376	fAfter = new TH1D(Form("FMD%d%c_Ana_Eloss", d, r),
	377	Form("Energy loss in FMD%d%c (sharing corrected)", d, r),
	378	1000, 0, 25);
	379	fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
	380	fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
	381	fBefore->SetFillColor(kRed+1);
	382	fBefore->SetFillStyle(3001);
	383	// fBefore->Sumw2();
	384	fBefore->SetDirectory(0);
	385	fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
	386	fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
	387	fAfter->SetFillColor(kBlue+1);
	388	fAfter->SetFillStyle(3001);
	389	// fAfter->Sumw2();
	390	fAfter->SetDirectory(0);
	391
	392	fHits = new TH1D(Form("FMD%d%c_Hits", d, r),
	393	Form("Number of hits in FMD%d%c", d, r),
	394	200, 0, 200000);
	395	fHits->SetDirectory(0);
	396	fHits->SetXTitle("# of hits");
	397	fHits->SetFillColor(kGreen+1);
	398	}
	399	//____________________________________________________________________
	400	AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
	401	: TObject(o),
	402	fDet(o.fDet),
	403	fRing(o.fRing),
	404	fBefore(o.fBefore),
	405	fAfter(o.fAfter),
	406	fHits(o.fHits),
	407	fNHits(o.fNHits)
	408	{}
	409
	410	//____________________________________________________________________
	411	AliFMDSharingFilter::RingHistos&
	412	AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
	413	{
	414	fDet = o.fDet;
	415	fRing = o.fRing;
	416
	417	if (fBefore) delete fBefore;
	418	if (fAfter) delete fAfter;
	419	if (fHits) delete fHits;
	420
	421	fBefore = static_cast<TH1D*>(o.fBefore->Clone());
	422	fAfter = static_cast<TH1D*>(o.fAfter->Clone());
	423	fHits = static_cast<TH1D*>(o.fHits->Clone());
	424
	425	return *this;
	426	}
	427	//____________________________________________________________________
	428	AliFMDSharingFilter::RingHistos::~RingHistos()
	429	{
	430	if (fBefore) delete fBefore;
	431	if (fAfter) delete fAfter;
	432	if (fHits) delete fHits;
	433	}
	434	//____________________________________________________________________
	435	void
	436	AliFMDSharingFilter::RingHistos::Finish()
	437	{
	438	fHits->Fill(fNHits);
	439	}
	440
	441	//____________________________________________________________________
	442	void
	443	AliFMDSharingFilter::RingHistos::Output(TList* dir)
	444	{
	445	TList* d = new TList;
	446	d->SetName(Form("FMD%d%c", fDet, fRing));
	447	d->Add(fBefore);
	448	d->Add(fAfter);
	449	d->Add(fHits);
	450	dir->Add(d);
	451	}
	452
	453	//____________________________________________________________________
	454	//
	455	// EOF
	456	//