[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMCCorrectionsTask.cxx

// 
// Calculate the corrections in the central regions
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODCentralMult 
// 
// Histograms 
//   
// Corrections used 
// 
#include "AliCentralMCCorrectionsTask.h"
#include "AliCentralCorrectionManager.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliAODForwardMult.h"
#include "AliCentralCorrSecondaryMap.h"
#include "AliCentralCorrAcceptance.h"
#include "AliForwardUtil.h"
#include "AliMultiplicity.h"
#include <TH1.h>
#include <TH2D.h>
#include <TDirectory.h>
#include <TList.h>
#include <TROOT.h>
#include <iostream>

//====================================================================
namespace {
  const char* GetEventName(Bool_t tr, Bool_t vtx) 
  {
    return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  }
}

//====================================================================
AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
  : AliBaseMCCorrectionsTask(),
    fTrackDensity(),
    fSecCorr(0), 
    fAccCorr(0),
    fNPhiBins(20),
    fEffectiveCorr(true),
    fEtaCut(1.9),
    fCorrCut(0.6)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
}

//____________________________________________________________________
AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
  : AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())),
    fTrackDensity("trackDensity"),
    fSecCorr(0), 
    fAccCorr(0),
    fNPhiBins(20),
    fEffectiveCorr(true),
    fEtaCut(1.9),
    fCorrCut(0.6)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
}

//____________________________________________________________________
AliBaseMCCorrectionsTask::VtxBin*
AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high)
{
  return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins);
}

//____________________________________________________________________
Bool_t
AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd, 
					const AliMCEvent& mc, 
					AliBaseMCCorrectionsTask::VtxBin& bin,
					Double_t          vz)
{
  AliCentralMCCorrectionsTask::VtxBin& vb = 
    static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin);

  // Now process our input data and store in own ESD object 
  fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary);
  
  // Get the ESD object
  const AliMultiplicity* spdmult = esd.GetMultiplicity();
 
  // Count number of tracklets per bin 
  for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) 
    vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
  //...and then the unused clusters in layer 1 
  for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
    Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
    vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j));
  }
  
  // Count events  
  bin.fCounts->Fill(0.5);
  
  return true;
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::CreateCorrections(TList* results)
{
  fSecCorr = new AliCentralCorrSecondaryMap;
  fSecCorr->SetVertexAxis(fVtxAxis);

  fAccCorr = new AliCentralCorrAcceptance;
  fAccCorr->SetVertexAxis(fVtxAxis);

  results->Add(fSecCorr);
  results->Add(fAccCorr);
}

//____________________________________________________________________
Bool_t 
AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin* 
					    bin,  UShort_t iVz) 
{
  
  AliCentralMCCorrectionsTask::VtxBin* vb = 
    static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin);
  vb->Terminate(fList, fResults, iVz, fEffectiveCorr, 
		fEtaCut, fCorrCut, fSecCorr,fAccCorr);
  return true;
}
					    

//____________________________________________________________________
void
AliCentralMCCorrectionsTask::Print(Option_t* option) const
{
  AliBaseMCCorrectionsTask::Print(option);
  std::cout << "  # of phi bins:    " << fNPhiBins << "\n"
	    << "  Effective corr.:  " << fEffectiveCorr << "\n"
	    << "  Eta cut-off:      " << fEtaCut << "\n"
	    << "  Acceptance cut:   " << fCorrCut 
	    << std::endl;
  gROOT->IncreaseDirLevel();
  fTrackDensity.Print(option);
  gROOT->DecreaseDirLevel();
}

//====================================================================
AliCentralMCCorrectionsTask::VtxBin::VtxBin()
  : AliBaseMCCorrectionsTask::VtxBin(),
    fHits(0), 
    fClusters(0)
{
}
//____________________________________________________________________
AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low, 
					    Double_t     high, 
					    const TAxis& axis,
					    UShort_t     nPhi)
  : AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi),
    fHits(0), 
    fClusters(0)
{
  fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
  fHits->SetTitle("Hits");
  fHits->SetDirectory(0);

  fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
  fClusters->SetTitle("Clusters");
  fClusters->SetDirectory(0);
}


