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

#include "AliFMDDensityCalculator.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
// #include "AliFMDAnaParameters.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TProfile.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

ClassImp(AliFMDDensityCalculator)
#if 0
; // For Emacs
#endif 

//____________________________________________________________________
AliFMDDensityCalculator::AliFMDDensityCalculator()
  : TNamed(), 
    fRingHistos(),
    fMultCut(0),
    fSumOfWeights(0),
    fWeightedSum(0),
    fCorrections(0),
    fMaxParticles(5),
    fAccI(0),
    fAccO(0),
    fDebug(0)
{}

//____________________________________________________________________
AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
  : TNamed("fmdDensityCalculator", title), 
    fRingHistos(), 
    fMultCut(0),
    fSumOfWeights(0),
    fWeightedSum(0),
    fCorrections(0),
    fMaxParticles(5),
    fAccI(0),
    fAccO(0),
    fDebug(0)
{
  fRingHistos.SetName(GetName());
  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'));
  fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
			   200, 0, 20);
  fSumOfWeights->SetFillColor(kRed+1);
  fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
  fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
			   200, 0, 20);
  fWeightedSum->SetFillColor(kBlue+1);
  fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
  fCorrections  = new TH1D("corrections", "Distribution of corrections", 
			   100, 0, 10);
  fCorrections->SetFillColor(kBlue+1);
  fCorrections->SetXTitle("correction");

  fAccI = GenerateAcceptanceCorrection('I');
  fAccO = GenerateAcceptanceCorrection('O');
}

//____________________________________________________________________
AliFMDDensityCalculator::AliFMDDensityCalculator(const 
						 AliFMDDensityCalculator& o)
  : TNamed(o), 
    fRingHistos(), 
    fMultCut(o.fMultCut),
    fSumOfWeights(o.fSumOfWeights),
    fWeightedSum(o.fWeightedSum),
    fCorrections(o.fCorrections),
    fMaxParticles(o.fMaxParticles),
    fAccI(o.fAccI),
    fAccO(o.fAccO),
    fDebug(o.fDebug)
{
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

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

//____________________________________________________________________
AliFMDDensityCalculator&
AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
{
  TNamed::operator=(o);

  fMultCut      = o.fMultCut;
  fDebug        = o.fDebug;
  fMaxParticles = o.fMaxParticles;
  fAccI         = o.fAccI;
  fAccO         = o.fAccO;

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

//____________________________________________________________________
AliFMDDensityCalculator::RingHistos*
AliFMDDensityCalculator::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));
}
    
//____________________________________________________________________
Double_t
AliFMDDensityCalculator::GetMultCut() const
{
  if (fMultCut > 0) return fMultCut;

  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
  AliFMDCorrELossFit* fits = fcm.GetELossFit();
  return fits->GetLowCut();
}
  
//____________________________________________________________________
Bool_t
AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
				   AliForwardUtil::Histos& hists,
				   UShort_t                vtxbin, 
				   Bool_t                  lowFlux)
{
  for (UShort_t d=1; d<=3; d++) { 
    UShort_t nr = (d == 1 ? 1 : 2);
    for (UShort_t q=0; q<nr; q++) { 
      Char_t      r = (q == 0 ? 'I' : 'O');
      UShort_t    ns= (q == 0 ?  20 :  40);
      UShort_t    nt= (q == 0 ? 512 : 256);
      TH2D*       h = hists.Get(d,r);
      RingHistos* rh= GetRingHistos(d,r);

      for (UShort_t s=0; s<ns; s++) { 
	for (UShort_t t=0; t<nt; t++) {
	  Float_t mult = fmd.Multiplicity(d,r,s,t);
	  
	  if (mult == 0 || mult > 20) continue;

	  Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
	  Float_t eta = fmd.Eta(d,r,s,t);
	  
	  Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
	  rh->fEvsN->Fill(mult,n);
	  rh->fEtaVsN->Fill(eta, n);

	  Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
	  fCorrections->Fill(c);
	  if (c > 0) n /= c;
	  rh->fEvsM->Fill(mult,n);
	  rh->fEtaVsM->Fill(eta, n);
	  rh->fCorr->Fill(eta, c);

	  h->Fill(eta,phi,n);
	  rh->fDensity->Fill(eta,phi,n);
	} // for t
      } // for s 
    } // for q
  } // for d
  
  return kTRUE;
}

