[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMult.C

/**
 * @file   UnfoldMult.C
 * @author Valentina Zaccolo
 * @date   Tue Mar 05 12:56:56 2013
 * 
 * @brief  Bayesian Method Unfolding Script
 *
 * @ingroup pwglf_forward_scripts
 * @ingroup pwglf_forward_multdist
 */
#include <iostream>
#include <TList.h>
#include <TFile.h>
#include <TObject.h>
#include <TDirectory.h>
#include "TH1D.h"
#include "TH2D.h"
#include <TROOT.h>
#include <TSystem.h>
#include "RooUnfoldBayes.h"
#include "RooUnfoldResponse.h"
#include <TMath.h>
/** 
 * @defgroup pwglf_forward_scripts_unfoldmult UnfoldMult stuff 
 *
 * @ingroup pwglf_forward_multdist
 */
/** 
 * Get an object from a directory 
 * 
 * @param d     Directory 
 * @param name  Name of object 
 * 
 * @return Pointer to object or null if not found 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TObject* GetObject(TDirectory* d, const char* name)
{
  TObject* o = d->Get(name);
  if (!o) { 
    Error("GetCollection", "%s not found in %s", name, d->GetName());
    return 0;
  }
  return o;
}
/** 
 * Get an object from a collection
 * 
 * @param d     Collection 
 * @param name  Name of object
 * 
 * @return Pointer to object or null if not found 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TObject* GetObject(TCollection* d, const char* name)
{
  TObject* o = d->FindObject(name);
  if (!o) { 
    Error("GetCollection", "%s not found in %s", name, d->GetName());
    d->ls();
    return 0;
  }
  return o;
}
/** 
 * Get a collection from a directory
 * 
 * @param d     Parent directory 
 * @param name  Name of collection
 * 
 * @return Pointer to collection or null 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TCollection* GetCollection(TDirectory* d, const char* name)
{
  return static_cast<TCollection*>(GetObject(d, name));
}
/** 
 * Get a collection from a collection
 * 
 * @param d     Parent collection
 * @param name  Name of collection
 * 
 * @return Pointer to collection or null 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TCollection* GetCollection(TCollection* d, const char* name)
{
  return static_cast<TCollection*>(GetObject(d, name));
}
/** 
 * Get a 1D-histogram from a directory 
 * 
 * @param d     Parent directory 
 * @param name  Name of histogram 
 * 
 * @return Pointer to histogram or null 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TH1* GetH1(TDirectory* d, const char* name)
{
  return static_cast<TH1*>(GetObject(d, name));
}
/** 
 * Get a 1D-histogram from a collection 
 * 
 * @param d     Parent collectoin 
 * @param name  Name of histogram 
 * 
 * @return Pointer to histogram or null 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TH1* GetH1(TCollection* d, const char* name)
{
  return static_cast<TH1*>(GetObject(d, name));
}
/** 
 * Get a 2D-histogram from a directory 
 * 
 * @param d     Parent directory 
 * @param name  Name of histogram 
 * 
 * @return Pointer to histogram or null 
 */
TH2* GetH2(TDirectory* d, const char* name)
{
  return static_cast<TH2*>(GetObject(d, name));
}
/** 
 * Get a 2D-histogram from a collection
 * 
 * @param d     Parent collection
 * @param name  Name of histogram 
 * 
 * @return Pointer to histogram or null 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TH2* GetH2(TCollection* d, const char* name)
{
  return static_cast<TH2*>(GetObject(d, name));
}
/** 
 * Open a file 
 * 
 * @param filename      File name 
 * @param readNotWrite  If true, open read-only
 * 
 * @return Pointer to file, or null
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TFile* OpenFile(const char* filename, Bool_t readNotWrite)
{
  TFile* ret = TFile::Open(filename, readNotWrite ? "READ" : "RECREATE");
  if (!ret) { 
    Error("OpenFile", "Failed to open the file \"%s\"", filename);
    return 0;
  }
  return ret;
}
/** 
 * A type definition to make it easier to switch to TDirectory structure 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
typedef TCollection Dir;

/** 
 * Get histograms for @f$\eta@f$-bin from @a lowLim to @a highLim and
 * process it .
 * 
 * @param lowLim   Lower @f$\eta@f$ limit
 * @param highLim  Upper @f$\eta@f$ limit 
 * @param resp     Directory containing response matrices 
 * @param data     Directory containing data histogram 
 * @param out      Output directory  
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
void GetHistos(Double_t    lowLim, 
	       Double_t    highLim, 
	       Dir*        resp, 
	       Dir*        data, 
	       TDirectory* out);
/** 
 * Do the actual unfolding and corrections 
 * 
 * @param response   Response matrix 
 * @param measured   Measured distribution
 * @param mcHist     MC truth distribution
 * @param esdHist    MC distribution after ESD event selection
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist); 
/** 
 * Caclulate the trigger bias correction 
 * 
 * @param mcHist     MC truth distribution
 * @param esdHist    MC distribution after ESD event selection
 * 
 * @return  Histogram of trigger bias correction
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist);
/** 
 * Apply the trigger bias correction to the unfolded result
 * 
 * @param unfoldBefore Unfolded measured distribution
 * @param triggerBias  Trigger bias correction 
 * @param mcHist       MC truth distribution
 */
void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1* mcHist);


