[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis / scripts / ViewEDCor.C

//____________________________________________________________________
AliFMDAnaCalibEnergyDistribution*
ViewEDCor(const char* filename)
{
  // gSystem->Load("libANALYSIS");
  // gSystem->Load("libANALYSISalice");
  // gSystem->Load("libPWG0base");
  // gSystem->Load("libPWGLFforward");

  AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();
  param->SetEnergy(900.F);
  param->SetCollisionSystem("pp");
  param->SetMagField(+5.F);
  param->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);

  TFile* file = TFile::Open(filename, "READ");
  if (!file) { 
    Error("ViewEDCor", "Cannot open file %s", filename);
    return 0;
  }

  AliFMDAnaCalibEnergyDistribution* bg = 
    static_cast<AliFMDAnaCalibEnergyDistribution*>(file->Get("energydistributions"));
  
  return bg;
}


//____________________________________________________________________
struct RingHistos 
{
  TH1D* MakeHisto(const char* name, const char* title, 
		  const TAxis* etaAxis)
  {
    TH1D* ret = new TH1D(Form(name, fD, fR), 
			 Form(title, fD, fR), 
			 etaAxis->GetNbins(), 
			 etaAxis->GetXmin(), 
			 etaAxis->GetXmax());
    ret->SetDirectory(0);
    ret->SetLineColor(fC);
    ret->SetMarkerColor(fC);
    ret->SetMarkerStyle(20+2*fD+(fR == 'I' || fR == 'i' ? 0 : 1));
    ret->SetFillColor(fC);
    ret->SetFillStyle(3001);
    ret->SetXTitle("#eta");
    ret->SetYTitle(title);
    ret->SetStats(0);
    return ret;
  }
  
  
  RingHistos(UShort_t d, Char_t r, const TAxis* etaAxis, Int_t c) 
    : fD(d), fR(r), fC(c) 
  {
    fHists.AddAt(MakeHisto("c%d%c",     "Constant FMD%d%c",      etaAxis),0);
    fHists.AddAt(MakeHisto("mpv%d%c",   "MPV FMD%d%c",           etaAxis),1);
    fHists.AddAt(MakeHisto("w%d%c",     "Width FMD%d%c",         etaAxis),2);
    fHists.AddAt(MakeHisto("alpha%d%c", "#alpha FMD%d%c",        etaAxis),3);
    fHists.AddAt(MakeHisto("beta%d%c",  "#beta FMD%d%c",         etaAxis),4);
    fHists.AddAt(MakeHisto("chi%d%c",   "#chi^{2}/#nu FMD%d%c",  etaAxis),5);
    fHists.AddAt(MakeHisto("acc%d%c",   "Acceptance in FMD%d%c", etaAxis),6);
  }
  void Fill(Int_t iEta, Int_t p, Double_t v, Double_t e=0) 
  {
    TH1D* h = Get(p);
    if (!h) return;
    h->SetBinContent(iEta, v);
    h->SetBinError(iEta, e);
  }
  TH1D* Get(Int_t p) const
  {
    return static_cast<TH1D*>(fHists.At(p));
  }
  UShort_t fD;
  Char_t   fR;
  Int_t    fC;
  TObjArray fHists;
};

//____________________________________________________________________
struct All 
{
  All(const TAxis* etaAxis)
  {
    fFMD1i = new RingHistos(1,'i', etaAxis, kRed+1);
    fFMD2i = new RingHistos(2,'i', etaAxis, kGreen+1);
    fFMD2o = new RingHistos(2,'o', etaAxis, kGreen+7);
    fFMD3i = new RingHistos(3,'i', etaAxis, kBlue+1);
    fFMD3o = new RingHistos(3,'o', etaAxis, kBlue+7);
  }
  RingHistos* Get(UShort_t d, Char_t r) const
  {
    switch (d) { 
    case 1:    return fFMD1i;
    case 2:    return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
    case 3:    return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
    }
    return 0;
  }
  RingHistos* fFMD1i; 
  RingHistos* fFMD2i; 
  RingHistos* fFMD2o; 
  RingHistos* fFMD3i; 
  RingHistos* fFMD3o; 
};