//_____________________________________________________________________
Float_t 
AliFMDDensityCalculator::NParticles(Float_t  mult, 
				    UShort_t d, 
				    Char_t   r, 
				    UShort_t /*s*/, 
				    UShort_t /*t*/, 
				    UShort_t /*v*/, 
				    Float_t  eta,
				    Bool_t   lowFlux) const
{
  if (mult <= GetMultCut()) return 0;
  if (lowFlux) return 1;
  
  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
  AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
  if (!fit) { 
    AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
    return 0;
  }
  
  Int_t    m   = fit->FindMaxWeight();
  if (m < 1) { 
    AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
    return 0;
  }
  UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
  Double_t ret = fit->EvaluateWeighted(mult, n);

  if (fDebug > 10) {
    AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
  }

  fWeightedSum->Fill(ret);
  fSumOfWeights->Fill(ret);

  return ret;
#if 0
  Float_t mpv  = pars->GetMPV(d,r,eta);
  Float_t w    = pars->GetSigma(d,r,eta);
  Float_t w2   = pars->Get2MIPWeight(d,r,eta);
  Float_t w3   = pars->Get3MIPWeight(d,r,eta);
  Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
  Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
  
  Float_t sum  = (TMath::Landau(mult,mpv,w,kTRUE) +
		  w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
		  w3  * TMath::Landau(mult,mpv3,3*w,kTRUE));
  Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
		  2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
		  3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
  
  fWeightedSum->Fill(wsum);
  fSumOfWeights->Fill(sum);
  return (sum > 0) ? wsum / sum : 1;
#endif
}

//_____________________________________________________________________
Float_t 
AliFMDDensityCalculator::Correction(UShort_t d, 
				    Char_t   r, 
				    UShort_t /*s*/, 
				    UShort_t t, 
				    UShort_t /*v*/, 
				    Float_t  eta,
				    Bool_t   lowFlux) const
{
  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();

  Float_t correction = AcceptanceCorrection(r,t);
  if (lowFlux) { 
    // TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
    TH1D* dblHitCor = 0;
    if (fcm.GetDoubleHit()) 
      dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);

    if (dblHitCor) {
      Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
      if (dblC > 0) correction *= dblC;
    }
    else {
      AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
    }
  }
  
#if 0
  TH1F* deadCor = pars->GetFMDDeadCorrection(v);
  if (deadCor) { 
    Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
    if (deadC > 0) 
      correction *= deadC; 
  }
#endif  

  return correction;
}

//_____________________________________________________________________
TH1D*
AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
{
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };
  const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
  const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
  Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
  Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
  Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
  Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
  Float_t   basearc    = 2 * TMath::Pi() / nSec;
  Double_t  rad        = maxR - minR;
  Float_t   segment    = rad / nStrips;
  Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);

  // Numbers used to find end-point of strip.
  // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
  Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
  Float_t dx           = c2[0] - c1[0];
  Float_t dy           = c2[1] - c1[1];
  Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);

  TH1D* ret = new TH1D(Form("acc%c", r), 
		       Form("Acceptance correction for FMDx%c", r), 
		       nStrips, -.5, nStrips-.5);
  ret->SetXTitle("Strip");
  ret->SetYTitle("#varphi acceptance");
  ret->SetDirectory(0);
  ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
  ret->SetFillStyle(3001);

  for (Int_t t = 0; t < nStrips; t++) { 
    Float_t   radius     = minR + t * segment;
    
    // If the radius of the strip is smaller than the radius corresponding 
    // to the first corner we have a full strip length 
    if (radius <= cr) {
      ret->SetBinContent(t+1, 1);
      continue;
    }

    // Next, we should find the end-point of the strip - that is, 
    // the coordinates where the line from c1 to c2 intersects a circle 
    // with radius given by the strip. 
    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
    // Calculate the determinant 
    Float_t det = radius * radius * dr * dr - D*D;

    if (det <=  0) { 
      // <0 means No intersection
      // =0 means Exactly tangent
      ret->SetBinContent(t+1, 1);
      continue;
    }

    // Calculate end-point and the corresponding opening angle 
    Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
    Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
    Float_t th  = TMath::ATan2(x, y);

    ret->SetBinContent(t+1, th / basearc);
  }
  return ret;
}

//_____________________________________________________________________
Float_t 
AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
{
  TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
  return acc->GetBinContent(t+1);

#if 0
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };
  const Double_t* c1   = (r == 'I' ? ic1      : oc1);
  const Double_t* c2   = (r == 'I' ? ic2      : oc2);
  Double_t  minR       = (r == 'I' ?   4.5213 :  15.4);
  Double_t  maxR       = (r == 'I' ?  17.2    :  28.0);
  Int_t     nStrips    = (r == 'I' ? 512      : 256);
  Int_t     nSec       = (r == 'I' ?  20      :  40);
  Float_t   basearc    = 2 * TMath::Pi() / nSec;
  Double_t  rad        = maxR - minR;
  Float_t   segment    = rad / nStrips;
  Float_t   radius     = minR + t * segment;

  // Old method - calculate full strip area and take ratio to extended
  // strip area
  Float_t   baselen    = basearc * radius;
  Float_t   slope      = (c1[1] - c2[1]) / (c1[0] - c2[0]);
  Float_t   constant   = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
  Float_t   d          = (TMath::Power(TMath::Abs(radius*slope),2) + 
			  TMath::Power(radius,2) - TMath::Power(constant,2));

  // If below corners return 1
  if (d >= 0) return 1;
 
  Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
			 (1+TMath::Power(slope, 2)));
  Float_t   y         = slope*x + constant;

  // If x is larger than corner x or y less than corner y, we have full
  // length strip
  if(x >= c1[0] || y <= c1[1]) return 1;

  //One sector since theta is by definition half-hybrid
  Float_t   theta     = TMath::ATan2(x,y);
  Float_t   arclen    = radius * theta;
  
  // Calculate the area of a strip with no cut
  Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
  Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
  Float_t   basearea  = basearea1 - basearea2;

  // Calculate the area of a strip with cut
  Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
  Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
  Float_t   area      = area1 - area2;
  
  // Acceptance is ratio 
  return area/basearea;