/** 
 * Main Loop to call eta range
 * 
 * @param respFile File containing response matrices 
 * @param dataFile File containing measured distributions
 * @param outFile  Output file 
 *
 * @ingroup pwglf_forward_scripts_unfoldmult
 */
void UnfoldMult(const char* respFile="forward_response.root", 
		const char* dataFile="forward_multiplicy.root", 
		const char* outFile="unfolded.root") 
{
  TFile* output = OpenFile(outFile,false);
  TFile* resp   = OpenFile(respFile, true);
  TFile* data   = OpenFile(dataFile, true);
  if (!output || !resp || !data) return;

  Dir* topResp  = GetCollection(resp, "ResponseMatricesSums");
  Dir* topData  = GetCollection(data, "MultSums");
  if (!topData || !topResp) return;
  
  Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. };
  Double_t* limit = limits;

  while ((*limit) > 0.1) {
    if ((*limit) >5.) GetHistos(-3.4,      (*limit), topResp, topData, output);
    else              GetHistos(-(*limit), (*limit), topResp, topData, output);
    limit++; 
    output->cd();
  }
  output->Close();
  resp->Close();
  data->Close();
}

//____________________________________________________________________
void GetHistos(Double_t    lowLim, 
	       Double_t    highLim, 
	       Dir*        topResp, 
	       Dir*        topData,
	       TDirectory* out) 
{
  // Format the name of the bin 
  TString sLow( TString::Format("%+03d",int(10*lowLim)));
  TString sHigh(TString::Format("%+03d",int(10*highLim)));
  sLow .ReplaceAll("+", "plus");
  sLow .ReplaceAll("-", "minus");
  sHigh.ReplaceAll("+", "plus");
  sHigh.ReplaceAll("-", "minus");
  TString     name = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
  TDirectory* dir = out->mkdir(name);
  dir->cd();

  // Get our collections 
  Dir* listResp = GetCollection(topResp, name);
  Dir* listData = GetCollection(topData,name);
  if(!listResp || !listData) return;

  // Get the MC histograms 
  TH2* response = GetH2(listResp, "responseMatrix");
  TH1* mcHist   = GetH1(listResp, "fMCNSD");
  TH1* esdHist  = GetH1(listResp, "fESDNSD");
  if (!response || !mcHist || !esdHist) return;

  // Get the data histogram 
  TH1* measured = GetH1(listData, "mult");
  if(!measured) return; 

  // Now do the unfolding 
  ProcessUnfold(response, measured, mcHist, esdHist);
} 


