[u/mrichter/AliRoot.git] / PWGLF / FORWARD / corrections / Empirical / GenerateEmpirical.C

/** 
 * Open a file 
 * 
 * @param filename Name of file 
 * 
 * @return Return value 
 */
TFile* OpenFile(const char* filename)
{
  TFile* file = TFile::Open(filename, "READ");
  if (!file) 
    Error("OpenFile", "Failed to open the file %s for reading", filename);
  return file;
}
/** 
 * Get an object from parent list 
 * 
 * @param l    Parent list 
 * @param name Name of object to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TObject* GetObject(const TList* l, const char* name)
{
  if (!l) {
    Warning("GetObject", "No parent list given");
    return 0;
  }
  TObject* o = l->FindObject(name);
  if (!o) { 
    Warning("GetObject", "Object named %s not found in list %s", 
	    name, l->GetName());
    return 0;
  }
  return o;
}
/** 
 * Get an object from parent list 
 * 
 * @param l    Parent list 
 * @param name Name of object to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TObject* GetObject(const TDirectory* l, const char* name)
{
  if (!l) {
    Warning("GetObject", "No parent directory given");
    return 0;
  }
  TObject* o = l->Get(name);
  if (!o) { 
    Warning("GetObject", "Object named %s not found in directory %s", 
	    name, l->GetName());
    return 0;
  }
  return o;
}
Bool_t CheckType(const TObject* o, const TClass* cl)
{
  if (!o) return false;
  if (!o->IsA()->InheritsFrom(cl)) { 
    Warning("CheckType", "Object %s is not a %s, but a %s", 
	    o->GetName(), cl->GetName(), o->ClassName());
    return false;
  }
  return true;
}
/** 
 * Get a list from another list 
 * 
 * @param l    Parent list 
 * @param name Name of list to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TList* GetList(const TList* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TList::Class())) return 0;
  return static_cast<TList*>(o);
}
/** 
 * Get a list from another list 
 * 
 * @param l    Parent list 
 * @param name Name of list to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TList* GetList(const TDirectory* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TList::Class())) return 0;
  return static_cast<TList*>(o);
}
/** 
 * Get a histogram from a list 
 * 
 * @param l    Parent list 
 * @param name Name of histogram to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TH1* GetHist(const TList* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TH1::Class())) return 0;
  return static_cast<TH1*>(o);
} 	  
/** 
 * Get a histogram from a list 
 * 
 * @param l    Parent list 
 * @param name Name of histogram to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TH1* GetHist(const TDirectory* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TH1::Class())) return 0;
  return static_cast<TH1*>(o);
} 	  
/** 
 * Get a histogram from a list 
 * 
 * @param l    Parent list 
 * @param name Name of histogram to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TGraphErrors* GetGraph(const TList* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TGraphErrors::Class())) return 0;
  return static_cast<TGraphErrors*>(o);
} 	  
/** 
 * Get a histogram from a list 
 * 
 * @param l    Parent list 
 * @param name Name of histogram to find 
 * 
 * @return 0 in case of problems, object pointer otherwise 
 */
TGraphErrors* GetGraph(const TDirectory* l, const char* name)
{
  TObject* o = GetObject(l, name);
  if (!CheckType(o, TGraphErrors::Class())) return 0;
  return static_cast<TGraphErrors*>(o);
} 	  

TGraphErrors* Ratio(const TGraphErrors* num, const TGraphErrors* denom)
{
  TGraphErrors* r = static_cast<TGraphErrors*>(num->Clone("tmp"));
  for (Int_t k = 0; k < r->GetN(); k++) { 
    Double_t xn = r->GetX()[k];
    Double_t yn = r->GetY()[k];
    Double_t en = r->GetErrorY(k);
    // Double_t xd = denom->GetX()[k];
    Double_t yd = denom->GetY()[k];
    Double_t ed = denom->GetErrorY(k);
    Double_t s  = (yn-yd)/yd;
    Double_t e  = 0;
    r->SetPoint(k, xn, s);
    r->SetPointError(k, 0, e);
  }
  return r;
}