//____________________________________________________________________
TList*
AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
{
  TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l);

  d->Add(fHits);
  d->Add(fClusters);

  return d;
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input, 
					       TList* output, 
					       UShort_t iVz, 
					       Bool_t effectiveCorr,
					       Double_t etaCut,
					       Double_t accCut,
					       AliCentralCorrSecondaryMap* map,
					       AliCentralCorrAcceptance* acorr)
{
  TList* out = new TList;
  out->SetName(GetName());
  out->SetOwner();
  output->Add(out);

  TList* l = static_cast<TList*>(input->FindObject(GetName()));
  if (!l) { 
    AliError(Form("List %s not found in %s", GetName(), input->GetName()));
    return;
  }

  // Get the sums 
  TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
  TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
  TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
  if (!hits || !prim) {
    AliError(Form("Missing histograms: %p, %p", hits, prim));
    return;
  }

  // Clone cluster and hit map
  TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
  TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
  secMapEff->SetTitle("2^{nd} map from clusters");
  secMapEff->SetDirectory(0);
  secMapHit->SetTitle("2^{nd} map from MC hits");
  secMapHit->SetDirectory(0);

  // Divide cluster and hit map with number of primaries 
  secMapEff->Divide(prim);
  secMapHit->Divide(prim);

  // Create acceptance histograms 
  TH1D* accEff = new TH1D("accEff",
			  "Acceptance correction for SPD (from 2^{nd} map)" ,
			  fPrimary->GetXaxis()->GetNbins(), 
			  fPrimary->GetXaxis()->GetXmin(), 
			  fPrimary->GetXaxis()->GetXmax());
  TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
  accHit->SetTitle("Acceptance correction for SPD (from clusters)");

  // Diagnostics histogra, 
  TH2*  dia    = static_cast<TH2D*>(clus->Clone("diagnostics"));
  dia->SetTitle("Scaled cluster density");

  // Total number of channels along phi and # of eta bins
  Int_t nTotal = secMapHit->GetNbinsY();
  Int_t nEta   = secMapHit->GetNbinsX();

  for(Int_t xx = 1; xx <= nEta; xx++) {
    Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
    Bool_t   ins = TMath::Abs(eta) <= etaCut;
    Double_t mm  = 0;
    if (ins) {
      // Find the maximum cluster signal in this phi range 
      for (Int_t yy = 1; yy <= nTotal; yy++) { 
	Double_t c = clus->GetBinContent(xx,yy);
	mm         = TMath::Max(mm, c);
      }
    }
    // Count number of phi bins with enough clusters or high enough
    // correction.
    Int_t nOKEff    = 0;
    Int_t nOKHit    = 0;
    for(Int_t yy = 1; yy <=nTotal; yy++) {
      if (!ins) { // Not inside Eta cut
	secMapEff->SetBinContent(xx,yy,0.); 
	secMapEff->SetBinError(xx,yy,0.); 
	secMapHit->SetBinContent(xx,yy,0.); 
	secMapHit->SetBinError(xx,yy,0.); 
	dia->SetBinContent(xx,yy,0);
	continue;
      }

      // Check if the background correction is large enough, or zero map
      if(secMapEff->GetBinContent(xx,yy) > 0.9) {
	// acc->Fill(h->GetXaxis()->GetBinCenter(xx));
	nOKEff++;
      }
      else {
	secMapEff->SetBinContent(xx,yy,0.); 
	secMapEff->SetBinError(xx,yy,0.); 
      }

      // Check if the number of cluster is large enough, or zero map
      Double_t c = clus->GetBinContent(xx,yy);
      Double_t s = (mm < 1e-6) ? 0 : c / mm;
      dia->SetBinContent(xx,yy,s);
      if (s >= accCut) {
	nOKHit++;
      }
      else {
	secMapHit->SetBinContent(xx,yy,0);
	secMapHit->SetBinError(xx,yy,0);
      }
    }

    // Calculate acceptance as ratio of bins with enough clusters and
    // total number of phi bins.
    Double_t accXX = float(nOKHit) / nTotal;
    if (accXX < 0.2) accXX = 0;
    accHit->SetBinContent(xx, accXX);

    // Calculate acceptance as ratio of bins with large enough
    // correction and total number of phi bins.
    accXX = float(nOKEff) / nTotal;
    if (accXX < 0.2) accXX = 0;
    accEff->SetBinContent(xx, accXX);
  }

  TH2D* secMap    = (effectiveCorr ? secMapEff : secMapHit);
  TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
  TH1D* acc       = (effectiveCorr ? accEff    : accHit);
  TH1D* accAlt    = (effectiveCorr ? accHit    : accEff);
  out->Add(secMap->Clone("secMap"));
  out->Add(secMapAlt->Clone());
  out->Add(acc->Clone("acc"));
  out->Add(accAlt->Clone());
  out->Add(dia->Clone());

  map->SetCorrection(iVz, secMap);
  acorr->SetCorrection(iVz, acc);
}