#endif
}

//____________________________________________________________________
void
AliFMDDensityCalculator::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
AliFMDDensityCalculator::DefineOutput(TList* dir)
{
  TList* d = new TList;
  d->SetName(GetName());
  dir->Add(d);
  d->Add(fWeightedSum);
  d->Add(fSumOfWeights);
  d->Add(fCorrections);
  d->Add(fAccI);
  d->Add(fAccO);

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->Output(d);
  }
}
//____________________________________________________________________
void
AliFMDDensityCalculator::Print(Option_t*) const
{
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << "AliFMDDensityCalculator: " << GetName() << '\n'
	    << ind << " Multiplicity cut:       " << fMultCut << '\n'
	    << ind << " Max(particles):         " << fMaxParticles 
	    << std::endl;
}

//====================================================================
AliFMDDensityCalculator::RingHistos::RingHistos()
  : AliForwardUtil::RingHistos(),
    fEvsN(0), 
    fEvsM(0), 
    fEtaVsN(0),
    fEtaVsM(0),
    fCorr(0),
    fDensity(0)
{}

//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r),
    fEvsN(0), 
    fEvsM(0),
    fEtaVsN(0),
    fEtaVsM(0),
    fCorr(0),
    fDensity(0)
{
  fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), 
		   Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
			"N_{ch} in %s", fName.Data()), 
		   2500, -.5, 24.5, 2500, -.5, 24.5);
  fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), 
		   Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
			"N_{ch} in %s", fName.Data()), 
		   2500, -.5, 24.5, 2500, -.5, 24.5);
  fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
  fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
  fEvsN->Sumw2();
  fEvsN->SetDirectory(0);
  fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
  fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
  fEvsM->Sumw2();
  fEvsM->SetDirectory(0);

  fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
			 Form("Average inclusive N_{ch} vs #eta (uncorrected) "
			      "in %s", fName.Data()), 200, -4, 6);
  fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
			 Form("Average inclusive N_{ch} vs #eta (corrected) "
			      "in %s", fName.Data()), 200, -4, 6);
  fEtaVsN->SetXTitle("#eta");
  fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
  fEtaVsN->SetDirectory(0);
  fEtaVsN->SetLineColor(Color());
  fEtaVsN->SetFillColor(Color());
  fEtaVsM->SetXTitle("#eta");
  fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
  fEtaVsM->SetDirectory(0);
  fEtaVsM->SetLineColor(Color());
  fEtaVsM->SetFillColor(Color());


  fCorr = new TProfile(Form("%s_corr", fName.Data()),
			 Form("Average correction in %s", fName.Data()), 
		       200, -4, 6);
  fCorr->SetXTitle("#eta");
  fCorr->SetYTitle("#LT correction#GT");
  fCorr->SetDirectory(0);
  fCorr->SetLineColor(Color());
  fCorr->SetFillColor(Color());

  fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), 
		      Form("Inclusive N_{ch} density in %s", fName.Data()), 
		      200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
		      0, 2*TMath::Pi());
  fDensity->SetDirectory(0);
  fDensity->SetXTitle("#eta");
  fDensity->SetYTitle("#phi [radians]");
  fDensity->SetZTitle("Inclusive N_{ch} density");
}
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
  : AliForwardUtil::RingHistos(o), 
    fEvsN(o.fEvsN), 
    fEvsM(o.fEvsM),
    fEtaVsN(o.fEtaVsN),
    fEtaVsM(o.fEtaVsM),
    fCorr(o.fCorr),
    fDensity(o.fDensity)
{}

//____________________________________________________________________
AliFMDDensityCalculator::RingHistos&
AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
{
  AliForwardUtil::RingHistos::operator=(o);
  
  if (fEvsN)    delete  fEvsN;
  if (fEvsM)    delete  fEvsM;
  if (fEtaVsN)  delete  fEtaVsN;
  if (fEtaVsM)  delete  fEtaVsM;
  if (fCorr)    delete  fCorr;
  if (fDensity) delete fDensity;
  
  fEvsN    = static_cast<TH2D*>(o.fEvsN->Clone());
  fEvsM    = static_cast<TH2D*>(o.fEvsM->Clone());
  fEtaVsN  = static_cast<TProfile*>(o.fEtaVsN->Clone());
  fEtaVsM  = static_cast<TProfile*>(o.fEtaVsM->Clone());
  fCorr    = static_cast<TProfile*>(o.fCorr->Clone());
  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
  
  return *this;
}
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::~RingHistos()
{
  if (fEvsN)    delete fEvsN;
  if (fEvsM)    delete fEvsM;
  if (fEtaVsN)  delete fEtaVsN;
  if (fEtaVsM)  delete fEtaVsM;
  if (fCorr)    delete fCorr;
  if (fDensity) delete fDensity;
}

