[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:
  //__________________________________________________________________
  DrawESD(Int_t n=1000, 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->Sumw2();
    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 
    if (mult/fCorr > 0.01) fMult->Fill(mult/fCorr);
    return kTRUE;
  }
  //__________________________________________________________________
  TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
  {
    TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
    hist->GetXaxis()->SetRangeUser(0, 4);
    hist->Fit(l, "0Q", "", min, max);
    Double_t mpv   = l->GetParameter(1);
    Double_t empv  = l->GetParError(1);
    Double_t sigma = l->GetParameter(2);
    l->SetRange(mpv-empv, mpv+3*sigma);
    hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
    std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
	      << " MPV: " << l->GetParameter(1) 
	      << " +/- "  << l->GetParError(1) 
	      << " Var: " << l->GetParameter(2) 
	      << " +/- "  << l->GetParError(2)
	      << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF() 
	      << std::endl;
    mpv   = l->GetParameter(1);
    sigma = l->GetParameter(2);
    min   = mpv - sigma * 2; // (n==1 ? 3 : 2);
    max   = mpv + sigma * 3;
    // l->SetRange(min, max);
    l->Draw("same");
    return  l;
  }
  //__________________________________________________________________
  void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
  {
    hist->GetXaxis()->SetRangeUser(min, 4);
    mean = hist->GetMean();
    var  = hist->GetRMS();
  }

  //__________________________________________________________________
  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,4);
    fMult->Scale(1. / fMult->GetEntries());
    fMult->SetStats(kFALSE);
    fMult->SetFillColor(2);
    fMult->SetFillStyle(3001);
    fMult->Draw("hist e");

    Double_t mean, rms;
    MaxInRange(fMult, 0.2, mean, rms);
    Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
    Double_t x2   = mean+rms/2; // 1.3; // 1.7; // fCorr;
    TF1*     l1   = FitPeak(1, fMult, x1, x2);
    x2            = TMath::Max(mean+rms/2, x2);
    MaxInRange(fMult, x2, mean, rms);
    Double_t x3   = mean + rms;
    TF1*     l2   = FitPeak(2, fMult, x2, x3);
    Double_t diff = l2->GetParameter(1)-l1->GetParameter(1);
    TF1*     f    = new TF1("user", "landau(0)+landau(3)", x1, x3);

    fMult->GetXaxis()->SetRangeUser(0, 4);
    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),
		     l2->GetParameter(0)/10, 
		     l2->GetParameter(1) + diff, 
		     l2->GetParameter(2));
    fMult->Fit(f, "0Q", "", x1, x3);
    fMult->Fit(f, "ME0", "E1", x1, x3);
    fMult->DrawClone("same hist");

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

    TLegend* l = new TLegend(0.6, 0.6, .89, .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");


#if 1
    c = new TCanvas("c2", "Landaus");
    c->SetLogy();
    fMult->DrawClone("axis");
    f->Draw("same");
    Double_t* p1 = f->GetParameters();
    Double_t* p2 = &(p1[3]);
    TF1* ll1 = new TF1("ll1", "landau", 0, 4);
    ll1->SetParameters(p1);
    ll1->SetLineColor(3);
    ll1->Draw("same");
    TF1* ll2 = new TF1("ll2", "landau", 0, 4);
    ll2->SetParameters(p2);
    ll2->SetLineColor(4);
    ll2->Draw("same");

    Double_t y1  = fMult->GetMinimum() * 1.1;
    Double_t y2  = fMult->GetMaximum() * .9;
    Double_t xc1 = p1[1]-3*p1[2];
    Double_t xc2 = p2[1]-2*p2[2];
    Double_t xc3 = p2[1]-2*p2[2]+diff;
    TLine* c1 = new TLine(xc1, y1, xc1, y2);
    c1->Draw("same");
    TLine* c2 = new TLine(xc2, y1, xc2, y2);
    c2->Draw("same");
    TLine* c3 = new TLine(xc3, y1, xc3, y2);
    c3->Draw("same");

    l = new TLegend(0.6, 0.6, .89, .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");
#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:
2b893216	49	//__________________________________________________________________
f48d9b11	50	DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5)
2b893216	51	: fCorr(1) // 0.68377 / 1.1)
	52	{
	53	AddLoad(kESD);
	54	fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
69893a66	55	fMult->Sumw2();
2b893216	56	fMult->SetXTitle("Strip Multiplicity");
	57	}
	58	//__________________________________________________________________
	59	/** Begining of event
	60	@param ev Event number
	61	@return @c false on error */
	62	Bool_t Begin(Int_t ev)
	63	{
	64	return AliFMDInput::Begin(ev);
	65	}
	66	//__________________________________________________________________
	67	Bool_t ProcessESD(UShort_t /* det */,
	68	Char_t /* ring */,
	69	UShort_t /* sec */,
	70	UShort_t /* strip */,
	71	Float_t /* eta */,
	72	Float_t mult)
	73	{
	74	// Cache the energy loss
f48d9b11	75	if (mult/fCorr > 0.01) fMult->Fill(mult/fCorr);
2b893216	76	return kTRUE;
2b893216	77	}
f48d9b11	78	//__________________________________________________________________
	79	TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
	80	{
	81	TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
	82	hist->GetXaxis()->SetRangeUser(0, 4);
	83	hist->Fit(l, "0Q", "", min, max);
	84	Double_t mpv = l->GetParameter(1);
	85	Double_t empv = l->GetParError(1);
	86	Double_t sigma = l->GetParameter(2);
	87	l->SetRange(mpv-empv, mpv+3*sigma);
	88	hist->Fit(l, "EMQ0", "", mpv-3empv, mpv+3sigma);
	89	std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
	90	<< " MPV: " << l->GetParameter(1)
	91	<< " +/- " << l->GetParError(1)
	92	<< " Var: " << l->GetParameter(2)
	93	<< " +/- " << l->GetParError(2)
	94	<< " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
	95	<< std::endl;
	96	mpv = l->GetParameter(1);
	97	sigma = l->GetParameter(2);
	98	min = mpv - sigma * 2; // (n==1 ? 3 : 2);
	99	max = mpv + sigma * 3;
	100	// l->SetRange(min, max);
	101	l->Draw("same");
	102	return l;
	103	}
	104	//__________________________________________________________________
	105	void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
	106	{
	107	hist->GetXaxis()->SetRangeUser(min, 4);
	108	mean = hist->GetMean();
	109	var = hist->GetRMS();
	110	}
	111
2b893216	112	//__________________________________________________________________
	113	Bool_t Finish()
	114	{
	115	gStyle->SetPalette(1);
	116	gStyle->SetOptTitle(0);
	117	gStyle->SetOptFit(1111111);
	118	gStyle->SetCanvasColor(0);
	119	gStyle->SetCanvasBorderSize(0);
	120	gStyle->SetPadColor(0);
	121	gStyle->SetPadBorderSize(0);
	122
	123	TCanvas* c = new TCanvas("c", "C");
	124	c->cd();
	125	c->SetLogy();
f48d9b11	126	fMult->GetXaxis()->SetRangeUser(0,4);
69893a66	127	fMult->Scale(1. / fMult->GetEntries());
2b893216	128	fMult->SetStats(kFALSE);
	129	fMult->SetFillColor(2);
	130	fMult->SetFillStyle(3001);
f48d9b11	131	fMult->Draw("hist e");
	132
	133	Double_t mean, rms;
	134	MaxInRange(fMult, 0.2, mean, rms);
	135	Double_t x1 = mean-rms/2; // .75; // .8; // .65 / fCorr;
	136	Double_t x2 = mean+rms/2; // 1.3; // 1.7; // fCorr;
	137	TF1* l1 = FitPeak(1, fMult, x1, x2);
	138	x2 = TMath::Max(mean+rms/2, x2);
	139	MaxInRange(fMult, x2, mean, rms);
	140	Double_t x3 = mean + rms;
	141	TF1* l2 = FitPeak(2, fMult, x2, x3);
	142	Double_t diff = l2->GetParameter(1)-l1->GetParameter(1);
	143	TF1* f = new TF1("user", "landau(0)+landau(3)", x1, x3);
2b893216	144
2b893216	145	fMult->GetXaxis()->SetRangeUser(0, 4);
2b893216	146	f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
	147	"A_{2}", "Mpv_{2}", "#sigma_{2}",
	148	"A_{3}", "Mpv_{3}", "#sigma_{3}");
	149	f->SetParameters(l1->GetParameter(0),
	150	l1->GetParameter(1),
	151	l1->GetParameter(2),
	152	l2->GetParameter(0),
	153	l2->GetParameter(1),
	154	l2->GetParameter(2),
f48d9b11	155	l2->GetParameter(0)/10,
	156	l2->GetParameter(1) + diff,
	157	l2->GetParameter(2));
	158	fMult->Fit(f, "0Q", "", x1, x3);
	159	fMult->Fit(f, "ME0", "E1", x1, x3);
	160	fMult->DrawClone("same hist");
2b893216	161
2b893216	162	l1->SetLineColor(3);
	163	l1->SetLineWidth(2);
	164	l1->SetRange(0, 4);
	165	l1->Draw("same");
2b893216	166	l2->SetLineColor(4);
	167	l2->SetLineWidth(2);
	168	l2->SetRange(0, 4);
	169	l2->Draw("same");
2b893216	170	f->SetLineWidth(2);
	171	f->SetRange(0, 4);
	172	f->Draw("same");
	173
f48d9b11	174	TLegend* l = new TLegend(0.6, 0.6, .89, .89);
2b893216	175	l->AddEntry(l1, "1 particle Landau", "l");
	176	l->AddEntry(l2, "2 particle Landau", "l");
	177	l->AddEntry(f, "1+2 particle Landau", "l");
	178	l->SetFillColor(0);
	179	l->Draw("same");
	180
	181
f48d9b11	182	#if 1
2b893216	183	c = new TCanvas("c2", "Landaus");
	184	c->SetLogy();
	185	fMult->DrawClone("axis");
	186	f->Draw("same");
f48d9b11	187	Double_t* p1 = f->GetParameters();
f48d9b11	188	Double_t* p2 = &(p1[3]);
2b893216	189	TF1* ll1 = new TF1("ll1", "landau", 0, 4);
f48d9b11	190	ll1->SetParameters(p1);
2b893216	191	ll1->SetLineColor(3);
	192	ll1->Draw("same");
	193	TF1* ll2 = new TF1("ll2", "landau", 0, 4);
f48d9b11	194	ll2->SetParameters(p2);
2b893216	195	ll2->SetLineColor(4);
	196	ll2->Draw("same");
	197
	198	Double_t y1 = fMult->GetMinimum() * 1.1;
	199	Double_t y2 = fMult->GetMaximum() * .9;
f48d9b11	200	Double_t xc1 = p1[1]-3*p1[2];
	201	Double_t xc2 = p2[1]-2*p2[2];
	202	Double_t xc3 = p2[1]-2*p2[2]+diff;
2b893216	203	TLine* c1 = new TLine(xc1, y1, xc1, y2);
	204	c1->Draw("same");
	205	TLine* c2 = new TLine(xc2, y1, xc2, y2);
	206	c2->Draw("same");
f48d9b11	207	TLine* c3 = new TLine(xc3, y1, xc3, y2);
f48d9b11	208	c3->Draw("same");
2b893216	209
f48d9b11	210	l = new TLegend(0.6, 0.6, .89, .89);
2b893216	211	l->AddEntry(ll1, "1 particle Landau", "l");
	212	l->AddEntry(ll2, "2 particle Landau", "l");
	213	l->AddEntry(f, "1+2 particle Landau", "l");
	214	l->SetFillColor(0);
	215	l->Draw("same");
2b893216	216	#endif
2b893216	217
f48d9b11	218
2b893216	219	return kTRUE;
	220	}
	221
	222	ClassDef(DrawESD,0);
	223
	224	};
	225
	226	//____________________________________________________________________
	227	//
	228	// EOF
	229	//