[u/mrichter/AliRoot.git] / FMD / scripts / DrawESD.C

//____________________________________________________________________
//
// $Id$
//
// Script that contains a class to draw eloss from hits, versus ADC
// counts from digits, using the AliFMDInputHits class in the util library. 
//
// It draws the energy loss versus the p/(mq^2).  It can be overlayed
// with the Bethe-Bloc curve to show how the simulation behaves
// relative to the expected. 
//
// Use the script `Compile.C' to compile this class using ACLic. 
//
#include <TH1D.h>
#include <AliFMDHit.h>
#include <AliFMDDigit.h>
#include <AliFMDInput.h>
#include <AliFMDUShortMap.h>
#include <AliFMDFloatMap.h>
#include <AliFMDRecPoint.h>
#include <AliESDFMD.h>
#include <AliLog.h>
#include <iostream>
#include <TStyle.h>
#include <TArrayF.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TF1.h>
#include <TSpectrum.h>
#include <TLegend.h>
#include <TLine.h>

/** @class DrawESD
    @brief Draw digit ADC versus Rec point mult
    @code 
    Root> .L Compile.C
    Root> Compile("DrawESD.C")
    Root> DrawESD c
    Root> c.Run();
    @endcode
    @ingroup FMD_script
 */
class DrawESD : public AliFMDInput
{
private:
  TH1D* fMult; // Histogram 
  const Double_t fCorr;
public:
  //__________________________________________________________________
  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
  {
    TArrayF bins(n+1);
    Float_t dp   = n / TMath::Log10(max / min);
    Float_t pmin = TMath::Log10(min);
    bins[0]      = min;
    for (Int_t i = 1; i < n+1; i++) {
      Float_t p = pmin + i / dp;
      bins[i]   = TMath::Power(10, p);
    }
    return bins;
  }
  //__________________________________________________________________
  DrawESD(Int_t n=420, Double_t mmin=-0.5, Double_t mmax=20.5) 
    : fCorr(1) // 0.68377 / 1.1)
  { 
    AddLoad(kESD);
    fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
    fMult->SetXTitle("Strip Multiplicity");
  }
  //__________________________________________________________________
  /** Begining of event
      @param ev Event number
      @return @c false on error */
  Bool_t Begin(Int_t ev) 
  {
    return AliFMDInput::Begin(ev);
  }
  //__________________________________________________________________
  Bool_t ProcessESD(UShort_t /* det */, 
		    Char_t   /* ring */, 
		    UShort_t /* sec */, 
		    UShort_t /* strip */, 
		    Float_t  /* eta */, 
		    Float_t  mult)
  {
    // Cache the energy loss 
    fMult->Fill(mult/fCorr);
    return kTRUE;
  }
  //__________________________________________________________________
  Bool_t Finish()
  {
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetOptFit(1111111);
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderSize(0);
    gStyle->SetPadColor(0);
    gStyle->SetPadBorderSize(0);
    
    TCanvas* c = new TCanvas("c", "C");
    c->cd();
    c->SetLogy();
    // fMult->GetXaxis()->SetRangeUser(0.2,4);
    // fMult->Sumw2();
    // fMult->Scale(1. / fMult->GetMaximum());
    fMult->SetStats(kFALSE);
    fMult->SetFillColor(2);
    fMult->SetFillStyle(3001);
    fMult->Draw();
    Double_t x1 = .75; // .8;  // .65 / fCorr;
    Double_t x2 = 1.3; // 1.7; // fCorr;
    Double_t x3 = 2.5; // 2.7; // fCorr;
    Double_t x4 = 3.7; // fCorr;
    TF1* l1 = new TF1("landau1", "landau", x1, x2);
    TF1* l2 = new TF1("landau2", "landau", x2, x3);
    // TF1* l3 = new TF1("landau3", "landau", x3, x4);
    TF1* f  = new TF1("user", "landau(0)+landau(3)", x1, x3);

    fMult->SetStats(kTRUE);
    fMult->GetXaxis()->SetRangeUser(0, 4);
    fMult->Fit(l1, "E0L", "", x1, x2);
    fMult->Fit(l2, "E0L", "", x2, x3);
    // fMult->Fit(l3, "E0L", "", x3, x4);
    f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
		   "A_{2}", "Mpv_{2}", "#sigma_{2}",
		   "A_{3}", "Mpv_{3}", "#sigma_{3}");
    f->SetParameters(l1->GetParameter(0), 
		     l1->GetParameter(1), 
		     l1->GetParameter(2), 
		     l2->GetParameter(0), 
		     l2->GetParameter(1), 
		     l2->GetParameter(2),
		     l1->GetParameter(0), 
		     l1->GetParameter(1) * 3, 
		     l1->GetParameter(2) * 3);
    fMult->Fit(f, "E0L", "", x1, x3);
    fMult->Fit(f, "MEL", "E1", x1, x3);