//____________________________________________________________________
void
AliFMDDensityCalculator::RingHistos::Output(TList* dir)
{
  TList* d = DefineOutputList(dir);
  d->Add(fEvsN);
  d->Add(fEvsM);
  d->Add(fEtaVsN);
  d->Add(fEtaVsM);
  d->Add(fCorr);
  d->Add(fDensity);
}

//____________________________________________________________________
void
AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
{
  TList* l = GetOutputList(dir);
  if (!l) return; 

  TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
  if (density) density->Scale(1./nEvents);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7e4038b5	1	#include "AliFMDDensityCalculator.h"
	2	#include <AliESDFMD.h>
	3	#include <TAxis.h>
	4	#include <TList.h>
	5	#include <TMath.h>
0bd4b00f	6	#include "AliForwardCorrectionManager.h"
0bd4b00f	7	// #include "AliFMDAnaParameters.h"
7e4038b5	8	#include "AliLog.h"
7e4038b5	9	#include <TH2D.h>
0bd4b00f	10	#include <TProfile.h>
	11	#include <TROOT.h>
	12	#include <iostream>
	13	#include <iomanip>
7e4038b5	14
	15	ClassImp(AliFMDDensityCalculator)
	16	#if 0
	17	; // For Emacs
	18	#endif
	19
	20	//____________________________________________________________________
	21	AliFMDDensityCalculator::AliFMDDensityCalculator()
	22	: TNamed(),
	23	fRingHistos(),
0bd4b00f	24	fMultCut(0),
dd497217	25	fSumOfWeights(0),
dd497217	26	fWeightedSum(0),
ea3e5d95	27	fCorrections(0),
0bd4b00f	28	fMaxParticles(5),
	29	fAccI(0),
	30	fAccO(0),
ea3e5d95	31	fDebug(0)
7e4038b5	32	{}
	33
	34	//____________________________________________________________________
	35	AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
	36	: TNamed("fmdDensityCalculator", title),
	37	fRingHistos(),
0bd4b00f	38	fMultCut(0),
dd497217	39	fSumOfWeights(0),
dd497217	40	fWeightedSum(0),
ea3e5d95	41	fCorrections(0),
0bd4b00f	42	fMaxParticles(5),
	43	fAccI(0),
	44	fAccO(0),
ea3e5d95	45	fDebug(0)
7e4038b5	46	{
	47	fRingHistos.SetName(GetName());
	48	fRingHistos.Add(new RingHistos(1, 'I'));
	49	fRingHistos.Add(new RingHistos(2, 'I'));
	50	fRingHistos.Add(new RingHistos(2, 'O'));
	51	fRingHistos.Add(new RingHistos(3, 'I'));
	52	fRingHistos.Add(new RingHistos(3, 'O'));
dd497217	53	fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
dd497217	54	200, 0, 20);
0bd4b00f	55	fSumOfWeights->SetFillColor(kRed+1);
0bd4b00f	56	fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
dd497217	57	fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
dd497217	58	200, 0, 20);
0bd4b00f	59	fWeightedSum->SetFillColor(kBlue+1);
0bd4b00f	60	fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
dd497217	61	fCorrections = new TH1D("corrections", "Distribution of corrections",
dd497217	62	100, 0, 10);
0bd4b00f	63	fCorrections->SetFillColor(kBlue+1);
0bd4b00f	64	fCorrections->SetXTitle("correction");
dd497217	65
0bd4b00f	66	fAccI = GenerateAcceptanceCorrection('I');
0bd4b00f	67	fAccO = GenerateAcceptanceCorrection('O');
7e4038b5	68	}
	69
	70	//____________________________________________________________________
	71	AliFMDDensityCalculator::AliFMDDensityCalculator(const
	72	AliFMDDensityCalculator& o)
	73	: TNamed(o),
	74	fRingHistos(),
dd497217	75	fMultCut(o.fMultCut),
	76	fSumOfWeights(o.fSumOfWeights),
	77	fWeightedSum(o.fWeightedSum),
ea3e5d95	78	fCorrections(o.fCorrections),
0bd4b00f	79	fMaxParticles(o.fMaxParticles),
	80	fAccI(o.fAccI),
	81	fAccO(o.fAccO),
ea3e5d95	82	fDebug(o.fDebug)
7e4038b5	83	{
	84	TIter next(&o.fRingHistos);
	85	TObject* obj = 0;
	86	while ((obj = next())) fRingHistos.Add(obj);
	87	}
	88
	89	//____________________________________________________________________
	90	AliFMDDensityCalculator::~AliFMDDensityCalculator()
	91	{
	92	fRingHistos.Delete();
	93	}
	94
	95	//____________________________________________________________________
	96	AliFMDDensityCalculator&
	97	AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
	98	{
ea3e5d95	99	TNamed::operator=(o);
7e4038b5	100
0bd4b00f	101	fMultCut = o.fMultCut;
	102	fDebug = o.fDebug;
	103	fMaxParticles = o.fMaxParticles;
	104	fAccI = o.fAccI;
	105	fAccO = o.fAccO;
7e4038b5	106
	107	fRingHistos.Delete();
	108	TIter next(&o.fRingHistos);
	109	TObject* obj = 0;
	110	while ((obj = next())) fRingHistos.Add(obj);
	111
	112	return *this;
	113	}
	114
	115	//____________________________________________________________________
	116	AliFMDDensityCalculator::RingHistos*
	117	AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
	118	{
	119	Int_t idx = -1;
	120	switch (d) {
	121	case 1: idx = 0; break;
	122	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	123	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	124	}
	125	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
	126
	127	return static_cast<RingHistos*>(fRingHistos.At(idx));
	128	}
	129