//____________________________________________________________________
void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist) 
{
  Int_t nX   = response->GetNbinsX();
  Int_t nY   = response->GetNbinsY();
  TH1* projY = static_cast<TH1*>(response->ProjectionY("projY",1,nX,""));
  TH1* projX = static_cast<TH1*>(response->ProjectionX("projX",1,nY,""));
  projX->SetDirectory(0);
  projY->SetDirectory(0);
  
  TH2* responseTranspose = static_cast<TH2*>(response->Clone("response"));
  for(Int_t i = 1; i <= nX; i++) {
    for(Int_t j = 1; j <= nY; j++) {
      responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j));
      responseTranspose->SetBinError(j,i, response->GetBinError(i,j));
    }
  }
  
  RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX, 
							    responseTranspose, 
							    "", "");
  RooUnfold*         unfold         = new RooUnfoldBayes(responseObject, 
							 measured, 5); 
  TH1*               unfolded       = static_cast<TH1*>(unfold->Hreco());
  TH1*               triggerBias    = GetTriggerBias(mcHist, esdHist);
  DoCorrection(unfolded, TriggerBias, mcHist);

  Double_t scale_unf = 1/unfolded->Integral();
  unfolded->Scale(scale_unf);
  unfolded->Write("Unfolded");

  Double_t scale_raw = 1/measured->Integral();
  measured->Scale(scale_raw);
  measured->Write("Raw");
  triggerBias->Write("TriggerBiasCorr");
}


//____________________________________________________________________
TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist) 
{
  mcHist->Sumw2();
  esdHist->Sumw2();
  
  TH1* corr = new TH1D("corr", "corr", 40, -0.5, 39.5);  
  for (Int_t j=1; j<corr->GetNbinsX(); j++) {
    Double_t errSqr = 0.;
    Double_t cESD   = esdHist->GetBinContent(j);
    Double_t cMS    = mcHist->GetBinContent(j);
    Double_t c      = (cMC  == 0 ? 1 : 
		       cESD == 0 ? 1/cMC : cESD/cMC);
    corr->SetBinContent(j, c);
    corr->SetBinError(j, c*TMath::Sqrt(errSqr));
  }
  return corr;  
}

//____________________________________________________________________
void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1*) 
{
  for (Int_t k=1; k<=35; k++) {
    Double_t tb = triggerBias->GetBinContent(k);
    Double_t ub = unfoldBefore->GetBinContent(k);
    if (tb > 1e-5 && ub > 0) {
      unfoldBefore->SetBinContent(k, ub / tb);
      Double_t eub = UnfoldBefore->GetBinError(k);
      Double_t etb = TriggerBias->GetBinError(k);
      unfoldBefore->SetBinError(k, TMath::Sqrt(TMath::Power(eub/ub)+
					       TMath::Power(etb/tb))*ub/tb); 
    }
  }
}
//
// EOF
//
Commit	Line	Data
4f09623b	1	/**
671df6c9	2	* @file UnfoldMult.C
4f09623b	3	* @author Valentina Zaccolo
	4	* @date Tue Mar 05 12:56:56 2013
	5	*
	6	* @brief Bayesian Method Unfolding Script
	7	*
671df6c9	8	* @ingroup pwglf_forward_scripts
671df6c9	9	* @ingroup pwglf_forward_multdist
4f09623b	10	*/
	11	#include <iostream>
	12	#include <TList.h>
	13	#include <TFile.h>
	14	#include <TObject.h>
	15	#include <TDirectory.h>
	16	#include "TH1D.h"
	17	#include "TH2D.h"
	18	#include <TROOT.h>
	19	#include <TSystem.h>
	20	#include "RooUnfoldBayes.h"
	21	#include "RooUnfoldResponse.h"
	22	#include <TMath.h>
671df6c9	23	/**
	24	* @defgroup pwglf_forward_scripts_unfoldmult UnfoldMult stuff
	25	*
	26	* @ingroup pwglf_forward_multdist
	27	*/
