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

//____________________________________________________________________
//
// $Id: DrawSDigits.C 20907 2007-09-25 08:44:03Z cholm $
//
// 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 <TH2D.h>
#include <TProfile2D.h>
#include <TCanvas.h>
#include <TMath.h>
#include <AliFMDSDigit.h>
#include <AliFMDInput.h>
#include <AliFMDEdepMap.h>
#include <AliFMDGeometry.h>
#include <iostream>
#include <TStyle.h>
#include <TLatex.h>
#include <TArrayF.h>
#include <AliLog.h>

/** @class DrawSDigits
    @brief Draw hit energy loss versus digit ADC
    @code 
    Root> .L Compile.C
    Root> Compile("DrawSDigits.C")
    Root> DrawSDigits c
    Root> c.Run();
    @endcode
    @ingroup FMD_script
 */
class DrawSDigits : public AliFMDInput
{
private:
  TH1D* fAdc; // Histogram 
  TProfile2D* fPrimRatio[5]; // Histograms
public:
  void Idx2Det(UShort_t idx, UShort_t& d, Char_t& r) const
  {
    d = 0;
    r = '\0';
    switch (idx) { 
    case 0: d = 1; r = 'I'; break;
    case 1: d = 2; r = 'I'; break;
    case 2: d = 2; r = 'O'; break;
    case 3: d = 3; r = 'I'; break;
    case 4: d = 3; r = 'O'; break;
    }
  }
  Short_t Det2Idx(UShort_t d, Char_t r) const
  { 
    Short_t idx = -1;
    switch (d) { 
    case 1: idx = 0; break;
    case 2: idx = 1; break;
    case 3: idx = 3; break;
    }
    return (idx + ((r == 'O' || r == 'o') ? 1 : 0));
  }
  
    
  //__________________________________________________________________
  DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) 
    : AliFMDInput("galice.root")
  { 
    AddLoad(kSDigits);
    AddLoad(kGeometry);
    fAdc = new TH1D("adc", "ADC", m, amin, amax);
    fAdc->SetXTitle("ADC value");
    fAdc->Sumw2();

    Int_t   eN   = 50;
    Float_t eMin = -4;
    Float_t eMax = 6;
    Int_t   pN   = 40;
    Float_t pMin = -.5;
    Float_t pMax = 359.5;
    for (UShort_t i = 0; i < 5; i++) { 
      UShort_t d;
      Char_t   r;
      Idx2Det(i, d, r);
      fPrimRatio[i] = new TProfile2D(Form("primRatio_fmd%d%c", d, r), 
				     Form("Primary/Total - FMD%d%c", d, r), 
				     eN, eMin, eMax, pN, pMin, pMax);
      fPrimRatio[i]->SetXTitle("#eta");
      fPrimRatio[i]->SetYTitle("#varphi");
      fPrimRatio[i]->SetZTitle("M_{ch,primary}/M_{ch,total}");
      fPrimRatio[i]->GetXaxis()->SetTitleFont(132);
      fPrimRatio[i]->GetYaxis()->SetTitleFont(132);
      fPrimRatio[i]->GetZaxis()->SetTitleFont(132);
      fPrimRatio[i]->GetXaxis()->SetLabelFont(132);
      fPrimRatio[i]->GetYaxis()->SetLabelFont(132);
      fPrimRatio[i]->GetZaxis()->SetLabelFont(132);
      fPrimRatio[i]->GetXaxis()->SetTitleSize(.06);
      fPrimRatio[i]->GetYaxis()->SetTitleSize(.06);
      fPrimRatio[i]->GetZaxis()->SetTitleSize(.06);
      fPrimRatio[i]->GetXaxis()->SetLabelSize(.06);
      fPrimRatio[i]->GetYaxis()->SetLabelSize(.06);
      fPrimRatio[i]->GetZaxis()->SetLabelSize(.06);
      fPrimRatio[i]->GetXaxis()->SetNdivisions(10);
      fPrimRatio[i]->GetYaxis()->SetNdivisions(10);
      fPrimRatio[i]->GetZaxis()->SetNdivisions(10);
    }
  }
  Bool_t Init() 
  {
    Bool_t ret = AliFMDInput::Init();
    AliFMDGeometry::Instance()->Init();
    AliFMDGeometry::Instance()->InitTransformations();
    return ret;
  }
    
  //__________________________________________________________________
  Bool_t ProcessSDigit(AliFMDSDigit* digit)
  {
    if (!digit) return kTRUE;
    fAdc->Fill(digit->Counts());
    digit->Print("lp");
    if (digit->NParticles() == 0) return kTRUE;
    

    Short_t primIdx = Det2Idx(digit->Detector(), digit->Ring());
    if (primIdx < 0) return kTRUE;
    
    AliFMDGeometry* geom = AliFMDGeometry::Instance();
    Double_t x, y, z;
    geom->Detector2XYZ(digit->Detector(), 
		       digit->Ring(), 
		       digit->Sector(), 
		       digit->Strip(), 
		       x, y, z);
    // We should get the primary vertex and do 
    // z     += fCurrentVertex;
    Double_t phi   =  TMath::ATan2(y, x) * 180. / TMath::Pi();
    Double_t r     =  TMath::Sqrt(y * y + x * x);
    Double_t theta =  TMath::ATan2(r, z);
    Double_t eta   = -TMath::Log(TMath::Tan(theta / 2));
    Double_t ratio = digit->NPrimaries() / digit->NParticles();
    if (phi < 0) phi += 360;
    fPrimRatio[primIdx]->Fill(eta, phi, ratio);
    return kTRUE;
  }
  //__________________________________________________________________
  Bool_t Finish()
  {
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderSize(0);
    gStyle->SetPadColor(0);
    gStyle->SetPadBorderSize(0);
    
    TCanvas* c1 = new TCanvas("fmdSdigitSpectra", "FMD SDigit spectra");
    fAdc->SetStats(kFALSE);
    if (fAdc->GetEntries() != 0) 
      fAdc->Scale(1. / fAdc->GetEntries());
    fAdc->Draw();
    c1->Modified();
    c1->Update();
    c1->cd();

    TCanvas* c2 = new TCanvas("fmdSdigitPrims", "FMD SDigit # prim/# ntotal",
			      800, 800);
    c2->SetTopMargin(0);
    c2->SetRightMargin(0);
    c2->Divide(2, 3, 0, 0);
    for (UShort_t i = 0; i < 5; i++) { 
      Int_t        np = i+1+(i == 0 ? 0 : 1);
      TVirtualPad* p   = c2->cd(np);
      p->SetFillColor(0);
      p->SetTopMargin((np <= 2) ? 0.05 : 0);
      p->SetLeftMargin((np % 2 == 1) ? 0.20 : 0);
      p->SetRightMargin((np % 2 == 1) ? 0 : 0.20);
      p->SetBottomMargin((np >= 5) ? 0.15 : 0);
      fPrimRatio[i]->SetStats(kFALSE);
      fPrimRatio[i]->Draw(Form("COL%c", (np % 2 == 1) ? ' ' : 'Z'));
      UShort_t d;
      Char_t   r;
      Idx2Det(i, d, r);
      TLatex* text = new TLatex(0, 180, Form("FMD%d%c", d, r));
      text->SetTextFont(132);
      text->SetTextAlign(22);
      text->SetTextSize(0.08);
      text->Draw();
    }
    c2->Modified();
    c2->Update();
    c2->cd();

    return kTRUE;
  }

  ClassDef(DrawSDigits,0);
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
83ad576a	1	//____________________________________________________________________
	2	//
	3	// $Id: DrawSDigits.C 20907 2007-09-25 08:44:03Z cholm $
	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 <TH2D.h>
	16	#include <TProfile2D.h>
	17	#include <TCanvas.h>
	18	#include <TMath.h>
	19	#include <AliFMDSDigit.h>
	20	#include <AliFMDInput.h>
	21	#include <AliFMDEdepMap.h>
	22	#include <AliFMDGeometry.h>
	23	#include <iostream>
	24	#include <TStyle.h>
	25	#include <TLatex.h>
	26	#include <TArrayF.h>
	27	#include <AliLog.h>
	28
	29	/** @class DrawSDigits
	30	@brief Draw hit energy loss versus digit ADC
	31	@code
	32	Root> .L Compile.C
	33	Root> Compile("DrawSDigits.C")
	34	Root> DrawSDigits c
	35	Root> c.Run();
	36	@endcode
	37	@ingroup FMD_script
	38	*/
	39	class DrawSDigits : public AliFMDInput
	40	{
	41	private:
	42	TH1D* fAdc; // Histogram
	43	TProfile2D* fPrimRatio[5]; // Histograms
	44	public:
	45	void Idx2Det(UShort_t idx, UShort_t& d, Char_t& r) const
	46	{
	47	d = 0;
	48	r = '\0';
	49	switch (idx) {
	50	case 0: d = 1; r = 'I'; break;
	51	case 1: d = 2; r = 'I'; break;
	52	case 2: d = 2; r = 'O'; break;
	53	case 3: d = 3; r = 'I'; break;
	54	case 4: d = 3; r = 'O'; break;
	55	}
	56	}
	57	Short_t Det2Idx(UShort_t d, Char_t r) const
	58	{
	59	Short_t idx = -1;
	60	switch (d) {
	61	case 1: idx = 0; break;
	62	case 2: idx = 1; break;
	63	case 3: idx = 3; break;
	64	}
65	return (idx + ((r == 'O' \|\| r == 'o') ? 1 : 0));
66	}
67
68
69	//__________________________________________________________________
70	DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5)
faf80567	71	: AliFMDInput("galice.root")
83ad576a	72	{
	73	AddLoad(kSDigits);
	74	AddLoad(kGeometry);
	75	fAdc = new TH1D("adc", "ADC", m, amin, amax);
	76	fAdc->SetXTitle("ADC value");
	77	fAdc->Sumw2();
	78
	79	Int_t eN = 50;
	80	Float_t eMin = -4;
	81	Float_t eMax = 6;
	82	Int_t pN = 40;
	83	Float_t pMin = -.5;
	84	Float_t pMax = 359.5;
	85	for (UShort_t i = 0; i < 5; i++) {
	86	UShort_t d;
	87	Char_t r;
	88	Idx2Det(i, d, r);
	89	fPrimRatio[i] = new TProfile2D(Form("primRatio_fmd%d%c", d, r),
	90	Form("Primary/Total - FMD%d%c", d, r),
	91	eN, eMin, eMax, pN, pMin, pMax);
	92	fPrimRatio[i]->SetXTitle("#eta");
	93	fPrimRatio[i]->SetYTitle("#varphi");
	94	fPrimRatio[i]->SetZTitle("M_{ch,primary}/M_{ch,total}");
	95	fPrimRatio[i]->GetXaxis()->SetTitleFont(132);
	96	fPrimRatio[i]->GetYaxis()->SetTitleFont(132);
	97	fPrimRatio[i]->GetZaxis()->SetTitleFont(132);
	98	fPrimRatio[i]->GetXaxis()->SetLabelFont(132);
	99	fPrimRatio[i]->GetYaxis()->SetLabelFont(132);
	100	fPrimRatio[i]->GetZaxis()->SetLabelFont(132);
	101	fPrimRatio[i]->GetXaxis()->SetTitleSize(.06);
	102	fPrimRatio[i]->GetYaxis()->SetTitleSize(.06);
	103	fPrimRatio[i]->GetZaxis()->SetTitleSize(.06);
	104	fPrimRatio[i]->GetXaxis()->SetLabelSize(.06);
	105	fPrimRatio[i]->GetYaxis()->SetLabelSize(.06);
	106	fPrimRatio[i]->GetZaxis()->SetLabelSize(.06);
	107	fPrimRatio[i]->GetXaxis()->SetNdivisions(10);
	108	fPrimRatio[i]->GetYaxis()->SetNdivisions(10);
	109	fPrimRatio[i]->GetZaxis()->SetNdivisions(10);
	110	}
	111	}
	112	Bool_t Init()
	113	{
	114	Bool_t ret = AliFMDInput::Init();
	115	AliFMDGeometry::Instance()->Init();
	116	AliFMDGeometry::Instance()->InitTransformations();
	117	return ret;
	118	}
	119
	120	//__________________________________________________________________
	121	Bool_t ProcessSDigit(AliFMDSDigit* digit)
	122	{
	123	if (!digit) return kTRUE;
	124	fAdc->Fill(digit->Counts());
faf80567	125	digit->Print("lp");
83ad576a	126	if (digit->NParticles() == 0) return kTRUE;
	127
	128
	129	Short_t primIdx = Det2Idx(digit->Detector(), digit->Ring());
	130	if (primIdx < 0) return kTRUE;
	131
	132	AliFMDGeometry* geom = AliFMDGeometry::Instance();
	133	Double_t x, y, z;
	134	geom->Detector2XYZ(digit->Detector(),
	135	digit->Ring(),
	136	digit->Sector(),
	137	digit->Strip(),
	138	x, y, z);
	139	// We should get the primary vertex and do
	140	// z += fCurrentVertex;
	141	Double_t phi = TMath::ATan2(y, x) * 180. / TMath::Pi();
	142	Double_t r = TMath::Sqrt(y * y + x * x);
	143	Double_t theta = TMath::ATan2(r, z);
	144	Double_t eta = -TMath::Log(TMath::Tan(theta / 2));
	145	Double_t ratio = digit->NPrimaries() / digit->NParticles();
	146	if (phi < 0) phi += 360;
	147	fPrimRatio[primIdx]->Fill(eta, phi, ratio);
83ad576a	148	return kTRUE;
	149	}
	150	//__________________________________________________________________
	151	Bool_t Finish()
	152	{
	153	gStyle->SetPalette(1);
	154	gStyle->SetOptTitle(0);
	155	gStyle->SetCanvasColor(0);
	156	gStyle->SetCanvasBorderSize(0);
	157	gStyle->SetPadColor(0);
	158	gStyle->SetPadBorderSize(0);
	159
	160	TCanvas* c1 = new TCanvas("fmdSdigitSpectra", "FMD SDigit spectra");
	161	fAdc->SetStats(kFALSE);
	162	if (fAdc->GetEntries() != 0)
	163	fAdc->Scale(1. / fAdc->GetEntries());
	164	fAdc->Draw();
	165	c1->Modified();
	166	c1->Update();
	167	c1->cd();
	168
	169	TCanvas* c2 = new TCanvas("fmdSdigitPrims", "FMD SDigit # prim/# ntotal",
	170	800, 800);
	171	c2->SetTopMargin(0);
	172	c2->SetRightMargin(0);
	173	c2->Divide(2, 3, 0, 0);
	174	for (UShort_t i = 0; i < 5; i++) {
	175	Int_t np = i+1+(i == 0 ? 0 : 1);
	176	TVirtualPad* p = c2->cd(np);
	177	p->SetFillColor(0);
	178	p->SetTopMargin((np <= 2) ? 0.05 : 0);
	179	p->SetLeftMargin((np % 2 == 1) ? 0.20 : 0);
	180	p->SetRightMargin((np % 2 == 1) ? 0 : 0.20);
	181	p->SetBottomMargin((np >= 5) ? 0.15 : 0);
	182	fPrimRatio[i]->SetStats(kFALSE);
	183	fPrimRatio[i]->Draw(Form("COL%c", (np % 2 == 1) ? ' ' : 'Z'));
	184	UShort_t d;
	185	Char_t r;
	186	Idx2Det(i, d, r);
	187	TLatex* text = new TLatex(0, 180, Form("FMD%d%c", d, r));
	188	text->SetTextFont(132);
	189	text->SetTextAlign(22);
	190	text->SetTextSize(0.08);
	191	text->Draw();
	192	}
	193	c2->Modified();
	194	c2->Update();
	195	c2->cd();
	196
	197	return kTRUE;
	198	}
	199
	200	ClassDef(DrawSDigits,0);
	201	};
	202
	203	//____________________________________________________________________
	204	//
	205	// EOF
	206	//