0bd4b00f	130	//____________________________________________________________________
	131	Double_t
	132	AliFMDDensityCalculator::GetMultCut() const
	133	{
	134	if (fMultCut > 0) return fMultCut;
	135
	136	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	137	AliFMDCorrELossFit* fits = fcm.GetELossFit();
	138	return fits->GetLowCut();
	139	}
	140
7e4038b5	141	//____________________________________________________________________
	142	Bool_t
	143	AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
	144	AliForwardUtil::Histos& hists,
0bd4b00f	145	UShort_t vtxbin,
7e4038b5	146	Bool_t lowFlux)
	147	{
	148	for (UShort_t d=1; d<=3; d++) {
	149	UShort_t nr = (d == 1 ? 1 : 2);
	150	for (UShort_t q=0; q<nr; q++) {
	151	Char_t r = (q == 0 ? 'I' : 'O');
	152	UShort_t ns= (q == 0 ? 20 : 40);
	153	UShort_t nt= (q == 0 ? 512 : 256);
	154	TH2D* h = hists.Get(d,r);
	155	RingHistos* rh= GetRingHistos(d,r);
	156
	157	for (UShort_t s=0; s<ns; s++) {
	158	for (UShort_t t=0; t<nt; t++) {
	159	Float_t mult = fmd.Multiplicity(d,r,s,t);
	160
	161	if (mult == 0 \|\| mult > 20) continue;
	162
	163	Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
	164	Float_t eta = fmd.Eta(d,r,s,t);
	165
	166	Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
	167	rh->fEvsN->Fill(mult,n);
0bd4b00f	168	rh->fEtaVsN->Fill(eta, n);
7e4038b5	169
7e4038b5	170	Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
dd497217	171	fCorrections->Fill(c);
7e4038b5	172	if (c > 0) n /= c;
7e4038b5	173	rh->fEvsM->Fill(mult,n);
0bd4b00f	174	rh->fEtaVsM->Fill(eta, n);
0bd4b00f	175	rh->fCorr->Fill(eta, c);
7e4038b5	176
	177	h->Fill(eta,phi,n);
	178	rh->fDensity->Fill(eta,phi,n);
	179	} // for t
	180	} // for s
	181	} // for q
	182	} // for d
	183
	184	return kTRUE;
	185	}
	186
	187	//_____________________________________________________________________
	188	Float_t
	189	AliFMDDensityCalculator::NParticles(Float_t mult,
	190	UShort_t d,
	191	Char_t r,
	192	UShort_t /s/,
	193	UShort_t /t/,
0bd4b00f	194	UShort_t /v/,
7e4038b5	195	Float_t eta,
	196	Bool_t lowFlux) const
	197	{
0bd4b00f	198	if (mult <= GetMultCut()) return 0;
7e4038b5	199	if (lowFlux) return 1;
0bd4b00f	200
	201	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	202	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	203	AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
	204	if (!fit) {
	205	AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
	206	return 0;
	207	}
	208
	209	Int_t m = fit->FindMaxWeight();
	210	if (m < 1) {
	211	AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
	212	return 0;
	213	}
	214	UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
	215	Double_t ret = fit->EvaluateWeighted(mult, n);
	216
	217	if (fDebug > 10) {
	218	AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
	219	}
7e4038b5	220
0bd4b00f	221	fWeightedSum->Fill(ret);
0bd4b00f	222	fSumOfWeights->Fill(ret);
7e4038b5	223
0bd4b00f	224	return ret;
0bd4b00f	225	#if 0
7e4038b5	226	Float_t mpv = pars->GetMPV(d,r,eta);
	227	Float_t w = pars->GetSigma(d,r,eta);
	228	Float_t w2 = pars->Get2MIPWeight(d,r,eta);
	229	Float_t w3 = pars->Get3MIPWeight(d,r,eta);
	230	Float_t mpv2 = 2mpv+2w*TMath::Log(2);
	231	Float_t mpv3 = 3mpv+3w*TMath::Log(3);
	232
	233	Float_t sum = (TMath::Landau(mult,mpv,w,kTRUE) +
	234	w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) +
	235	w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
	236	Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
	237	2w2 TMath::Landau(mult,mpv2,2*w,kTRUE) +
	238	3w3 TMath::Landau(mult,mpv3,3*w,kTRUE));
	239
dd497217	240	fWeightedSum->Fill(wsum);
dd497217	241	fSumOfWeights->Fill(sum);
7e4038b5	242	return (sum > 0) ? wsum / sum : 1;
0bd4b00f	243	#endif
7e4038b5	244	}
	245
	246	//_____________________________________________________________________
	247	Float_t