4f09623b	28	/**
	29	* Get an object from a directory
	30	*
	31	* @param d Directory
	32	* @param name Name of object
	33	*
	34	* @return Pointer to object or null if not found
671df6c9	35	*
671df6c9	36	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	37	*/
	38	TObject* GetObject(TDirectory* d, const char* name)
	39	{
	40	TObject* o = d->Get(name);
	41	if (!o) {
	42	Error("GetCollection", "%s not found in %s", name, d->GetName());
	43	return 0;
	44	}
	45	return o;
	46	}
	47	/**
	48	* Get an object from a collection
	49	*
	50	* @param d Collection
	51	* @param name Name of object
	52	*
	53	* @return Pointer to object or null if not found
671df6c9	54	*
671df6c9	55	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	56	*/
	57	TObject* GetObject(TCollection* d, const char* name)
	58	{
	59	TObject* o = d->FindObject(name);
	60	if (!o) {
	61	Error("GetCollection", "%s not found in %s", name, d->GetName());
	62	d->ls();
	63	return 0;
	64	}
	65	return o;
	66	}
	67	/**
	68	* Get a collection from a directory
	69	*
	70	* @param d Parent directory
	71	* @param name Name of collection
	72	*
	73	* @return Pointer to collection or null
671df6c9	74	*
671df6c9	75	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	76	*/
	77	TCollection* GetCollection(TDirectory* d, const char* name)
	78	{
	79	return static_cast<TCollection*>(GetObject(d, name));
	80	}
	81	/**
	82	* Get a collection from a collection
	83	*
	84	* @param d Parent collection
	85	* @param name Name of collection
	86	*
	87	* @return Pointer to collection or null
671df6c9	88	*
671df6c9	89	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	90	*/
	91	TCollection* GetCollection(TCollection* d, const char* name)
	92	{
	93	return static_cast<TCollection*>(GetObject(d, name));
	94	}
	95	/**
	96	* Get a 1D-histogram from a directory
	97	*
	98	* @param d Parent directory
	99	* @param name Name of histogram
	100	*
	101	* @return Pointer to histogram or null
671df6c9	102	*
671df6c9	103	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	104	*/
	105	TH1* GetH1(TDirectory* d, const char* name)
	106	{
	107	return static_cast<TH1*>(GetObject(d, name));
	108	}
	109	/**
	110	* Get a 1D-histogram from a collection
	111	*
	112	* @param d Parent collectoin
	113	* @param name Name of histogram
	114	*
	115	* @return Pointer to histogram or null
671df6c9	116	*
671df6c9	117	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	118	*/
	119	TH1* GetH1(TCollection* d, const char* name)
	120	{
	121	return static_cast<TH1*>(GetObject(d, name));
	122	}
	123	/**
	124	* Get a 2D-histogram from a directory
	125	*
	126	* @param d Parent directory
	127	* @param name Name of histogram
	128	*
	129	* @return Pointer to histogram or null
	130	*/
	131	TH2* GetH2(TDirectory* d, const char* name)
	132	{
	133	return static_cast<TH2*>(GetObject(d, name));
	134	}
	135	/**
	136	* Get a 2D-histogram from a collection
	137	*
	138	* @param d Parent collection
	139	* @param name Name of histogram
	140	*
	141	* @return Pointer to histogram or null
671df6c9	142	*
671df6c9	143	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	144	*/
	145	TH2* GetH2(TCollection* d, const char* name)
	146	{
	147	return static_cast<TH2*>(GetObject(d, name));
	148	}
	149	/**
	150	* Open a file
	151	*
	152	* @param filename File name
	153	* @param readNotWrite If true, open read-only
	154	*
	155	* @return Pointer to file, or null
671df6c9	156	*
671df6c9	157	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	158	*/
	159	TFile* OpenFile(const char* filename, Bool_t readNotWrite)
	160	{
	161	TFile* ret = TFile::Open(filename, readNotWrite ? "READ" : "RECREATE");
	162	if (!ret) {
	163	Error("OpenFile", "Failed to open the file \"%s\"", filename);
	164	return 0;
	165	}
	166	return ret;
	167	}
