[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / FitELoss.C

/**
 * Test script to fit the energy loss spectra 
 * 
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
#ifndef __CINT__
# include "AliForwardUtil.h"
# include <TFile.h>
# include <TList.h>
# include <TH1.h>
# include <TF1.h>
# include <TFitResult.h>
# include <TMath.h>
# include <TStyle.h>
# include <TArrow.h>
# include <TCanvas.h>
#else
class TCanvas;
class TFile;
class TH1;
class TList;
class TF1;
#endif

//____________________________________________________________________
/** 
 * 
 * 
 * @param ef 
 * @param d 
 * @param r 
 * @param etabin 
 * 
 * @return 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
TH1* GetEDist(TList* ef, UShort_t d, Char_t r, UShort_t etabin)
{
  TList* dL = static_cast<TList*>(ef->FindObject(Form("FMD%d%c",d,r)));
  if (!dL) {
    Error("GetEDist", "Couldn't find list FMD%d%c",d,r);
    ef->ls();
    return 0;
  }
  if (etabin > 999) {
    TH1* hist = static_cast<TH1*>(dL->FindObject(Form("FMD%d%c_edist",d,r)));
    if (hist) {
      Error("GetEDist", "Couldn't find EDists histogram for FMD%d%c",d,r);
      return 0;
    }
  }
      
  TList* edL = static_cast<TList*>(dL->FindObject("EDists"));
  if (!edL) {
    Error("GetEDist", "Couldn't find list EDists for FMD%d%c",d,r);
    return 0;
  }
  
  TH1* hist = static_cast<TH1*>(edL->FindObject(Form("FMD%d%c_etabin%03d", 
						     d, r, etabin)));
  if (!hist) {
    Error("GetEDist", "Couldn't find histogra FMD%d%c_etabin%03d",
	  d,r, etabin);
    return 0;
  }
  
  return hist;
}

//____________________________________________________________________
/** 
 * 
 * 
 * @param ef 
 * @param d 
 * @param r 
 * @param eta 
 * 
 * @return  
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
*/
TH1* GetEDist(TList* ef, UShort_t d, Char_t r, Float_t eta)
{
  if (!ef) { 
    Error("GetEDist", "EF not set");
    return 0;
  }
  TAxis* etaAxis = static_cast<TAxis*>(ef->FindObject("etaAxis"));
  if (!etaAxis) { 
    Error("GetEDist", "Couldn't find eta axis in list");
    return 0;
  }

  UShort_t bin = etaAxis->FindBin(eta);
  if (bin <= 0 || bin > etaAxis->GetNbins()) { 
    Error("GetEDist", "eta=%f out of range [%f,%f] - getting ring histo", 
	  eta, etaAxis->GetXmin(), etaAxis->GetXmax());
    return GetEDist(ef, d, r, UShort_t(1000));
  }

  return GetEDist(ef, d, r, bin);
}

//____________________________________________________________________
/** 
 * 
 * 
 * @param file 
 * 
 * @return 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
TList* GetEF(TFile* file) 
{
  TList* forward = static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
  if (!forward) {
    Error("GetEF", "Failed to get list PWG2forwardDnDeta/Forward from file");
    return 0;
  }
  TList* ef = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
  if (!ef) {
    Error("GetEF", "Failed to get energy fitter list");
    return 0;
  }
  
  return ef;
}

//____________________________________________________________________
TList* ef = 0;

//____________________________________________________________________
/** 
 * 
 * 
 * 
 * @return 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
TList*  CheckEF()
{
  if (ef) return ef;
  
  TFile* file = TFile::Open("AnalysisResults.root", "READ");
  if (!file) { 
    Error("Fit1", "Failed to open file");
    return 0;
  }
  return GetEF(file);
}

//____________________________________________________________________
TCanvas* c = 0;

//____________________________________________________________________
/** 
 * 
 * 
 * 
 * @return 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
TCanvas* CheckC()
{
  if (c) return c;

  gStyle->SetOptFit(1111);
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderSize(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetCanvasDefW(800);
  gStyle->SetCanvasDefH(800);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleStyle(0);
  gStyle->SetStatBorderSize(1);
  gStyle->SetStatColor(0);
  gStyle->SetStatStyle(0);
  gStyle->SetStatX(0.95);
  gStyle->SetStatY(0.95);
  gStyle->SetStatW(0.15);
  gStyle->SetStatH(0.15);

  c = new TCanvas("c", "c");
  c->SetLogy();

  return c;
}

//____________________________________________________________________
/** 
 * 
 * 
 * @param f 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
void PrintFit(TF1* f)
{
  Int_t    nu     = f->GetNDF();
  Double_t chi2   = f->GetChisquare();
  Double_t chi2nu = (nu > 0 ? chi2/nu : 0);
  Printf("%s: chi^2/nu=%f/%d=%f [%f,%f]", 
	 f->GetName(), chi2, nu, chi2nu,
	 f->GetXmin(), f->GetXmax());
  for (Int_t i = 0; i < f->GetNpar(); i++) { 
    Double_t v = f->GetParameter(i);
    Double_t e = f->GetParError(i);
    Double_t r = 100*(v == 0 ? 1 : e / v);
    Printf("%32s = %14.7f +/- %14.7f (%5.1f%%)",f->GetParName(i),v,e,r);
  }
}
  
//____________________________________________________________________
/** 
 * 
 * 
 * @param n 
 * @param d 
 * @param r 
 * @param eta 
 *
 * @deprecated
 * This is a simple test script 
 *
 * @ingroup pwg2_forward_analysis_scripts
 */
void FitELoss(Int_t n, UShort_t d, Char_t r, Float_t eta)
{
  TList* ef1 = CheckEF();
  TCanvas* c1 = CheckC();
  if (!ef1) return;
  if (!c1)  return;

  TH1* dist = GetEDist(ef1, d, r, eta);
  if (!dist) return;

  AliForwardUtil::ELossFitter f(0.4, 10, 4);
  TF1* landau1 = new TF1(*f.Fit1Particle(dist, 0));
  landau1->SetName("Landau1");
  PrintFit(landau1);
  TF1* landaun = new TF1(*f.FitNParticle(dist, n, 0)); 
  landau1->SetName(Form("Landau%d", n));
  PrintFit(landaun);
  landau1->SetRange(0,10);
  landaun->SetRange(0,10);
  landau1->SetLineWidth(4);
  landaun->SetLineWidth(4);
  landau1->SetLineColor(kBlack);
  landaun->SetLineColor(kBlack);

  dist->GetListOfFunctions()->Clear();
  dist->GetListOfFunctions()->Add(landau1);
  dist->GetListOfFunctions()->Add(landaun);
  dist->GetListOfFunctions()->ls();
  dist->Draw();
  landau1->Draw("same");
  landaun->Draw("same");

  Double_t mp = landaun->GetParameter(1);
  Double_t xi = landaun->GetParameter(2);
  for (Int_t i  = 1; i <= n; i++) { 
    Double_t x  = i * (mp + xi * TMath::Log(i));
    Double_t y  = landaun->Eval(x);
    Double_t y1 = y < 0.05 ? 1 : 0.01;
    TArrow* a = new TArrow(x,y1,x,y,0.03,"|>");
    Info("FitSteer", "Delta_{p,%d}=%f", i, x);
    a->SetLineWidth(2);
    a->SetAngle(30);
    a->Draw();
  }

  c1->cd();
}
//
// EOF
//
Commit	Line	Data
ab0f914c	1	/**
	2	* Test script to fit the energy loss spectra
	3	*
1c762251	4	* @deprecated
	5	* This is a simple test script
	6	*
ab0f914c	7	* @ingroup pwg2_forward_analysis_scripts
ab0f914c	8	*/
c389303e	9	#ifndef __CINT__
	10	# include "AliForwardUtil.h"
	11	# include <TFile.h>
	12	# include <TList.h>
	13	# include <TH1.h>
	14	# include <TF1.h>
	15	# include <TFitResult.h>
	16	# include <TMath.h>
	17	# include <TStyle.h>
	18	# include <TArrow.h>
	19	# include <TCanvas.h>
	20	#else
	21	class TCanvas;
	22	class TFile;
	23	class TH1;
	24	class TList;
	25	class TF1;
	26	#endif
	27
	28	//____________________________________________________________________
1c762251	29	/**
	30	*
	31	*
	32	* @param ef
	33	* @param d
	34	* @param r
	35	* @param etabin
	36	*
	37	* @return
	38	*
	39	* @deprecated
	40	* This is a simple test script
	41	*
	42	* @ingroup pwg2_forward_analysis_scripts
	43	*/
c389303e	44	TH1* GetEDist(TList* ef, UShort_t d, Char_t r, UShort_t etabin)
	45	{
	46	TList* dL = static_cast<TList*>(ef->FindObject(Form("FMD%d%c",d,r)));
	47	if (!dL) {
	48	Error("GetEDist", "Couldn't find list FMD%d%c",d,r);
	49	ef->ls();
	50	return 0;
	51	}
	52	if (etabin > 999) {
	53	TH1* hist = static_cast<TH1*>(dL->FindObject(Form("FMD%d%c_edist",d,r)));
	54	if (hist) {
	55	Error("GetEDist", "Couldn't find EDists histogram for FMD%d%c",d,r);
	56	return 0;
	57	}
	58	}
	59
	60	TList* edL = static_cast<TList*>(dL->FindObject("EDists"));
	61	if (!edL) {
	62	Error("GetEDist", "Couldn't find list EDists for FMD%d%c",d,r);
	63	return 0;
	64	}
	65
	66	TH1* hist = static_cast<TH1*>(edL->FindObject(Form("FMD%d%c_etabin%03d",
	67	d, r, etabin)));
	68	if (!hist) {
	69	Error("GetEDist", "Couldn't find histogra FMD%d%c_etabin%03d",
	70	d,r, etabin);
	71	return 0;
	72	}
	73
	74	return hist;
	75	}
	76
	77	//____________________________________________________________________
1c762251	78	/**
	79	*
	80	*
	81	* @param ef
	82	* @param d
	83	* @param r
	84	* @param eta
	85	*
	86	* @return
	87	*
	88	* @deprecated
	89	* This is a simple test script
	90	*
	91	* @ingroup pwg2_forward_analysis_scripts
	92	*/
c389303e	93	TH1* GetEDist(TList* ef, UShort_t d, Char_t r, Float_t eta)
	94	{
	95	if (!ef) {
	96	Error("GetEDist", "EF not set");
	97	return 0;
	98	}
	99	TAxis* etaAxis = static_cast<TAxis*>(ef->FindObject("etaAxis"));
	100	if (!etaAxis) {
	101	Error("GetEDist", "Couldn't find eta axis in list");
	102	return 0;
	103	}
	104
	105	UShort_t bin = etaAxis->FindBin(eta);
	106	if (bin <= 0 \|\| bin > etaAxis->GetNbins()) {
	107	Error("GetEDist", "eta=%f out of range [%f,%f] - getting ring histo",
	108	eta, etaAxis->GetXmin(), etaAxis->GetXmax());
	109	return GetEDist(ef, d, r, UShort_t(1000));
	110	}
	111
	112	return GetEDist(ef, d, r, bin);
	113	}
	114
	115	//____________________________________________________________________
1c762251	116	/**
	117	*
	118	*
	119	* @param file
	120	*
	121	* @return
	122	*
	123	* @deprecated
	124	* This is a simple test script
	125	*
	126	* @ingroup pwg2_forward_analysis_scripts
	127	*/
c389303e	128	TList* GetEF(TFile* file)
	129	{
	130	TList* forward = static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
	131	if (!forward) {
	132	Error("GetEF", "Failed to get list PWG2forwardDnDeta/Forward from file");
	133	return 0;
	134	}
	135	TList* ef = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
	136	if (!ef) {
	137	Error("GetEF", "Failed to get energy fitter list");
	138	return 0;
	139	}
	140
	141	return ef;
	142	}
	143
	144	//____________________________________________________________________
	145	TList* ef = 0;
	146
	147	//____________________________________________________________________
1c762251	148	/**
	149	*
	150	*
	151	*
	152	* @return
	153	*
	154	* @deprecated
	155	* This is a simple test script
	156	*
	157	* @ingroup pwg2_forward_analysis_scripts
	158	*/
c389303e	159	TList* CheckEF()
	160	{
	161	if (ef) return ef;
	162
	163	TFile* file = TFile::Open("AnalysisResults.root", "READ");
	164	if (!file) {
	165	Error("Fit1", "Failed to open file");
	166	return 0;
	167	}
	168	return GetEF(file);
	169	}
	170
	171	//____________________________________________________________________
	172	TCanvas* c = 0;
	173
	174	//____________________________________________________________________
1c762251	175	/**
	176	*
	177	*
	178	*
	179	* @return
	180	*
	181	* @deprecated
	182	* This is a simple test script
	183	*
	184	* @ingroup pwg2_forward_analysis_scripts
	185	*/
c389303e	186	TCanvas* CheckC()
	187	{
	188	if (c) return c;
	189
	190	gStyle->SetOptFit(1111);
	191	gStyle->SetCanvasColor(0);
	192	gStyle->SetCanvasBorderSize(0);
	193	gStyle->SetCanvasBorderMode(0);
	194	gStyle->SetCanvasDefW(800);
	195	gStyle->SetCanvasDefH(800);
	196	gStyle->SetPadTopMargin(0.05);
	197	gStyle->SetPadRightMargin(0.05);
	198	gStyle->SetTitleBorderSize(0);
	199	gStyle->SetTitleFillColor(0);
	200	gStyle->SetTitleStyle(0);
	201	gStyle->SetStatBorderSize(1);
	202	gStyle->SetStatColor(0);
	203	gStyle->SetStatStyle(0);
	204	gStyle->SetStatX(0.95);
	205	gStyle->SetStatY(0.95);
	206	gStyle->SetStatW(0.15);
	207	gStyle->SetStatH(0.15);
	208
	209	c = new TCanvas("c", "c");
	210	c->SetLogy();
	211
	212	return c;
	213	}
	214
	215	//____________________________________________________________________
1c762251	216	/**
	217	*
	218	*
	219	* @param f
	220	*
	221	* @deprecated
	222	* This is a simple test script
	223	*
	224	* @ingroup pwg2_forward_analysis_scripts
	225	*/
c389303e	226	void PrintFit(TF1* f)
	227	{
	228	Int_t nu = f->GetNDF();
	229	Double_t chi2 = f->GetChisquare();
	230	Double_t chi2nu = (nu > 0 ? chi2/nu : 0);
	231	Printf("%s: chi^2/nu=%f/%d=%f [%f,%f]",
	232	f->GetName(), chi2, nu, chi2nu,
	233	f->GetXmin(), f->GetXmax());
	234	for (Int_t i = 0; i < f->GetNpar(); i++) {
	235	Double_t v = f->GetParameter(i);
	236	Double_t e = f->GetParError(i);
	237	Double_t r = 100*(v == 0 ? 1 : e / v);
	238	Printf("%32s = %14.7f +/- %14.7f (%5.1f%%)",f->GetParName(i),v,e,r);
	239	}
	240	}
	241
	242	//____________________________________________________________________
1c762251	243	/**
	244	*
	245	*
	246	* @param n
	247	* @param d
	248	* @param r
	249	* @param eta
	250	*
	251	* @deprecated
	252	* This is a simple test script
	253	*
	254	* @ingroup pwg2_forward_analysis_scripts
	255	*/
c389303e	256	void FitELoss(Int_t n, UShort_t d, Char_t r, Float_t eta)
	257	{
	258	TList* ef1 = CheckEF();
	259	TCanvas* c1 = CheckC();
	260	if (!ef1) return;
	261	if (!c1) return;
	262
	263	TH1* dist = GetEDist(ef1, d, r, eta);
	264	if (!dist) return;
	265
	266	AliForwardUtil::ELossFitter f(0.4, 10, 4);
	267	TF1* landau1 = new TF1(*f.Fit1Particle(dist, 0));
	268	landau1->SetName("Landau1");
	269	PrintFit(landau1);
	270	TF1* landaun = new TF1(*f.FitNParticle(dist, n, 0));
	271	landau1->SetName(Form("Landau%d", n));
	272	PrintFit(landaun);
	273	landau1->SetRange(0,10);
	274	landaun->SetRange(0,10);
	275	landau1->SetLineWidth(4);
	276	landaun->SetLineWidth(4);
	277	landau1->SetLineColor(kBlack);
	278	landaun->SetLineColor(kBlack);
	279
	280	dist->GetListOfFunctions()->Clear();
	281	dist->GetListOfFunctions()->Add(landau1);
	282	dist->GetListOfFunctions()->Add(landaun);
	283	dist->GetListOfFunctions()->ls();
	284	dist->Draw();
	285	landau1->Draw("same");
	286	landaun->Draw("same");
	287
	288	Double_t mp = landaun->GetParameter(1);
	289	Double_t xi = landaun->GetParameter(2);
	290	for (Int_t i = 1; i <= n; i++) {
	291	Double_t x = i * (mp + xi * TMath::Log(i));
	292	Double_t y = landaun->Eval(x);
	293	Double_t y1 = y < 0.05 ? 1 : 0.01;
	294	TArrow* a = new TArrow(x,y1,x,y,0.03,"\|>");
	295	Info("FitSteer", "Delta_{p,%d}=%f", i, x);
	296	a->SetLineWidth(2);
	297	a->SetAngle(30);
	298	a->Draw();
	299	}
	300
	301	c1->cd();
	302	}
1c762251	303	//
	304	// EOF
	305	//