0bd4b00f	248	AliFMDDensityCalculator::Correction(UShort_t d,
	249	Char_t r,
	250	UShort_t /s/,
	251	UShort_t t,
	252	UShort_t /v/,
	253	Float_t eta,
	254	Bool_t lowFlux) const
7e4038b5	255	{
0bd4b00f	256	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
0bd4b00f	257	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5	258
	259	Float_t correction = AcceptanceCorrection(r,t);
	260	if (lowFlux) {
0bd4b00f	261	// TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
	262	TH1D* dblHitCor = 0;
	263	if (fcm.GetDoubleHit())
	264	dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
	265
7e4038b5	266	if (dblHitCor) {
0bd4b00f	267	Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
0bd4b00f	268	if (dblC > 0) correction *= dblC;
7e4038b5	269	}
0bd4b00f	270	else {
7e4038b5	271	AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
0bd4b00f	272	}
7e4038b5	273	}
7e4038b5	274
0bd4b00f	275	#if 0
7e4038b5	276	TH1F* deadCor = pars->GetFMDDeadCorrection(v);
	277	if (deadCor) {
	278	Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
	279	if (deadC > 0)
	280	correction *= deadC;
	281	}
0bd4b00f	282	#endif
0bd4b00f	283
7e4038b5	284	return correction;
	285	}
	286
0bd4b00f	287	//_____________________________________________________________________
	288	TH1D*
	289	AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
	290	{
	291	const Double_t ic1[] = { 4.9895, 15.3560 };
	292	const Double_t ic2[] = { 1.8007, 17.2000 };
	293	const Double_t oc1[] = { 4.2231, 26.6638 };
	294	const Double_t oc2[] = { 1.8357, 27.9500 };
	295	const Double_t* c1 = (r == 'I' \|\| r == 'i' ? ic1 : oc1);
	296	const Double_t* c2 = (r == 'I' \|\| r == 'i' ? ic2 : oc2);
	297	Double_t minR = (r == 'I' \|\| r == 'i' ? 4.5213 : 15.4);
	298	Double_t maxR = (r == 'I' \|\| r == 'i' ? 17.2 : 28.0);
	299	Int_t nStrips = (r == 'I' \|\| r == 'i' ? 512 : 256);
	300	Int_t nSec = (r == 'I' \|\| r == 'i' ? 20 : 40);
	301	Float_t basearc = 2 * TMath::Pi() / nSec;
	302	Double_t rad = maxR - minR;
	303	Float_t segment = rad / nStrips;
	304	Float_t cr = TMath::Sqrt(c1[0]c1[0]+c1[1]c1[1]);
	305
	306	// Numbers used to find end-point of strip.
	307	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
	308	Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
	309	Float_t dx = c2[0] - c1[0];
	310	Float_t dy = c2[1] - c1[1];
	311	Float_t dr = TMath::Sqrt(dxdx+dydy);
	312
	313	TH1D* ret = new TH1D(Form("acc%c", r),
	314	Form("Acceptance correction for FMDx%c", r),
	315	nStrips, -.5, nStrips-.5);
	316	ret->SetXTitle("Strip");
	317	ret->SetYTitle("#varphi acceptance");
	318	ret->SetDirectory(0);
	319	ret->SetFillColor(r == 'I' \|\| r == 'i' ? kRed+1 : kBlue+1);
	320	ret->SetFillStyle(3001);
	321
	322	for (Int_t t = 0; t < nStrips; t++) {
	323	Float_t radius = minR + t * segment;
	324
	325	// If the radius of the strip is smaller than the radius corresponding
	326	// to the first corner we have a full strip length
	327	if (radius <= cr) {
	328	ret->SetBinContent(t+1, 1);
	329	continue;
	330	}
	331
	332	// Next, we should find the end-point of the strip - that is,
	333	// the coordinates where the line from c1 to c2 intersects a circle
	334	// with radius given by the strip.
	335	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
	336	// Calculate the determinant
	337	Float_t det = radius * radius * dr * dr - D*D;
	338
	339	if (det <= 0) {
	340	// <0 means No intersection
	341	// =0 means Exactly tangent
	342	ret->SetBinContent(t+1, 1);
	343	continue;
	344	}
	345
	346	// Calculate end-point and the corresponding opening angle
	347	Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
	348	Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
	349	Float_t th = TMath::ATan2(x, y);
	350
351	ret->SetBinContent(t+1, th / basearc);
352	}
353	return ret;
354	}
7e4038b5	355
	356	//_____________________________________________________________________
	357	Float_t
	358	AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
	359	{
0bd4b00f	360	TH1D* acc = (r == 'I' \|\| r == 'i' ? fAccI : fAccO);
	361	return acc->GetBinContent(t+1);
	362
	363	#if 0
	364	const Double_t ic1[] = { 4.9895, 15.3560 };
	365	const Double_t ic2[] = { 1.8007, 17.2000 };
	366	const Double_t oc1[] = { 4.2231, 26.6638 };
	367	const Double_t oc2[] = { 1.8357, 27.9500 };
	368	const Double_t* c1 = (r == 'I' ? ic1 : oc1);
	369	const Double_t* c2 = (r == 'I' ? ic2 : oc2);
	370	Double_t minR = (r == 'I' ? 4.5213 : 15.4);
	371	Double_t maxR = (r == 'I' ? 17.2 : 28.0);
	372	Int_t nStrips = (r == 'I' ? 512 : 256);
	373	Int_t nSec = (r == 'I' ? 20 : 40);
	374	Float_t basearc = 2 * TMath::Pi() / nSec;
	375	Double_t rad = maxR - minR;
	376	Float_t segment = rad / nStrips;
	377	Float_t radius = minR + t * segment;
	378
	379	// Old method - calculate full strip area and take ratio to extended
	380	// strip area
	381	Float_t baselen = basearc * radius;
	382	Float_t slope = (c1[1] - c2[1]) / (c1[0] - c2[0]);
	383	Float_t constant = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
	384	Float_t d = (TMath::Power(TMath::Abs(radius*slope),2) +
	385	TMath::Power(radius,2) - TMath::Power(constant,2));
	386
	387	// If below corners return 1
	388	if (d >= 0) return 1;
	389
	390	Float_t x = ((-TMath::Sqrt(d) - slope * constant) /
	391	(1+TMath::Power(slope, 2)));
	392	Float_t y = slope*x + constant;
	393
	394	// If x is larger than corner x or y less than corner y, we have full
	395	// length strip
	396	if(x >= c1[0] \|\| y <= c1[1]) return 1;
	397
	398	//One sector since theta is by definition half-hybrid
	399	Float_t theta = TMath::ATan2(x,y);
	400	Float_t arclen = radius * theta;
7e4038b5	401
0bd4b00f	402	// Calculate the area of a strip with no cut
	403	Float_t basearea1 = 0.5 * baselen * TMath::Power(radius,2);
	404	Float_t basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
7e4038b5	405	Float_t basearea = basearea1 - basearea2;
0bd4b00f	406
	407	// Calculate the area of a strip with cut
	408	Float_t area1 = 0.5 * arclen * TMath::Power(radius,2);
	409	Float_t area2 = 0.5 * arclen * TMath::Power((radius-segment),2);
7e4038b5	410	Float_t area = area1 - area2;
7e4038b5	411
0bd4b00f	412	// Acceptance is ratio
7e4038b5	413	return area/basearea;
0bd4b00f	414	#endif
7e4038b5	415	}
	416
	417	//____________________________________________________________________
	418	void