//
// EOF
//
Commit	Line	Data
600b313f	1	//
	2	// Calculate the corrections in the central regions
	3	//
	4	// Inputs:
	5	// - AliESDEvent
	6	//
	7	// Outputs:
	8	// - AliAODCentralMult
	9	//
	10	// Histograms
	11	//
	12	// Corrections used
	13	//
	14	#include "AliCentralMCCorrectionsTask.h"
c8b1a7db	15	#include "AliCentralCorrectionManager.h"
600b313f	16	#include "AliTriggerAnalysis.h"
	17	#include "AliPhysicsSelection.h"
	18	#include "AliLog.h"
	19	#include "AliHeader.h"
	20	#include "AliGenEventHeader.h"
	21	#include "AliESDEvent.h"
	22	#include "AliAODHandler.h"
	23	#include "AliMultiplicity.h"
	24	#include "AliInputEventHandler.h"
	25	#include "AliStack.h"
	26	#include "AliMCEvent.h"
	27	#include "AliAODForwardMult.h"
	28	#include "AliCentralCorrSecondaryMap.h"
ea6f96fb	29	#include "AliCentralCorrAcceptance.h"
f53fb4f6	30	#include "AliForwardUtil.h"
ea6f96fb	31	#include "AliMultiplicity.h"
600b313f	32	#include <TH1.h>
	33	#include <TH2D.h>
	34	#include <TDirectory.h>
	35	#include <TList.h>
	36	#include <TROOT.h>
	37	#include <iostream>
	38
	39	//====================================================================
	40	namespace {
	41	const char* GetEventName(Bool_t tr, Bool_t vtx)
	42	{
	43	return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
	44	}
600b313f	45	}
	46
	47	//====================================================================
	48	AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
c8b1a7db	49	: AliBaseMCCorrectionsTask(),
600b313f	50	fTrackDensity(),
c8b1a7db	51	fSecCorr(0),
c8b1a7db	52	fAccCorr(0),
ea6f96fb	53	fNPhiBins(20),
bfab35d9	54	fEffectiveCorr(true),
	55	fEtaCut(1.9),
	56	fCorrCut(0.6)
600b313f	57	{
	58	//
	59	// Constructor
	60	//
	61	// Parameters:
	62	// name Name of task
	63	//
5ca83fee	64	DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
600b313f	65	}
	66
	67	//____________________________________________________________________
	68	AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
c8b1a7db	69	: AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())),
600b313f	70	fTrackDensity("trackDensity"),
c8b1a7db	71	fSecCorr(0),
c8b1a7db	72	fAccCorr(0),
ea6f96fb	73	fNPhiBins(20),
bfab35d9	74	fEffectiveCorr(true),
	75	fEtaCut(1.9),
	76	fCorrCut(0.6)
600b313f	77	{
	78	//
	79	// Constructor
	80	//
	81	// Parameters:
	82	// name Name of task
	83	//
5ca83fee	84	DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
600b313f	85	}
	86
	87	//____________________________________________________________________
c8b1a7db	88	AliBaseMCCorrectionsTask::VtxBin*
c8b1a7db	89	AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high)
600b313f	90	{
c8b1a7db	91	return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins);
600b313f	92	}
	93
	94	//____________________________________________________________________
c8b1a7db	95	Bool_t
	96	AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd,
	97	const AliMCEvent& mc,
	98	AliBaseMCCorrectionsTask::VtxBin& bin,
	99	Double_t vz)
600b313f	100	{
c8b1a7db	101	AliCentralMCCorrectionsTask::VtxBin& vb =
c8b1a7db	102	static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin);
600b313f	103
600b313f	104	// Now process our input data and store in own ESD object
c8b1a7db	105	fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary);
c8b1a7db	106
ea6f96fb	107	// Get the ESD object
c8b1a7db	108	const AliMultiplicity* spdmult = esd.GetMultiplicity();
c8b1a7db	109
ea6f96fb	110	// Count number of tracklets per bin
ea6f96fb	111	for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
c8b1a7db	112	vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
ea6f96fb	113	//...and then the unused clusters in layer 1
ea6f96fb	114	for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
c8b1a7db	115	Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
c8b1a7db	116	vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j));
ea6f96fb	117	}
c8b1a7db	118
ea6f96fb	119	// Count events
c8b1a7db	120	bin.fCounts->Fill(0.5);
bcd522ff	121
c8b1a7db	122	return true;
600b313f	123	}
600b313f	124	//____________________________________________________________________
600b313f	125	void
c8b1a7db	126	AliCentralMCCorrectionsTask::CreateCorrections(TList* results)
600b313f	127	{
c8b1a7db	128	fSecCorr = new AliCentralCorrSecondaryMap;
c8b1a7db	129	fSecCorr->SetVertexAxis(fVtxAxis);
600b313f	130
c8b1a7db	131	fAccCorr = new AliCentralCorrAcceptance;
c8b1a7db	132	fAccCorr->SetVertexAxis(fVtxAxis);
ea6f96fb	133
c8b1a7db	134	results->Add(fSecCorr);
	135	results->Add(fAccCorr);
	136	}