671df6c9	168	/**
	169	* A type definition to make it easier to switch to TDirectory structure
	170	*
	171	* @ingroup pwglf_forward_scripts_unfoldmult
	172	*/
4f09623b	173	typedef TCollection Dir;
	174
	175	/**
	176	* Get histograms for @f$\eta@f$-bin from @a lowLim to @a highLim and
	177	* process it .
	178	*
	179	* @param lowLim Lower @f$\eta@f$ limit
	180	* @param highLim Upper @f$\eta@f$ limit
	181	* @param resp Directory containing response matrices
	182	* @param data Directory containing data histogram
	183	* @param out Output directory
671df6c9	184	*
671df6c9	185	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	186	*/
	187	void GetHistos(Double_t lowLim,
	188	Double_t highLim,
	189	Dir* resp,
	190	Dir* data,
	191	TDirectory* out);
	192	/**
	193	* Do the actual unfolding and corrections
	194	*
	195	* @param response Response matrix
	196	* @param measured Measured distribution
	197	* @param mcHist MC truth distribution
	198	* @param esdHist MC distribution after ESD event selection
671df6c9	199	*
671df6c9	200	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	201	*/
671df6c9	202	void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist);
4f09623b	203	/**
	204	* Caclulate the trigger bias correction
	205	*
	206	* @param mcHist MC truth distribution
	207	* @param esdHist MC distribution after ESD event selection
	208	*
	209	* @return Histogram of trigger bias correction
671df6c9	210	*
671df6c9	211	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	212	*/
	213	TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist);
	214	/**
	215	* Apply the trigger bias correction to the unfolded result
	216	*
	217	* @param unfoldBefore Unfolded measured distribution
	218	* @param triggerBias Trigger bias correction
	219	* @param mcHist MC truth distribution
	220	*/
	221	void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1* mcHist);
	222
	223
	224	/**
	225	* Main Loop to call eta range
	226	*
	227	* @param respFile File containing response matrices
	228	* @param dataFile File containing measured distributions
	229	* @param outFile Output file
671df6c9	230	*
671df6c9	231	* @ingroup pwglf_forward_scripts_unfoldmult
4f09623b	232	*/
	233	void UnfoldMult(const char* respFile="forward_response.root",
	234	const char* dataFile="forward_multiplicy.root",
	235	const char* outFile="unfolded.root")
	236	{
	237	TFile* output = OpenFile(outFile,false);
	238	TFile* resp = OpenFile(respFile, true);
	239	TFile* data = OpenFile(dataFile, true);
	240	if (!output \|\| !resp \|\| !data) return;
	241
	242	Dir* topResp = GetCollection(resp, "ResponseMatricesSums");
	243	Dir* topData = GetCollection(data, "MultSums");
	244	if (!topData \|\| !topResp) return;
	245
	246	Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. };
	247	Double_t* limit = limits;
	248
	249	while ((*limit) > 0.1) {
	250	if ((limit) >5.) GetHistos(-3.4, (limit), topResp, topData, output);
	251	else GetHistos(-(limit), (limit), topResp, topData, output);
	252	limit++;
	253	output->cd();
	254	}
	255	output->Close();
	256	resp->Close();
	257	data->Close();
	258	}
	259
	260	//____________________________________________________________________
	261	void GetHistos(Double_t lowLim,
	262	Double_t highLim,
	263	Dir* topResp,
	264	Dir* topData,
	265	TDirectory* out)
	266	{
	267	// Format the name of the bin
	268	TString sLow( TString::Format("%+03d",int(10*lowLim)));
	269	TString sHigh(TString::Format("%+03d",int(10*highLim)));
	270	sLow .ReplaceAll("+", "plus");
	271	sLow .ReplaceAll("-", "minus");
	272	sHigh.ReplaceAll("+", "plus");
	273	sHigh.ReplaceAll("-", "minus");
	274	TString name = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
	275	TDirectory* dir = out->mkdir(name);
	276	dir->cd();
	277
	278	// Get our collections
	279	Dir* listResp = GetCollection(topResp, name);
	280	Dir* listData = GetCollection(topData,name);
	281	if(!listResp \|\| !listData) return;
	282
	283	// Get the MC histograms
	284	TH2* response = GetH2(listResp, "responseMatrix");
	285	TH1* mcHist = GetH1(listResp, "fMCNSD");
	286	TH1* esdHist = GetH1(listResp, "fESDNSD");
	287	if (!response \|\| !mcHist \|\| !esdHist) return;
	288
	289	// Get the data histogram
	290	TH1* measured = GetH1(listData, "mult");
	291	if(!measured) return;
	292
	293	// Now do the unfolding
	294	ProcessUnfold(response, measured, mcHist, esdHist);
	295	}