9d99b0dd	419	AliFMDDensityCalculator::ScaleHistograms(TList* dir, Int_t nEvents)
7e4038b5	420	{
7e4038b5	421	if (nEvents <= 0) return;
9d99b0dd	422	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
9d99b0dd	423	if (!d) return;
7e4038b5	424
	425	TIter next(&fRingHistos);
	426	RingHistos* o = 0;
9d99b0dd	427	while ((o = static_cast<RingHistos*>(next())))
9d99b0dd	428	o->ScaleHistograms(d, nEvents);
7e4038b5	429	}
	430
	431	//____________________________________________________________________
	432	void
9d99b0dd	433	AliFMDDensityCalculator::DefineOutput(TList* dir)
7e4038b5	434	{
	435	TList* d = new TList;
	436	d->SetName(GetName());
	437	dir->Add(d);
dd497217	438	d->Add(fWeightedSum);
	439	d->Add(fSumOfWeights);
	440	d->Add(fCorrections);
0bd4b00f	441	d->Add(fAccI);
	442	d->Add(fAccO);
	443
7e4038b5	444	TIter next(&fRingHistos);
	445	RingHistos* o = 0;
	446	while ((o = static_cast<RingHistos*>(next()))) {
	447	o->Output(d);
	448	}
	449	}
0bd4b00f	450	//____________________________________________________________________
	451	void
	452	AliFMDDensityCalculator::Print(Option_t*) const
	453	{
	454	char ind[gROOT->GetDirLevel()+1];
	455	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	456	ind[gROOT->GetDirLevel()] = '\0';
	457	std::cout << ind << "AliFMDDensityCalculator: " << GetName() << '\n'
	458	<< ind << " Multiplicity cut: " << fMultCut << '\n'
	459	<< ind << " Max(particles): " << fMaxParticles
	460	<< std::endl;
	461	}
7e4038b5	462
	463	//====================================================================
	464	AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd	465	: AliForwardUtil::RingHistos(),
7e4038b5	466	fEvsN(0),
7e4038b5	467	fEvsM(0),
0bd4b00f	468	fEtaVsN(0),
	469	fEtaVsM(0),
	470	fCorr(0),
7e4038b5	471	fDensity(0)
	472	{}
	473
	474	//____________________________________________________________________
	475	AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd	476	: AliForwardUtil::RingHistos(d,r),
7e4038b5	477	fEvsN(0),
7e4038b5	478	fEvsM(0),
0bd4b00f	479	fEtaVsN(0),
	480	fEtaVsM(0),
	481	fCorr(0),