/** 
 * Generate the emperical correction for unknown material in the
 * AliROOT geometric description.  This is done by taking the ratio
 * of the nominal vertex results to the results from satelitte events. 
 *
 * This can be done in two ways 
 *
 * - Either by using FMD results only 
 * - Or by using the full combined FMD+V0+SPD results 
 * 
 * The correction is written to the file "EmpiricalCorrection.root",
 * which will contain a TGraphErrors for each defined centrality
 * class, plus a TGraphErrors object that contains the average over
 * all centrality classes.
 *
 * Other results should be DIVIDED by the correction obtained from
 * this script.
 * 
 * @param fmdfmd If true, use FMD results only 
 */
void 
GenerateEmpirical(const char* fmdDispVtx="forward_dndetaDispVtx.root",
		  const char* fmdNomVtx="forward_dndetaNominalVtxZDC.root",
		  const char* fullDispVtx="dndetaPreliminaryQM12.root")
{
  
  TString outName    = "EmpiricalCorrection.root";
  TFile* fdispvtxFMD = OpenFile(fmdDispVtx);
  TFile* fNominalFMD = OpenFile(fmdNomVtx);
  TFile* fdispvtxAll = OpenFile(fullDispVtx);
  if (!fdispvtxFMD || !fNominalFMD || !fdispvtxAll) return;

  TFile* out = TFile::Open(outName, "RECREATE");
  if (!out) { 
    Error("GenerateEmpirical", "Failed to open output file %s", 
	  outName.Data());
    return;
  }

  Int_t limits[] = {     0,        5,      10,         20,    30  };
  Int_t colors[] = {kRed+2, kGreen+2, kBlue+1, kMagenta+1, kBlack };

  TGraphErrors* afmdfmd  = 0;
  TGraphErrors* afmdfull = 0;

  // --- Do two things fmd-to-fmd and fmd-to-full
  for (Int_t iM = 0; iM < 2; iM++) {
    Bool_t       fmdfmd   = (iM == 0);
    TMultiGraph* mg       = new TMultiGraph();
    TMultiGraph* mgfmdnom = (iM == 0 ? new TMultiGraph() : 0);
    TMultiGraph* mgref    = new TMultiGraph();
    mgref->SetTitle(Form("Satellite, FMD%s", fmdfmd ? "" : "+V0+Tracklets"));
    if (mgfmdnom) mgfmdnom->SetTitle("Nominal FMD");

    Info("", "Doing %s/FMD", fmdfmd ? "FMD" : "Full");
    TGraphErrors* corrcent[4];    
    TGraphErrors* allsym     = 0;
    TGraphErrors* fmddispvtx = 0;
    TGraphErrors* fmdnominal = 0;

    TList* displist = GetList(fdispvtxFMD,"ForwardResults");
    TList* nomlist  = GetList(fNominalFMD,"ForwardResults");
    
    Int_t nMin = 1000;
    // --- Get each defined centrality and form the ratio ------------
    for(Int_t iC=0; iC<4; iC++) {
      Int_t cl = limits[iC];
      Int_t ch = limits[iC+1];
      TString base; base.Form("_cent_%03d_%03d", cl, ch);
      corrcent[iC] = new TGraphErrors;
      corrcent[iC]->SetName(Form("correction%s",base.Data()));
      corrcent[iC]->SetTitle(Form("%2d%% - %2d%%", cl, ch));
      corrcent[iC]->GetHistogram()->SetXTitle("#eta");
      corrcent[iC]->GetHistogram()->SetYTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
					     "dN_{ch}/d#eta|_{satelitte}}");  
      corrcent[iC]->SetMarkerColor(colors[iC]);
      corrcent[iC]->SetLineColor(colors[iC]);
      corrcent[iC]->SetMarkerStyle(fmdfmd ? 20 : 24);
      corrcent[iC]->SetMarkerSize(fmdfmd ? 1 : 1.2);
      corrcent[iC]->SetFillStyle(0);
      corrcent[iC]->SetFillColor(0);
      
      mg->Add(corrcent[iC]);
      TGraphErrors* allsym = GetGraph(fdispvtxAll,Form("graphErrors_cent_%d_%d",
						       cl, ch));
    
      TString folderName   = Form("cent%03d_%03d",cl, ch);
      TList*  dispcentlist = GetList(displist, folderName);
      TList*  nomcentlist  = GetList(nomlist, folderName);
      TString histName     = "dndetaForward_rebin05";
      TH1D*   hDisp        = GetHist(dispcentlist, histName);
      TH1D*   hNominal     = GetHist(nomcentlist, histName);
    
      // Make our temporary graph 
      if   (fmdfmd) fmddispvtx = new TGraphErrors(hDisp);
      else          fmddispvtx = static_cast<TGraphErrors*>(allsym->Clone());
      fmdnominal               = new TGraphErrors(hNominal);
      if (mgfmdnom) {
	TGraph* g = static_cast<TGraph*>(fmdnominal->Clone(Form("nominal%s", 
								base.Data())));
	g->SetMarkerColor(colors[iC]);
	g->SetLineColor(colors[iC]);
	g->SetMarkerStyle(21);
	g->SetTitle("Nominal FMD");
	mgfmdnom->Add(g);
      }
      TGraph* ref = static_cast<TGraph*>(fmddispvtx->Clone(Form("satellite%s",
								base.Data())));
      ref->SetTitle(Form("%2d - %2d", cl, ch));
      ref->SetMarkerColor(colors[iC]);
      ref->SetLineColor(colors[iC]);
      ref->SetMarkerStyle(fmdfmd ? 22 : 20);
      ref->SetMarkerSize(1);
      mgref->Add(ref);
    
      Int_t nPoints = 0;
    
      for(Int_t iN=0; iN<fmdnominal->GetN(); iN++) {
	Double_t eta        = fmdnominal->GetX()[iN];
	Double_t nommult    = fmdnominal->GetY()[iN];
	Double_t nommulterr = fmdnominal->GetErrorY(iN);
	Double_t etaerr     = fmdnominal->GetErrorX(iN);

	// Ignore empty bins 
	if(nommult < 0.0001) continue;

	// Find the corresponding bin from the dispaclaced vertex analysis
	for(Int_t iS=0; iS < fmddispvtx->GetN(); iS++) {
	  Double_t eta1        = fmddispvtx->GetX()[iS];
	  Double_t dispmult    = fmddispvtx->GetY()[iS];
	  Double_t dispmulterr = fmddispvtx->GetErrorY(iS);

	  if(TMath::Abs(eta-eta1) >= 0.001) continue;
	  
	  Double_t corr  = nommult/dispmult;
	  Double_t rd    = dispmulterr/dispmult;
	  Double_t rn    = nommulterr/nommult;
	  Double_t error = (1/corr)*TMath::Sqrt(TMath::Power(rd,2) + 
					      TMath::Power(rn,2));
	  corrcent[iC]->SetPoint(nPoints,eta,corr);
	  corrcent[iC]->SetPointError(nPoints,etaerr,error);
	  nPoints++;
	  
	}
	//std::cout<<eta<<"  "<<nommult<<std::endl;
      }
      nMin = TMath::Min(nPoints, nMin);
    }

    // --- Calulate the average --------------------------------------
    TGraphErrors* average = new TGraphErrors;
    average->SetName(Form("average"));
    average->SetTitle(Form("%2d%% - %2d%%", limits[0], limits[4]));
    average->SetMarkerColor(colors[4]);
    average->SetLineColor(colors[4]);
    average->SetMarkerStyle(fmdfmd ? 20 : 24);
    average->SetMarkerSize(fmdfmd ? 1 : 1.2);
    average->SetFillStyle(0);
    average->SetFillColor(0);
    mg->Add(average);
    
    for(Int_t iA = 0; iA < nMin; iA++) {
      Double_t mean   = 0;
      Double_t error2 = 0;
      Double_t sumw2  = 0;
      
      // Loop over centralities 
      for(Int_t iC=0; iC<4; iC++) {
	Double_t eta  = corrcent[iC]->GetX()[iA];
	Double_t corr = corrcent[iC]->GetY()[iA];
	Double_t err  = corrcent[iC]->GetErrorY(iA);
	Double_t err2 = err * err;
	Double_t w    = 1 / err2;
	sumw2         += w;
	mean          += w * corr;
      }
      mean           /= sumw2;
      Double_t error = TMath::Sqrt(1. / sumw2);
          average->SetPoint(iA,eta,mean);
      average->SetPointError(iA,0.125,error);    
    }
    if (fmdfmd) afmdfmd  = average;
    else        afmdfull = average;

    // --- Calculate ratios ------------------------------------------
    TMultiGraph* ratios = new TMultiGraph;
    ratios->SetName("ratios");
    for(Int_t iC=0; iC<4; iC++) {
      Int_t cl        = limits[iC];
      Int_t ch        = limits[iC+1];
      TGraphErrors* r = Ratio(corrcent[iC], average);
      r->SetName(Form("ratio%s",base.Data()));
      ratios->Add(r);
    }

    // --- Store the result ------------------------------------------
    TDirectory* d = out->mkdir(fmdfmd ? "fmdfmd" : "fmdfull", 
			       Form("Empirical corrections %s", 
				    fmdfmd ? "FMD/FMD" : "FULL/FMD"));
    d->cd();
    for(Int_t p=0; p<4; p++) corrcent[p]->Write();
    average->Write("average");
    mg->Write("all");
    ratios->Write("ratios");
    mgref->Write("satellite");
    out->cd();
    if (mgfmdnom) mgfmdnom->Write("nominal");
    out->ls();
  } // End method loop

  TGraphErrors* r = Ratio(afmdfull, afmdfmd);
  r->Write("ratio");

  out->Close();
}
//
// EOF
//
Commit	Line	Data
a1593a3b	1	/**
	2	* Open a file
	3	*
	4	* @param filename Name of file
	5	*
	6	* @return Return value
	7	*/
	8	TFile* OpenFile(const char* filename)
	9	{
	10	TFile* file = TFile::Open(filename, "READ");
	11	if (!file)
	12	Error("OpenFile", "Failed to open the file %s for reading", filename);
	13	return file;
	14	}
	15	/**
	16	* Get an object from parent list
	17	*
	18	* @param l Parent list
	19	* @param name Name of object to find
	20	*
	21	* @return 0 in case of problems, object pointer otherwise
	22	*/
	23	TObject* GetObject(const TList* l, const char* name)
	24	{
	25	if (!l) {
	26	Warning("GetObject", "No parent list given");
	27	return 0;
	28	}
	29	TObject* o = l->FindObject(name);
	30	if (!o) {
	31	Warning("GetObject", "Object named %s not found in list %s",
	32	name, l->GetName());
	33	return 0;
	34	}
	35	return o;
	36	}
	37	/**
	38	* Get an object from parent list
	39	*
	40	* @param l Parent list
	41	* @param name Name of object to find
	42	*
	43	* @return 0 in case of problems, object pointer otherwise
	44	*/
	45	TObject* GetObject(const TDirectory* l, const char* name)
	46	{
	47	if (!l) {
	48	Warning("GetObject", "No parent directory given");
	49	return 0;
	50	}
	51	TObject* o = l->Get(name);
	52	if (!o) {
	53	Warning("GetObject", "Object named %s not found in directory %s",
	54	name, l->GetName());
	55	return 0;
	56	}
	57	return o;
	58	}
	59	Bool_t CheckType(const TObject* o, const TClass* cl)
	60	{
	61	if (!o) return false;
	62	if (!o->IsA()->InheritsFrom(cl)) {
	63	Warning("CheckType", "Object %s is not a %s, but a %s",
	64	o->GetName(), cl->GetName(), o->ClassName());
65	return false;
66	}
67	return true;
68	}
69	/**
70	* Get a list from another list
71	*
72	* @param l Parent list
73	* @param name Name of list to find
74	*
75	* @return 0 in case of problems, object pointer otherwise
76	*/
77	TList* GetList(const TList* l, const char* name)
78	{
79	TObject* o = GetObject(l, name);
80	if (!CheckType(o, TList::Class())) return 0;
81	return static_cast<TList*>(o);
82	}
83	/**
84	* Get a list from another list
85	*
86	* @param l Parent list
87	* @param name Name of list to find
88	*
89	* @return 0 in case of problems, object pointer otherwise
90	*/
91	TList* GetList(const TDirectory* l, const char* name)
92	{
93	TObject* o = GetObject(l, name);
94	if (!CheckType(o, TList::Class())) return 0;
95	return static_cast<TList*>(o);
96	}
97	/**
98	* Get a histogram from a list
99	*
100	* @param l Parent list
101	* @param name Name of histogram to find
102	*
103	* @return 0 in case of problems, object pointer otherwise
104	*/
105	TH1* GetHist(const TList* l, const char* name)
106	{
107	TObject* o = GetObject(l, name);
108	if (!CheckType(o, TH1::Class())) return 0;
109	return static_cast<TH1*>(o);
110	}
111	/**
112	* Get a histogram from a list
113	*
114	* @param l Parent list
115	* @param name Name of histogram to find
116	*
117	* @return 0 in case of problems, object pointer otherwise
118	*/
119	TH1* GetHist(const TDirectory* l, const char* name)
120	{
121	TObject* o = GetObject(l, name);
122	if (!CheckType(o, TH1::Class())) return 0;
123	return static_cast<TH1*>(o);
124	}
125	/**
126	* Get a histogram from a list
127	*
128	* @param l Parent list
129	* @param name Name of histogram to find
130	*
131	* @return 0 in case of problems, object pointer otherwise
132	*/
133	TGraphErrors* GetGraph(const TList* l, const char* name)
134	{
135	TObject* o = GetObject(l, name);
136	if (!CheckType(o, TGraphErrors::Class())) return 0;
137	return static_cast<TGraphErrors*>(o);
138	}
139	/**
140	* Get a histogram from a list
141	*
142	* @param l Parent list
143	* @param name Name of histogram to find
144	*
145	* @return 0 in case of problems, object pointer otherwise
146	*/
147	TGraphErrors* GetGraph(const TDirectory* l, const char* name)
148	{
149	TObject* o = GetObject(l, name);
150	if (!CheckType(o, TGraphErrors::Class())) return 0;
151	return static_cast<TGraphErrors*>(o);
152	}
153
154	TGraphErrors* Ratio(const TGraphErrors* num, const TGraphErrors* denom)
155	{
156	TGraphErrors* r = static_cast<TGraphErrors*>(num->Clone("tmp"));
157	for (Int_t k = 0; k < r->GetN(); k++) {
158	Double_t xn = r->GetX()[k];
159	Double_t yn = r->GetY()[k];
160	Double_t en = r->GetErrorY(k);
161	// Double_t xd = denom->GetX()[k];
162	Double_t yd = denom->GetY()[k];
163	Double_t ed = denom->GetErrorY(k);
164	Double_t s = (yn-yd)/yd;
165	Double_t e = 0;
166	r->SetPoint(k, xn, s);
167	r->SetPointError(k, 0, e);
168	}
169	return r;
170	}
171
172	/**
173	* Generate the emperical correction for unknown material in the
174	* AliROOT geometric description. This is done by taking the ratio
175	* of the nominal vertex results to the results from satelitte events.
176	*
177	* This can be done in two ways
178	*
179	* - Either by using FMD results only
180	* - Or by using the full combined FMD+V0+SPD results
181	*
182	* The correction is written to the file "EmpiricalCorrection.root",
183	* which will contain a TGraphErrors for each defined centrality
184	* class, plus a TGraphErrors object that contains the average over
185	* all centrality classes.
186	*
187	* Other results should be DIVIDED by the correction obtained from
188	* this script.
189	*
190	* @param fmdfmd If true, use FMD results only
191	*/
192	void
193	GenerateEmpirical(const char* fmdDispVtx="forward_dndetaDispVtx.root",
194	const char* fmdNomVtx="forward_dndetaNominalVtxZDC.root",
195	const char* fullDispVtx="dndetaPreliminaryQM12.root")
196	{
197
198	TString outName = "EmpiricalCorrection.root";
199	TFile* fdispvtxFMD = OpenFile(fmdDispVtx);
200	TFile* fNominalFMD = OpenFile(fmdNomVtx);
201	TFile* fdispvtxAll = OpenFile(fullDispVtx);
202	if (!fdispvtxFMD \|\| !fNominalFMD \|\| !fdispvtxAll) return;
203
204	TFile* out = TFile::Open(outName, "RECREATE");
205	if (!out) {
206	Error("GenerateEmpirical", "Failed to open output file %s",
207	outName.Data());
208	return;
209	}
210
211	Int_t limits[] = { 0, 5, 10, 20, 30 };
212	Int_t colors[] = {kRed+2, kGreen+2, kBlue+1, kMagenta+1, kBlack };
213
214	TGraphErrors* afmdfmd = 0;
215	TGraphErrors* afmdfull = 0;
216
217	// --- Do two things fmd-to-fmd and fmd-to-full
218	for (Int_t iM = 0; iM < 2; iM++) {
219	Bool_t fmdfmd = (iM == 0);
220	TMultiGraph* mg = new TMultiGraph();
221	TMultiGraph* mgfmdnom = (iM == 0 ? new TMultiGraph() : 0);
222	TMultiGraph* mgref = new TMultiGraph();
223	mgref->SetTitle(Form("Satellite, FMD%s", fmdfmd ? "" : "+V0+Tracklets"));
224	if (mgfmdnom) mgfmdnom->SetTitle("Nominal FMD");
225
226	Info("", "Doing %s/FMD", fmdfmd ? "FMD" : "Full");
227	TGraphErrors* corrcent[4];
228	TGraphErrors* allsym = 0;
229	TGraphErrors* fmddispvtx = 0;
230	TGraphErrors* fmdnominal = 0;
231
232	TList* displist = GetList(fdispvtxFMD,"ForwardResults");
233	TList* nomlist = GetList(fNominalFMD,"ForwardResults");
234
235	Int_t nMin = 1000;
236	// --- Get each defined centrality and form the ratio ------------
237	for(Int_t iC=0; iC<4; iC++) {
238	Int_t cl = limits[iC];
239	Int_t ch = limits[iC+1];
240	TString base; base.Form("_cent_%03d_%03d", cl, ch);
241	corrcent[iC] = new TGraphErrors;
242	corrcent[iC]->SetName(Form("correction%s",base.Data()));
243	corrcent[iC]->SetTitle(Form("%2d%% - %2d%%", cl, ch));
244	corrcent[iC]->GetHistogram()->SetXTitle("#eta");
245	corrcent[iC]->GetHistogram()->SetYTitle("#frac{dN_{ch}/d#eta\|_{nominal}}{"
246	"dN_{ch}/d#eta\|_{satelitte}}");
247	corrcent[iC]->SetMarkerColor(colors[iC]);
248	corrcent[iC]->SetLineColor(colors[iC]);
249	corrcent[iC]->SetMarkerStyle(fmdfmd ? 20 : 24);
250	corrcent[iC]->SetMarkerSize(fmdfmd ? 1 : 1.2);
251	corrcent[iC]->SetFillStyle(0);
252	corrcent[iC]->SetFillColor(0);
253
254	mg->Add(corrcent[iC]);
255	TGraphErrors* allsym = GetGraph(fdispvtxAll,Form("graphErrors_cent_%d_%d",
256	cl, ch));
257
258	TString folderName = Form("cent%03d_%03d",cl, ch);
259	TList* dispcentlist = GetList(displist, folderName);
260	TList* nomcentlist = GetList(nomlist, folderName);
261	TString histName = "dndetaForward_rebin05";
262	TH1D* hDisp = GetHist(dispcentlist, histName);
263	TH1D* hNominal = GetHist(nomcentlist, histName);
264
265	// Make our temporary graph
266	if (fmdfmd) fmddispvtx = new TGraphErrors(hDisp);
267	else fmddispvtx = static_cast<TGraphErrors*>(allsym->Clone());
268	fmdnominal = new TGraphErrors(hNominal);
269	if (mgfmdnom) {
270	TGraph* g = static_cast<TGraph*>(fmdnominal->Clone(Form("nominal%s",
271	base.Data())));
272	g->SetMarkerColor(colors[iC]);
273	g->SetLineColor(colors[iC]);
274	g->SetMarkerStyle(21);
275	g->SetTitle("Nominal FMD");
276	mgfmdnom->Add(g);
277	}
278	TGraph* ref = static_cast<TGraph*>(fmddispvtx->Clone(Form("satellite%s",
279	base.Data())));
280	ref->SetTitle(Form("%2d - %2d", cl, ch));
281	ref->SetMarkerColor(colors[iC]);
282	ref->SetLineColor(colors[iC]);
283	ref->SetMarkerStyle(fmdfmd ? 22 : 20);
284	ref->SetMarkerSize(1);
285	mgref->Add(ref);
286
287	Int_t nPoints = 0;
288
289	for(Int_t iN=0; iN<fmdnominal->GetN(); iN++) {
290	Double_t eta = fmdnominal->GetX()[iN];
291	Double_t nommult = fmdnominal->GetY()[iN];
292	Double_t nommulterr = fmdnominal->GetErrorY(iN);
293	Double_t etaerr = fmdnominal->GetErrorX(iN);
294
295	// Ignore empty bins
296	if(nommult < 0.0001) continue;
297
298	// Find the corresponding bin from the dispaclaced vertex analysis
299	for(Int_t iS=0; iS < fmddispvtx->GetN(); iS++) {
300	Double_t eta1 = fmddispvtx->GetX()[iS];
301	Double_t dispmult = fmddispvtx->GetY()[iS];
302	Double_t dispmulterr = fmddispvtx->GetErrorY(iS);
303
304	if(TMath::Abs(eta-eta1) >= 0.001) continue;
305
306	Double_t corr = nommult/dispmult;
307	Double_t rd = dispmulterr/dispmult;
308	Double_t rn = nommulterr/nommult;
309	Double_t error = (1/corr)*TMath::Sqrt(TMath::Power(rd,2) +
310	TMath::Power(rn,2));
311	corrcent[iC]->SetPoint(nPoints,eta,corr);
312	corrcent[iC]->SetPointError(nPoints,etaerr,error);
313	nPoints++;
314
315	}
316	//std::cout<<eta<<" "<<nommult<<std::endl;
317	}
318	nMin = TMath::Min(nPoints, nMin);
319	}
320
321	// --- Calulate the average --------------------------------------
322	TGraphErrors* average = new TGraphErrors;
323	average->SetName(Form("average"));
324	average->SetTitle(Form("%2d%% - %2d%%", limits[0], limits[4]));
325	average->SetMarkerColor(colors[4]);
326	average->SetLineColor(colors[4]);
327	average->SetMarkerStyle(fmdfmd ? 20 : 24);
328	average->SetMarkerSize(fmdfmd ? 1 : 1.2);
329	average->SetFillStyle(0);
330	average->SetFillColor(0);
331	mg->Add(average);
332
333	for(Int_t iA = 0; iA < nMin; iA++) {
334	Double_t mean = 0;
335	Double_t error2 = 0;
336	Double_t sumw2 = 0;
337
338	// Loop over centralities
339	for(Int_t iC=0; iC<4; iC++) {
340	Double_t eta = corrcent[iC]->GetX()[iA];
341	Double_t corr = corrcent[iC]->GetY()[iA];
342	Double_t err = corrcent[iC]->GetErrorY(iA);
343	Double_t err2 = err * err;
344	Double_t w = 1 / err2;
345	sumw2 += w;
346	mean += w * corr;
347	}
348	mean /= sumw2;
349	Double_t error = TMath::Sqrt(1. / sumw2);
350	average->SetPoint(iA,eta,mean);
351	average->SetPointError(iA,0.125,error);
352	}
353	if (fmdfmd) afmdfmd = average;
354	else afmdfull = average;
355
356	// --- Calculate ratios ------------------------------------------
357	TMultiGraph* ratios = new TMultiGraph;
358	ratios->SetName("ratios");
359	for(Int_t iC=0; iC<4; iC++) {
360	Int_t cl = limits[iC];
361	Int_t ch = limits[iC+1];
362	TGraphErrors* r = Ratio(corrcent[iC], average);
363	r->SetName(Form("ratio%s",base.Data()));
364	ratios->Add(r);
365	}
366
367	// --- Store the result ------------------------------------------
368	TDirectory* d = out->mkdir(fmdfmd ? "fmdfmd" : "fmdfull",
369	Form("Empirical corrections %s",
370	fmdfmd ? "FMD/FMD" : "FULL/FMD"));
371	d->cd();
372	for(Int_t p=0; p<4; p++) corrcent[p]->Write();
373	average->Write("average");
374	mg->Write("all");
375	ratios->Write("ratios");
376	mgref->Write("satellite");
377	out->cd();
378	if (mgfmdnom) mgfmdnom->Write("nominal");
379	out->ls();
380	} // End method loop
381
382	TGraphErrors* r = Ratio(afmdfull, afmdfmd);
383	r->Write("ratio");
384
385	out->Close();
386	}
387	//
388	// EOF
389	//
390