600b313f	137
c8b1a7db	138	//____________________________________________________________________
	139	Bool_t
	140	AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin*
	141	bin, UShort_t iVz)
	142	{
	143
	144	AliCentralMCCorrectionsTask::VtxBin* vb =
	145	static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin);
	146	vb->Terminate(fList, fResults, iVz, fEffectiveCorr,
	147	fEtaCut, fCorrCut, fSecCorr,fAccCorr);
	148	return true;
600b313f	149	}
c8b1a7db	150
600b313f	151
	152	//____________________________________________________________________
	153	void
	154	AliCentralMCCorrectionsTask::Print(Option_t* option) const
	155	{
c8b1a7db	156	AliBaseMCCorrectionsTask::Print(option);
c8b1a7db	157	std::cout << " # of phi bins: " << fNPhiBins << "\n"
bfab35d9	158	<< " Effective corr.: " << fEffectiveCorr << "\n"
	159	<< " Eta cut-off: " << fEtaCut << "\n"
	160	<< " Acceptance cut: " << fCorrCut
600b313f	161	<< std::endl;
600b313f	162	gROOT->IncreaseDirLevel();
600b313f	163	fTrackDensity.Print(option);
	164	gROOT->DecreaseDirLevel();
	165	}
	166
	167	//====================================================================
600b313f	168	AliCentralMCCorrectionsTask::VtxBin::VtxBin()
c8b1a7db	169	: AliBaseMCCorrectionsTask::VtxBin(),
	170	fHits(0),
	171	fClusters(0)
600b313f	172	{
	173	}
	174	//____________________________________________________________________
	175	AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
	176	Double_t high,
	177	const TAxis& axis,
	178	UShort_t nPhi)
c8b1a7db	179	: AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi),
600b313f	180	fHits(0),
c8b1a7db	181	fClusters(0)
600b313f	182	{
600b313f	183	fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
	184	fHits->SetTitle("Hits");
	185	fHits->SetDirectory(0);
	186
ea6f96fb	187	fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
	188	fClusters->SetTitle("Clusters");
	189	fClusters->SetDirectory(0);
600b313f	190	}
600b313f	191
600b313f	192
600b313f	193	//____________________________________________________________________
c8b1a7db	194	TList*
5934a3e3	195	AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
600b313f	196	{
c8b1a7db	197	TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l);
600b313f	198
600b313f	199	d->Add(fHits);
ea6f96fb	200	d->Add(fClusters);
c8b1a7db	201
c8b1a7db	202	return d;
600b313f	203	}
	204	//____________________________________________________________________
	205	void
5934a3e3	206	AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input,
bfab35d9	207	TList* output,
	208	UShort_t iVz,
	209	Bool_t effectiveCorr,
	210	Double_t etaCut,
	211	Double_t accCut,
	212	AliCentralCorrSecondaryMap* map,
	213	AliCentralCorrAcceptance* acorr)