7e4038b5	482	fDensity(0)
7e4038b5	483	{
9d99b0dd	484	fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()),
0bd4b00f	485	Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
9d99b0dd	486	"N_{ch} in %s", fName.Data()),
0bd4b00f	487	2500, -.5, 24.5, 2500, -.5, 24.5);
9d99b0dd	488	fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()),
0bd4b00f	489	Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
	490	"N_{ch} in %s", fName.Data()),
	491	2500, -.5, 24.5, 2500, -.5, 24.5);
7e4038b5	492	fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
	493	fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
	494	fEvsN->Sumw2();
	495	fEvsN->SetDirectory(0);
0bd4b00f	496	fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
0bd4b00f	497	fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
7e4038b5	498	fEvsM->Sumw2();
	499	fEvsM->SetDirectory(0);
	500
0bd4b00f	501	fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
	502	Form("Average inclusive N_{ch} vs #eta (uncorrected) "
	503	"in %s", fName.Data()), 200, -4, 6);
	504	fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
	505	Form("Average inclusive N_{ch} vs #eta (corrected) "
	506	"in %s", fName.Data()), 200, -4, 6);
	507	fEtaVsN->SetXTitle("#eta");
	508	fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
	509	fEtaVsN->SetDirectory(0);
	510	fEtaVsN->SetLineColor(Color());
	511	fEtaVsN->SetFillColor(Color());
	512	fEtaVsM->SetXTitle("#eta");
	513	fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
	514	fEtaVsM->SetDirectory(0);
	515	fEtaVsM->SetLineColor(Color());
	516	fEtaVsM->SetFillColor(Color());
	517
	518
	519	fCorr = new TProfile(Form("%s_corr", fName.Data()),
	520	Form("Average correction in %s", fName.Data()),
	521	200, -4, 6);
	522	fCorr->SetXTitle("#eta");
	523	fCorr->SetYTitle("#LT correction#GT");
	524	fCorr->SetDirectory(0);
	525	fCorr->SetLineColor(Color());
	526	fCorr->SetFillColor(Color());
	527
9d99b0dd	528	fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()),
0bd4b00f	529	Form("Inclusive N_{ch} density in %s", fName.Data()),
7e4038b5	530	200, -4, 6, (r == 'I' \|\| r == 'i' ? 20 : 40),
	531	0, 2*TMath::Pi());
	532	fDensity->SetDirectory(0);
	533	fDensity->SetXTitle("#eta");
	534	fDensity->SetYTitle("#phi [radians]");
	535	fDensity->SetZTitle("Inclusive N_{ch} density");
	536	}
	537	//____________________________________________________________________
	538	AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd	539	: AliForwardUtil::RingHistos(o),
7e4038b5	540	fEvsN(o.fEvsN),
7e4038b5	541	fEvsM(o.fEvsM),
0bd4b00f	542	fEtaVsN(o.fEtaVsN),
	543	fEtaVsM(o.fEtaVsM),
	544	fCorr(o.fCorr),
7e4038b5	545	fDensity(o.fDensity)
	546	{}
	547
	548	//____________________________________________________________________
	549	AliFMDDensityCalculator::RingHistos&
	550	AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
	551	{
9d99b0dd	552	AliForwardUtil::RingHistos::operator=(o);
7e4038b5	553
	554	if (fEvsN) delete fEvsN;
	555	if (fEvsM) delete fEvsM;
0bd4b00f	556	if (fEtaVsN) delete fEtaVsN;
	557	if (fEtaVsM) delete fEtaVsM;
	558	if (fCorr) delete fCorr;
7e4038b5	559	if (fDensity) delete fDensity;
	560
	561	fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
	562	fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
0bd4b00f	563	fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
	564	fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
	565	fCorr = static_cast<TProfile*>(o.fCorr->Clone());
7e4038b5	566	fDensity = static_cast<TH2D*>(o.fDensity->Clone());
	567
	568	return *this;
	569	}
	570	//____________________________________________________________________
	571	AliFMDDensityCalculator::RingHistos::~RingHistos()
	572	{
	573	if (fEvsN) delete fEvsN;
	574	if (fEvsM) delete fEvsM;
0bd4b00f	575	if (fEtaVsN) delete fEtaVsN;
	576	if (fEtaVsM) delete fEtaVsM;
	577	if (fCorr) delete fCorr;
7e4038b5	578	if (fDensity) delete fDensity;
	579	}
	580
	581	//____________________________________________________________________
	582	void
	583	AliFMDDensityCalculator::RingHistos::Output(TList* dir)
	584	{
9d99b0dd	585	TList* d = DefineOutputList(dir);
7e4038b5	586	d->Add(fEvsN);
7e4038b5	587	d->Add(fEvsM);
0bd4b00f	588	d->Add(fEtaVsN);
	589	d->Add(fEtaVsM);
	590	d->Add(fCorr);
7e4038b5	591	d->Add(fDensity);
9d99b0dd	592	}
	593
	594	//____________________________________________________________________
	595	void
	596	AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
	597	{
	598	TList* l = GetOutputList(dir);
	599	if (!l) return;
	600
	601	TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
	602	if (density) density->Scale(1./nEvents);
7e4038b5	603	}
	604
	605	//____________________________________________________________________
	606	//
	607	// EOF
	608	//
	609
	610
	611