296
297
298
299	//____________________________________________________________________
300	void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist)
301	{
302	Int_t nX = response->GetNbinsX();
303	Int_t nY = response->GetNbinsY();
304	TH1* projY = static_cast<TH1*>(response->ProjectionY("projY",1,nX,""));
305	TH1* projX = static_cast<TH1*>(response->ProjectionX("projX",1,nY,""));
306	projX->SetDirectory(0);
307	projY->SetDirectory(0);
308
309	TH2* responseTranspose = static_cast<TH2*>(response->Clone("response"));
310	for(Int_t i = 1; i <= nX; i++) {
311	for(Int_t j = 1; j <= nY; j++) {
312	responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j));
313	responseTranspose->SetBinError(j,i, response->GetBinError(i,j));
314	}
315	}
316
317	RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX,
318	responseTranspose,
319	"", "");
320	RooUnfold* unfold = new RooUnfoldBayes(responseObject,
321	measured, 5);
322	TH1* unfolded = static_cast<TH1*>(unfold->Hreco());
323	TH1* triggerBias = GetTriggerBias(mcHist, esdHist);
324	DoCorrection(unfolded, TriggerBias, mcHist);
325
326	Double_t scale_unf = 1/unfolded->Integral();
327	unfolded->Scale(scale_unf);
328	unfolded->Write("Unfolded");
329
330	Double_t scale_raw = 1/measured->Integral();
331	measured->Scale(scale_raw);
332	measured->Write("Raw");
333	triggerBias->Write("TriggerBiasCorr");
334	}
335
336
337
338	//____________________________________________________________________
339	TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist)
340	{
341	mcHist->Sumw2();
342	esdHist->Sumw2();
343
344	TH1* corr = new TH1D("corr", "corr", 40, -0.5, 39.5);
345	for (Int_t j=1; j<corr->GetNbinsX(); j++) {
346	Double_t errSqr = 0.;
347	Double_t cESD = esdHist->GetBinContent(j);
348	Double_t cMS = mcHist->GetBinContent(j);
349	Double_t c = (cMC == 0 ? 1 :
350	cESD == 0 ? 1/cMC : cESD/cMC);
351	corr->SetBinContent(j, c);
352	corr->SetBinError(j, c*TMath::Sqrt(errSqr));
353	}
354	return corr;
355	}
356
357	//____________________________________________________________________
358	void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1*)
359	{
360	for (Int_t k=1; k<=35; k++) {
361	Double_t tb = triggerBias->GetBinContent(k);
362	Double_t ub = unfoldBefore->GetBinContent(k);
363	if (tb > 1e-5 && ub > 0) {
364	unfoldBefore->SetBinContent(k, ub / tb);
365	Double_t eub = UnfoldBefore->GetBinError(k);
366	Double_t etb = TriggerBias->GetBinError(k);
367	unfoldBefore->SetBinError(k, TMath::Sqrt(TMath::Power(eub/ub)+
368	TMath::Power(etb/tb))*ub/tb);
369	}
370	}
371	}
372	//
373	// EOF
374	//
375
376
377