600b313f	214	{
	215	TList* out = new TList;
	216	out->SetName(GetName());
	217	out->SetOwner();
	218	output->Add(out);
bfab35d9	219
600b313f	220	TList* l = static_cast<TList*>(input->FindObject(GetName()));
	221	if (!l) {
	222	AliError(Form("List %s not found in %s", GetName(), input->GetName()));
	223	return;
	224	}
	225
bfab35d9	226	// Get the sums
600b313f	227	TH2D* hits = static_cast<TH2D*>(l->FindObject("hits"));
ea6f96fb	228	TH2D* clus = static_cast<TH2D*>(l->FindObject("clusters"));
600b313f	229	TH2D* prim = static_cast<TH2D*>(l->FindObject("primary"));
	230	if (!hits \|\| !prim) {
	231	AliError(Form("Missing histograms: %p, %p", hits, prim));
	232	return;
	233	}
	234
bfab35d9	235	// Clone cluster and hit map
	236	TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
	237	TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
	238	secMapEff->SetTitle("2^{nd} map from clusters");
	239	secMapEff->SetDirectory(0);
	240	secMapHit->SetTitle("2^{nd} map from MC hits");
	241	secMapHit->SetDirectory(0);
	242
	243	// Divide cluster and hit map with number of primaries
	244	secMapEff->Divide(prim);
	245	secMapHit->Divide(prim);
	246
	247	// Create acceptance histograms
	248	TH1D* accEff = new TH1D("accEff",
	249	"Acceptance correction for SPD (from 2^{nd} map)" ,
	250	fPrimary->GetXaxis()->GetNbins(),
	251	fPrimary->GetXaxis()->GetXmin(),
	252	fPrimary->GetXaxis()->GetXmax());
	253	TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
	254	accHit->SetTitle("Acceptance correction for SPD (from clusters)");
	255
	256	// Diagnostics histogra,
	257	TH2* dia = static_cast<TH2D*>(clus->Clone("diagnostics"));
	258	dia->SetTitle("Scaled cluster density");
	259
	260	// Total number of channels along phi and # of eta bins
	261	Int_t nTotal = secMapHit->GetNbinsY();
	262	Int_t nEta = secMapHit->GetNbinsX();
	263
	264	for(Int_t xx = 1; xx <= nEta; xx++) {
	265	Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
	266	Bool_t ins = TMath::Abs(eta) <= etaCut;
	267	Double_t mm = 0;
	268	if (ins) {
	269	// Find the maximum cluster signal in this phi range
	270	for (Int_t yy = 1; yy <= nTotal; yy++) {
	271	Double_t c = clus->GetBinContent(xx,yy);
	272	mm = TMath::Max(mm, c);
	273	}
	274	}
	275	// Count number of phi bins with enough clusters or high enough
	276	// correction.
	277	Int_t nOKEff = 0;
	278	Int_t nOKHit = 0;
	279	for(Int_t yy = 1; yy <=nTotal; yy++) {
	280	if (!ins) { // Not inside Eta cut
	281	secMapEff->SetBinContent(xx,yy,0.);
	282	secMapEff->SetBinError(xx,yy,0.);
	283	secMapHit->SetBinContent(xx,yy,0.);
	284	secMapHit->SetBinError(xx,yy,0.);
	285	dia->SetBinContent(xx,yy,0);
	286	continue;
	287	}
	288
	289	// Check if the background correction is large enough, or zero map
	290	if(secMapEff->GetBinContent(xx,yy) > 0.9) {
	291	// acc->Fill(h->GetXaxis()->GetBinCenter(xx));
	292	nOKEff++;
	293	}
	294	else {
	295	secMapEff->SetBinContent(xx,yy,0.);
	296	secMapEff->SetBinError(xx,yy,0.);
	297	}
	298
299	// Check if the number of cluster is large enough, or zero map
300	Double_t c = clus->GetBinContent(xx,yy);
301	Double_t s = (mm < 1e-6) ? 0 : c / mm;
302	dia->SetBinContent(xx,yy,s);
303	if (s >= accCut) {
304	nOKHit++;
ea6f96fb	305	}
ea6f96fb	306	else {
bfab35d9	307	secMapHit->SetBinContent(xx,yy,0);
bfab35d9	308	secMapHit->SetBinError(xx,yy,0);
ea6f96fb	309	}
ea6f96fb	310	}
bfab35d9	311
	312	// Calculate acceptance as ratio of bins with enough clusters and
	313	// total number of phi bins.
	314	Double_t accXX = float(nOKHit) / nTotal;
	315	if (accXX < 0.2) accXX = 0;
	316	accHit->SetBinContent(xx, accXX);
	317
	318	// Calculate acceptance as ratio of bins with large enough
	319	// correction and total number of phi bins.
	320	accXX = float(nOKEff) / nTotal;
	321	if (accXX < 0.2) accXX = 0;
	322	accEff->SetBinContent(xx, accXX);
ea6f96fb	323	}
600b313f	324
bfab35d9	325	TH2D* secMap = (effectiveCorr ? secMapEff : secMapHit);
	326	TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
	327	TH1D* acc = (effectiveCorr ? accEff : accHit);
	328	TH1D* accAlt = (effectiveCorr ? accHit : accEff);
	329	out->Add(secMap->Clone("secMap"));
	330	out->Add(secMapAlt->Clone());
	331	out->Add(acc->Clone("acc"));
	332	out->Add(accAlt->Clone());
	333	out->Add(dia->Clone());
	334
	335	map->SetCorrection(iVz, secMap);
ea6f96fb	336	acorr->SetCorrection(iVz, acc);
600b313f	337	}
	338
	339	//
	340	// EOF
	341	//