//____________________________________________________________________
PlotEDCor(const char* filename)
{
  AliFMDAnaCalibEnergyDistribution* ed = ViewEDCor(filename);

  AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();

  TAxis* etaAxis = param->GetBackgroundCorrection(1,'I',4)->GetXaxis();
  All* all = new All(etaAxis);

  for (Int_t iEta = 1; iEta < etaAxis->GetNbins(); iEta++) { 
    Float_t eta = etaAxis->GetBinCenter(iEta);
    Info("PlotEDCor", "eta bin %3d, eta=%f", iEta, eta);

    for (UShort_t d=1; d <= 3; d++) { 
      UShort_t nq = d == 1 ? 1 : 2;
      for (UShort_t q = 0; q < nq; q++) { 
	Char_t r = q == 0 ? 'I' : 'O';
	RingHistos* hists = all->Get(d, r);
	if (!hists) continue;
	
	TH1F* hed = ed->GetEnergyDistribution(d, r, eta);
	if (!hed) continue;

	if (hed->GetEntries() > 0) 
	  hists->Fill(iEta, 6, 1);
	
	TF1* fed = hed->GetFunction("FMDfitFunc");
	if (!fed) continue;

	Int_t npar = fed->GetNpar();
	for (Int_t i = 0; i < npar; i++) { 
	  Double_t par = fed->GetParameter(i);
	  Double_t err = fed->GetParError(i);
	  hists->Fill(iEta, i, par, err);
	}
	Double_t chi2 = fed->GetChisquare();
	Double_t ndf  = fed->GetNDF();
	if (ndf != 0) 
	  hists->Fill(iEta, 5, chi2/ndf);
      }
    }
  }
	  
  TCanvas* c = new TCanvas("c", "c", 800, 800);
  c->SetFillColor(0);
  c->SetBorderMode(0);
  c->SetBorderSize(0);
  c->Divide(1,7,0,0);
  
  THStack* stack = 0;
  
  const char* yTitles[] = { "Constant", "#Delta_{p}", "w", 
			    "#alpha", "#beta", "#chi^2/NDF", 
			    "Acceptance" };
  for (Int_t i = 0; i < 7; i++) { 
    TVirtualPad* p = c->cd(i+1);
    p->SetFillColor(0);
    p->SetFillStyle(0);
    THStack* stack = new THStack;
    stack->Add(all->Get(1,'i')->Get(i));
    stack->Add(all->Get(2,'i')->Get(i));
    stack->Add(all->Get(2,'o')->Get(i));
    stack->Add(all->Get(3,'i')->Get(i));
    stack->Add(all->Get(3,'o')->Get(i));
    stack->Draw("nostack");
    stack->GetHistogram()->SetYTitle(yTitles[i]);
    stack->GetHistogram()->SetXTitle("#eta");
    TAxis* yaxis = stack->GetHistogram()->GetYaxis();
    yaxis->SetTitleSize(0.15);
    yaxis->SetLabelSize(0.08);
    yaxis->SetTitleOffset(0.35);
    yaxis->SetNdivisions(10);
    TAxis* xaxis = stack->GetHistogram()->GetXaxis();
    xaxis->SetTitleSize(0.15);
    xaxis->SetLabelSize(0.08);
    xaxis->SetTitleOffset(0.35);
    xaxis->SetNdivisions(320);
    p->cd();
  }
  c->cd();
  c->Print("energyDist.png");
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
45535892	1	//____________________________________________________________________
	2	AliFMDAnaCalibEnergyDistribution*
	3	ViewEDCor(const char* filename)
	4	{
	5	// gSystem->Load("libANALYSIS");
	6	// gSystem->Load("libANALYSISalice");
	7	// gSystem->Load("libPWG0base");
bd6f5206	8	// gSystem->Load("libPWGLFforward");
45535892	9
	10	AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();
	11	param->SetEnergy(900.F);
	12	param->SetCollisionSystem("pp");
	13	param->SetMagField(+5.F);
	14	param->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
	15
	16	TFile* file = TFile::Open(filename, "READ");
	17	if (!file) {
	18	Error("ViewEDCor", "Cannot open file %s", filename);
	19	return 0;
	20	}
	21
	22	AliFMDAnaCalibEnergyDistribution* bg =
	23	static_cast<AliFMDAnaCalibEnergyDistribution*>(file->Get("energydistributions"));
	24
	25	return bg;
	26	}
	27
	28
	29	//____________________________________________________________________
	30	struct RingHistos
	31	{
	32	TH1D* MakeHisto(const char* name, const char* title,
	33	const TAxis* etaAxis)
	34	{
	35	TH1D* ret = new TH1D(Form(name, fD, fR),
	36	Form(title, fD, fR),
	37	etaAxis->GetNbins(),
	38	etaAxis->GetXmin(),
	39	etaAxis->GetXmax());
	40	ret->SetDirectory(0);
	41	ret->SetLineColor(fC);
	42	ret->SetMarkerColor(fC);
	43	ret->SetMarkerStyle(20+2*fD+(fR == 'I' \|\| fR == 'i' ? 0 : 1));
	44	ret->SetFillColor(fC);
	45	ret->SetFillStyle(3001);
	46	ret->SetXTitle("#eta");
	47	ret->SetYTitle(title);
	48	ret->SetStats(0);
	49	return ret;
	50	}
	51
	52
	53	RingHistos(UShort_t d, Char_t r, const TAxis* etaAxis, Int_t c)
	54	: fD(d), fR(r), fC(c)
	55	{
	56	fHists.AddAt(MakeHisto("c%d%c", "Constant FMD%d%c", etaAxis),0);
	57	fHists.AddAt(MakeHisto("mpv%d%c", "MPV FMD%d%c", etaAxis),1);
	58	fHists.AddAt(MakeHisto("w%d%c", "Width FMD%d%c", etaAxis),2);
	59	fHists.AddAt(MakeHisto("alpha%d%c", "#alpha FMD%d%c", etaAxis),3);
	60	fHists.AddAt(MakeHisto("beta%d%c", "#beta FMD%d%c", etaAxis),4);
	61	fHists.AddAt(MakeHisto("chi%d%c", "#chi^{2}/#nu FMD%d%c", etaAxis),5);
	62	fHists.AddAt(MakeHisto("acc%d%c", "Acceptance in FMD%d%c", etaAxis),6);
	63	}
	64	void Fill(Int_t iEta, Int_t p, Double_t v, Double_t e=0)
	65	{
	66	TH1D* h = Get(p);
	67	if (!h) return;
	68	h->SetBinContent(iEta, v);
	69	h->SetBinError(iEta, e);
	70	}
	71	TH1D* Get(Int_t p) const
	72	{
73	return static_cast<TH1D*>(fHists.At(p));
74	}
75	UShort_t fD;
76	Char_t fR;
77	Int_t fC;
78	TObjArray fHists;
79	};
80
81	//____________________________________________________________________
82	struct All
83	{
84	All(const TAxis* etaAxis)
85	{
86	fFMD1i = new RingHistos(1,'i', etaAxis, kRed+1);
87	fFMD2i = new RingHistos(2,'i', etaAxis, kGreen+1);
88	fFMD2o = new RingHistos(2,'o', etaAxis, kGreen+7);
89	fFMD3i = new RingHistos(3,'i', etaAxis, kBlue+1);
90	fFMD3o = new RingHistos(3,'o', etaAxis, kBlue+7);
91	}
92	RingHistos* Get(UShort_t d, Char_t r) const
93	{
94	switch (d) {
95	case 1: return fFMD1i;
96	case 2: return (r == 'I' \|\| r == 'i' ? fFMD2i : fFMD2o);
97	case 3: return (r == 'I' \|\| r == 'i' ? fFMD3i : fFMD3o);
98	}
99	return 0;
100	}
101	RingHistos* fFMD1i;
102	RingHistos* fFMD2i;
103	RingHistos* fFMD2o;
104	RingHistos* fFMD3i;
105	RingHistos* fFMD3o;
106	};
107
108	//____________________________________________________________________
109	PlotEDCor(const char* filename)
110	{
111	AliFMDAnaCalibEnergyDistribution* ed = ViewEDCor(filename);
112
113	AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();
114
115	TAxis* etaAxis = param->GetBackgroundCorrection(1,'I',4)->GetXaxis();
116	All* all = new All(etaAxis);
117
118	for (Int_t iEta = 1; iEta < etaAxis->GetNbins(); iEta++) {
119	Float_t eta = etaAxis->GetBinCenter(iEta);
120	Info("PlotEDCor", "eta bin %3d, eta=%f", iEta, eta);
121
122	for (UShort_t d=1; d <= 3; d++) {
123	UShort_t nq = d == 1 ? 1 : 2;
124	for (UShort_t q = 0; q < nq; q++) {
125	Char_t r = q == 0 ? 'I' : 'O';
126	RingHistos* hists = all->Get(d, r);
127	if (!hists) continue;
128
129	TH1F* hed = ed->GetEnergyDistribution(d, r, eta);
130	if (!hed) continue;
131
132	if (hed->GetEntries() > 0)
133	hists->Fill(iEta, 6, 1);
134
135	TF1* fed = hed->GetFunction("FMDfitFunc");
136	if (!fed) continue;
137
138	Int_t npar = fed->GetNpar();
139	for (Int_t i = 0; i < npar; i++) {
140	Double_t par = fed->GetParameter(i);
141	Double_t err = fed->GetParError(i);
142	hists->Fill(iEta, i, par, err);
143	}
144	Double_t chi2 = fed->GetChisquare();
145	Double_t ndf = fed->GetNDF();
146	if (ndf != 0)
147	hists->Fill(iEta, 5, chi2/ndf);
148	}
149	}
150	}
151
152	TCanvas* c = new TCanvas("c", "c", 800, 800);
153	c->SetFillColor(0);
154	c->SetBorderMode(0);
155	c->SetBorderSize(0);
156	c->Divide(1,7,0,0);
157
158	THStack* stack = 0;
159
160	const char* yTitles[] = { "Constant", "#Delta_{p}", "w",
161	"#alpha", "#beta", "#chi^2/NDF",
162	"Acceptance" };
163	for (Int_t i = 0; i < 7; i++) {
164	TVirtualPad* p = c->cd(i+1);
165	p->SetFillColor(0);
166	p->SetFillStyle(0);
167	THStack* stack = new THStack;
168	stack->Add(all->Get(1,'i')->Get(i));
169	stack->Add(all->Get(2,'i')->Get(i));
170	stack->Add(all->Get(2,'o')->Get(i));
171	stack->Add(all->Get(3,'i')->Get(i));
172	stack->Add(all->Get(3,'o')->Get(i));
173	stack->Draw("nostack");
174	stack->GetHistogram()->SetYTitle(yTitles[i]);
175	stack->GetHistogram()->SetXTitle("#eta");
176	TAxis* yaxis = stack->GetHistogram()->GetYaxis();
177	yaxis->SetTitleSize(0.15);
178	yaxis->SetLabelSize(0.08);
179	yaxis->SetTitleOffset(0.35);
180	yaxis->SetNdivisions(10);
181	TAxis* xaxis = stack->GetHistogram()->GetXaxis();
182	xaxis->SetTitleSize(0.15);
183	xaxis->SetLabelSize(0.08);
184	xaxis->SetTitleOffset(0.35);
185	xaxis->SetNdivisions(320);
186	p->cd();
187	}
188	c->cd();
189	c->Print("energyDist.png");
190	}
191	//____________________________________________________________________
192	//
193	// EOF
194	//
195