    TH1* h = static_cast<TH1*>(fMult->DrawClone("same hist"));

    // l1->SetLineStyle(2);
    l1->SetLineColor(3);
    l1->SetLineWidth(2);
    l1->SetRange(0, 4);
    l1->Draw("same");
    // l2->SetLineStyle(3);
    l2->SetLineColor(4);
    l2->SetLineWidth(2);
    l2->SetRange(0, 4);
    l2->Draw("same");
    // l3->SetLineStyle(3);
    // l3->SetLineWidth(2);
    // l3->SetRange(0, 5);
    // l3->Draw("same");
    f->SetLineWidth(2);
    f->SetRange(0, 4);
    f->Draw("same");

    TLegend* l = new TLegend(0.2, 0.6, .6, .89);
    l->AddEntry(l1, "1 particle Landau", "l");
    l->AddEntry(l2, "2 particle Landau", "l");
    l->AddEntry(f,  "1+2 particle Landau", "l");
    l->SetFillColor(0);
    l->Draw("same");


    c = new TCanvas("c2", "Landaus");
    c->SetLogy();
    fMult->DrawClone("axis");
    f->Draw("same");
    TF1* ll1 = new TF1("ll1", "landau", 0, 4);
    ll1->SetParameters(f->GetParameter(0), 
		       f->GetParameter(1), 
		       f->GetParameter(2));
    ll1->SetLineColor(3);
    ll1->Draw("same");
    TF1* ll2 = new TF1("ll2", "landau", 0, 4);
    ll2->SetParameters(f->GetParameter(3), 
		       f->GetParameter(4), 
		       f->GetParameter(5));
    ll2->SetLineColor(4);
    ll2->Draw("same");

    Double_t y1  = fMult->GetMinimum() * 1.1;
    Double_t y2  = fMult->GetMaximum() * .9;
    Double_t xc1 = f->GetParameter(1)-3*f->GetParameter(2);
    Double_t xc2 = f->GetParameter(4)-2*f->GetParameter(5);
    TLine* c1 = new TLine(xc1, y1, xc1, y2);
    c1->Draw("same");
    TLine* c2 = new TLine(xc2, y1, xc2, y2);
    c2->Draw("same");

    l = new TLegend(0.2, 0.6, .6, .89);
    l->AddEntry(ll1, "1 particle Landau", "l");
    l->AddEntry(ll2, "2 particle Landau", "l");
    l->AddEntry(f,  "1+2 particle Landau", "l");
    l->SetFillColor(0);
    l->Draw("same");

#if 0
    c = new TCanvas("s", "Spectrum");
    TSpectrum* s = new TSpectrum(16);
    h->GetXaxis()->SetRangeUser(0.3, 20);
    s->Search(h, .15, "", 0.01);
    c->Update();
    TH1* b = s->Background(h);
    b->SetFillColor(4);
    b->SetFillStyle(3001);
    b->Draw("same");
#endif

