[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 "AliForwardCorrectionManager.h"
// #include "AliFMDAnaParameters.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)
{}

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

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

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

//____________________________________________________________________
AliFMDSharingFilter&
AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
{
  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
{
  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)
{
  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
{
  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
{
  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
{
  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) 
{
  // 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
{
  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(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->ScaleHistograms(d, nEvents);
}

//____________________________________________________________________
void
AliFMDSharingFilter::DefineOutput(TList* dir)
{
  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
{
  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)
{}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r), 
    fBefore(0), 
    fAfter(0),
    fHits(0),
    fNHits(0)
{
  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)
{}

//____________________________________________________________________
AliFMDSharingFilter::RingHistos&
AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
{
  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()
{
  if (fBefore) delete fBefore;
  if (fAfter)  delete fAfter;
  if (fHits)   delete fHits;
}
//____________________________________________________________________
void
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 = DefineOutputList(dir);
  d->Add(fBefore);
  d->Add(fAfter);
  d->Add(fHits);
  dir->Add(d);
}

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