    return kTRUE;
  }

  ClassDef(DrawESD,0);
  
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
2b893216	1	//____________________________________________________________________
	2	//
	3	// $Id$
	4	//
	5	// Script that contains a class to draw eloss from hits, versus ADC
	6	// counts from digits, using the AliFMDInputHits class in the util library.
	7	//
	8	// It draws the energy loss versus the p/(mq^2). It can be overlayed
	9	// with the Bethe-Bloc curve to show how the simulation behaves
	10	// relative to the expected.
	11	//
	12	// Use the script `Compile.C' to compile this class using ACLic.
	13	//
	14	#include <TH1D.h>
	15	#include <AliFMDHit.h>
	16	#include <AliFMDDigit.h>
	17	#include <AliFMDInput.h>
	18	#include <AliFMDUShortMap.h>
	19	#include <AliFMDFloatMap.h>
	20	#include <AliFMDRecPoint.h>
	21	#include <AliESDFMD.h>
	22	#include <AliLog.h>
	23	#include <iostream>
	24	#include <TStyle.h>
	25	#include <TArrayF.h>
	26	#include <TCanvas.h>
	27	#include <TMath.h>
	28	#include <TF1.h>
	29	#include <TSpectrum.h>
	30	#include <TLegend.h>
	31	#include <TLine.h>
	32
	33	/** @class DrawESD
	34	@brief Draw digit ADC versus Rec point mult
	35	@code
	36	Root> .L Compile.C
	37	Root> Compile("DrawESD.C")
	38	Root> DrawESD c
	39	Root> c.Run();
	40	@endcode
	41	@ingroup FMD_script
	42	*/
	43	class DrawESD : public AliFMDInput
	44	{
	45	private:
	46	TH1D* fMult; // Histogram
	47	const Double_t fCorr;
	48	public:
	49	//__________________________________________________________________
	50	TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
	51	{
	52	TArrayF bins(n+1);
	53	Float_t dp = n / TMath::Log10(max / min);
	54	Float_t pmin = TMath::Log10(min);
	55	bins[0] = min;
	56	for (Int_t i = 1; i < n+1; i++) {
	57	Float_t p = pmin + i / dp;
	58	bins[i] = TMath::Power(10, p);
	59	}
	60	return bins;
	61	}
	62	//__________________________________________________________________
	63	DrawESD(Int_t n=420, Double_t mmin=-0.5, Double_t mmax=20.5)
	64	: fCorr(1) // 0.68377 / 1.1)
65	{
66	AddLoad(kESD);
67	fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
68	fMult->SetXTitle("Strip Multiplicity");
69	}
70	//__________________________________________________________________
71	/** Begining of event
72	@param ev Event number
73	@return @c false on error */
74	Bool_t Begin(Int_t ev)
75	{
76	return AliFMDInput::Begin(ev);
77	}
78	//__________________________________________________________________
79	Bool_t ProcessESD(UShort_t /* det */,
80	Char_t /* ring */,
81	UShort_t /* sec */,
82	UShort_t /* strip */,
83	Float_t /* eta */,
84	Float_t mult)
85	{
86	// Cache the energy loss
87	fMult->Fill(mult/fCorr);
88	return kTRUE;
89	}
90	//__________________________________________________________________
91	Bool_t Finish()
92	{
93	gStyle->SetPalette(1);
94	gStyle->SetOptTitle(0);
95	gStyle->SetOptFit(1111111);
96	gStyle->SetCanvasColor(0);
97	gStyle->SetCanvasBorderSize(0);
98	gStyle->SetPadColor(0);
99	gStyle->SetPadBorderSize(0);
100
101	TCanvas* c = new TCanvas("c", "C");
102	c->cd();
103	c->SetLogy();
104	// fMult->GetXaxis()->SetRangeUser(0.2,4);
105	// fMult->Sumw2();
106	// fMult->Scale(1. / fMult->GetMaximum());
107	fMult->SetStats(kFALSE);
108	fMult->SetFillColor(2);
109	fMult->SetFillStyle(3001);
110	fMult->Draw();
111	Double_t x1 = .75; // .8; // .65 / fCorr;
112	Double_t x2 = 1.3; // 1.7; // fCorr;
113	Double_t x3 = 2.5; // 2.7; // fCorr;
114	Double_t x4 = 3.7; // fCorr;
115	TF1* l1 = new TF1("landau1", "landau", x1, x2);
116	TF1* l2 = new TF1("landau2", "landau", x2, x3);
117	// TF1* l3 = new TF1("landau3", "landau", x3, x4);
118	TF1* f = new TF1("user", "landau(0)+landau(3)", x1, x3);
119
120	fMult->SetStats(kTRUE);
121	fMult->GetXaxis()->SetRangeUser(0, 4);
122	fMult->Fit(l1, "E0L", "", x1, x2);
123	fMult->Fit(l2, "E0L", "", x2, x3);
124	// fMult->Fit(l3, "E0L", "", x3, x4);
125	f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
126	"A_{2}", "Mpv_{2}", "#sigma_{2}",
127	"A_{3}", "Mpv_{3}", "#sigma_{3}");
128	f->SetParameters(l1->GetParameter(0),
129	l1->GetParameter(1),
130	l1->GetParameter(2),
131	l2->GetParameter(0),
132	l2->GetParameter(1),
133	l2->GetParameter(2),
134	l1->GetParameter(0),
135	l1->GetParameter(1) * 3,
136	l1->GetParameter(2) * 3);
137	fMult->Fit(f, "E0L", "", x1, x3);
138	fMult->Fit(f, "MEL", "E1", x1, x3);
139
140	TH1* h = static_cast<TH1*>(fMult->DrawClone("same hist"));
141
142	// l1->SetLineStyle(2);
143	l1->SetLineColor(3);
144	l1->SetLineWidth(2);
145	l1->SetRange(0, 4);
146	l1->Draw("same");
147	// l2->SetLineStyle(3);
148	l2->SetLineColor(4);
149	l2->SetLineWidth(2);
150	l2->SetRange(0, 4);
151	l2->Draw("same");
152	// l3->SetLineStyle(3);
153	// l3->SetLineWidth(2);
154	// l3->SetRange(0, 5);
155	// l3->Draw("same");
156	f->SetLineWidth(2);
157	f->SetRange(0, 4);
158	f->Draw("same");
159
160	TLegend* l = new TLegend(0.2, 0.6, .6, .89);
161	l->AddEntry(l1, "1 particle Landau", "l");
162	l->AddEntry(l2, "2 particle Landau", "l");
163	l->AddEntry(f, "1+2 particle Landau", "l");
164	l->SetFillColor(0);
165	l->Draw("same");
166
167
168	c = new TCanvas("c2", "Landaus");
169	c->SetLogy();
170	fMult->DrawClone("axis");
171	f->Draw("same");
172	TF1* ll1 = new TF1("ll1", "landau", 0, 4);
173	ll1->SetParameters(f->GetParameter(0),
174	f->GetParameter(1),
175	f->GetParameter(2));
176	ll1->SetLineColor(3);
177	ll1->Draw("same");
178	TF1* ll2 = new TF1("ll2", "landau", 0, 4);
179	ll2->SetParameters(f->GetParameter(3),
180	f->GetParameter(4),
181	f->GetParameter(5));
182	ll2->SetLineColor(4);
183	ll2->Draw("same");
184
185	Double_t y1 = fMult->GetMinimum() * 1.1;
186	Double_t y2 = fMult->GetMaximum() * .9;
187	Double_t xc1 = f->GetParameter(1)-3*f->GetParameter(2);
188	Double_t xc2 = f->GetParameter(4)-2*f->GetParameter(5);
189	TLine* c1 = new TLine(xc1, y1, xc1, y2);
190	c1->Draw("same");
191	TLine* c2 = new TLine(xc2, y1, xc2, y2);
192	c2->Draw("same");
193
194	l = new TLegend(0.2, 0.6, .6, .89);
195	l->AddEntry(ll1, "1 particle Landau", "l");
196	l->AddEntry(ll2, "2 particle Landau", "l");
197	l->AddEntry(f, "1+2 particle Landau", "l");
198	l->SetFillColor(0);
199	l->Draw("same");
200
201	#if 0
202	c = new TCanvas("s", "Spectrum");
203	TSpectrum* s = new TSpectrum(16);
204	h->GetXaxis()->SetRangeUser(0.3, 20);
205	s->Search(h, .15, "", 0.01);
206	c->Update();
207	TH1* b = s->Background(h);
208	b->SetFillColor(4);
209	b->SetFillStyle(3001);
210	b->Draw("same");
211	#endif
212
213	return kTRUE;
214	}
215
216	ClassDef(DrawESD,0);
217
218	};
219
220	//____________________________________________________________________
221	//
222	// EOF
223	//