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

#include <AliFMDHit.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TParticle.h>
#include <TBrowser.h>
#include <TArrayF.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TLegend.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TROOT.h>
#include <AliFMDInput.h>
#include <AliFMDAnaRing.h>
#include <AliFMDGeometry.h>
#include <AliCDBManager.h>
#include <iostream>

//====================================================================
class AliFMDAnaBgRing : public AliFMDAnaRing
{
public:
  AliFMDAnaBgRing(UShort_t det, Char_t ring, 
		  Float_t deta, Float_t etaMin, Float_t etaMax)
    : AliFMDAnaRing(det, ring, 0, 0, 0),
      fTotal(Form("total_%s", fName), "Total number of hits", 
	     Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax, 
	     fNSeq, 0, 2*TMath::Pi()),
      fPrimary(Form("primary_%s", fName), "Number of hits from primary", 
	       Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax, 
	       fNSeq, 0, 2*TMath::Pi()),
      fRatio(Form("ratio_%s", fName), "Correction factors", 
	       Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax, 
	     fNSeq, 0, 2*TMath::Pi()),
      fInverse(Form("inverse_%s", fName), "Inverse Correction factors", 
	       Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax, 
	       fNSeq, 0, 2*TMath::Pi())
  {
    fTotal.SetDirectory(0);
    fTotal.Sumw2();
    fTotal.SetXTitle("#eta");
    fTotal.SetYTitle("#varphi");
    fTotal.SetZTitle("# of hits");
    fPrimary.SetDirectory(0);
    fPrimary.Sumw2();
    fPrimary.SetXTitle("#eta");
    fPrimary.SetYTitle("#varphi");
    fPrimary.SetZTitle("# of primary");
    fRatio.SetDirectory(0);
    fRatio.Sumw2();
    fRatio.SetXTitle("#eta");
    fRatio.SetYTitle("#varphi");
    fRatio.SetZTitle("# of primary/# of hits");
    fInverse.SetDirectory(0);
    fInverse.Sumw2();
    fInverse.SetXTitle("#eta");
    fInverse.SetYTitle("#varphi");
    fInverse.SetZTitle("# of hits/# of primary");
  }
  void Fill(Float_t phi, Float_t eta, Float_t mult)
  {
    fTotal.Fill(eta, phi);
    if (mult >= 1) fPrimary.Fill(eta,phi);
  }
  void Finish()
  {
    fRatio.Add(&fPrimary);
    fRatio.Divide(&fTotal);
    fInverse.Add(&fTotal);
    fInverse.Divide(&fPrimary);
  }
  void Browse(TBrowser* b) 
  {
    // AliFMDAnaRing::Browse(b);
    b->Add(&fTotal);
    b->Add(&fPrimary);
    b->Add(&fRatio);
    b->Add(&fInverse);
  }
  TH2D fTotal;
  TH2D fPrimary;
  TH2D fRatio;
  TH2D fInverse;
};

//====================================================================
class AliFMDAnaBg : public AliFMDInput
{
public:
  //__________________________________________________________________
  /** Constructor */
  AliFMDAnaBg(Float_t deta=.1) 
    : fEtaMin(-3.7), 
      fEtaMax(5.3),
      fTotal("total", "Total number of hits", 
	     Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax, 
	     40, 0, 2*TMath::Pi()),
      fPrimary("primary", "Number of hits from primary", 
	       Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax, 
	       40, 0, 2*TMath::Pi()),
      fRatio("ratio", "Correction factors", 
	     Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax, 
	     40, 0, 2*TMath::Pi())
  {
    fTotal.SetDirectory(0);
    fTotal.Sumw2();
    fTotal.SetXTitle("#eta");
    fTotal.SetYTitle("#varphi");
    fTotal.SetZTitle("# of hits");
    fPrimary.SetDirectory(0);
    fPrimary.Sumw2();
    fPrimary.SetXTitle("#eta");
    fPrimary.SetYTitle("#varphi");
    fPrimary.SetZTitle("# of primary");
    fRatio.SetDirectory(0);
    fRatio.Sumw2();
    fRatio.SetXTitle("#eta");
    fRatio.SetYTitle("#varphi");
    fRatio.SetZTitle("# of primary/# of hits");

    AliCDBManager::Instance()->SetRun(0);

    AddLoad(kHits);
    AddLoad(kKinematics);
    AddLoad(kGeometry);
    for (int i = 0; i < 5; i++) { 
      UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
      Char_t   rng = (i == 0 ? 'I' :  (i % 2 == 1 ? 'I' : 'O'));
      fRing[i] = new AliFMDAnaBgRing(det,rng,deta,EtaMin(i)-.2,EtaMax(i)+.12);
    }
  }
  //__________________________________________________________________
  Int_t FindRing(UShort_t det, Char_t ring) const
  { 
    Int_t idx = -1;
    switch (det) { 
    case 1: idx = 0; break;
    case 2: idx = 1 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
    case 3: idx = 3 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
    }
    return idx;
  }  
  //__________________________________________________________________
  /** Initialize */
  Bool_t Init()
  {
    Bool_t ret = AliFMDInput::Init();
    for (size_t i = 0; i < 5; i++) fRing[i]->Init();
    AliFMDGeometry::Instance()->Init();
    AliFMDGeometry::Instance()->InitTransformations();
    return ret;
  }
  //__________________________________________________________________
  /** Called at the beginning of an event */
  Bool_t Begin(Int_t ev) 
  {
    for (size_t i = 0; i < 5; i++) fRing[i]->Begin();
    return AliFMDInput::Begin(ev);
  }
  //__________________________________________________________________
  /** Called at the beginning of an event */
  Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
  {
    Float_t          mult = p->IsPrimary() ? 1 : 0;
    UShort_t         det  = hit->Detector();
    Char_t           rng  = hit->Ring();
    UShort_t         sec  = hit->Sector();
    UShort_t         str  = hit->Strip();
    AliFMDAnaBgRing* frn  = fRing[FindRing(det, rng)];
    Double_t x, y, z;
    AliFMDGeometry::Instance()->Detector2XYZ(det,rng,sec,str,x,y,z);
    Double_t         phi  = TMath::ATan2(y,x);
    Double_t         r    = TMath::Sqrt(x * x + y * y);
    Double_t         the  = TMath::ATan2(r, z);
    Double_t         eta  = -TMath::Log(TMath::Tan(the/2));
    if (phi < 0) phi += 2*TMath::Pi();
    frn->Fill(phi, eta, mult);
    Fill(phi, eta, mult);
    return kTRUE;
  }

  //__________________________________________________________________
  /** Fill in one strip */
  void Fill(Float_t phi, Float_t eta, Float_t mult)
  {
    fTotal.Fill(eta, phi);
    if (mult >= 1) fPrimary.Fill(eta,phi);
  }
  //__________________________________________________________________
  /** Called at the end of an event */
  Bool_t End() 
  {
    for (UShort_t i = 0; i < 5; i++) fRing[i]->End();
    return AliFMDInput::End();
  }
  //__________________________________________________________________
  /** called at the end of a run */
  Bool_t Finish()
  {
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetOptStat(0);

    fRatio.Add(&fPrimary);
    fRatio.Divide(&fTotal);

    for (UShort_t i = 0; i < 5; i++) fRing[i]->Finish();

    TCanvas* c = new TCanvas("c1", "Result", 800, 600);
    c->cd();
    c->SetGridy();
    c->SetTopMargin(0.02);
    c->SetRightMargin(0.02);
    c->SetFillColor(0);
    c->SetBorderSize(0);
    c->SetBorderMode(0);
    c->SetHighLightColor(0);
    c->Divide(5,1,0,0);
    for (UShort_t i = 0; i < 5; i++) { 
      c->cd(i+1);
      fRing[i]->fRatio.Draw("colz");
    }

    c = new TCanvas("c2", "All", 800, 600);
    c->cd();
    c->SetGridy();
    c->SetTopMargin(0.02);
    c->SetRightMargin(0.02);
    c->SetFillColor(0);
    c->SetBorderSize(0);
    c->SetBorderMode(0);
    c->SetHighLightColor(0);
    c->Divide(5,1,0,0);
    fRatio.Draw("colz");

    // gROOT->GetListOfSpecials()->Add(this);
    return kTRUE;
  }
  //__________________________________________________________________
  Bool_t IsFolder() const { return kTRUE; }
  //__________________________________________________________________
  void Browse(TBrowser* b) 
  {
    for (size_t i = 0; i < 5; i++) 
      b->Add(fRing[i], fRing[i]->Name());
  }
  //__________________________________________________________________
  /** Get min pseudo rapidity of ring idx */
  Float_t EtaMin(Int_t idx) const { 
    switch (idx) { 
    case 0: return  3.61728; break; // FMD1I
    case 1: return  2.18122; break; // FMD2I
    case 2: return  1.79821; break; // FMD2O
    case 3: return -3.39913; break; // FMD3I
    case 4: return -2.28925; break; // FMD3O
    }
    return 0;
  }
  //__________________________________________________________________
  /** Get max pseudo rapidity of ring idx */
  Float_t EtaMax(Int_t idx) const { 
    switch (idx) { 
    case 0: return  5.02643; break; // FMD1I
    case 1: return  3.57899; break; // FMD2I
    case 2: return  2.39085; break; // FMD2O
    case 3: return -2.00644; break; // FMD3I
    case 4: return -1.70080; break; // FMD3O
    }
    return 0;
  }
  /** Flow rings */
  AliFMDAnaBgRing* fRing[5];
  Float_t fEtaMin;
  Float_t fEtaMax;
  TH2D fTotal;
  TH2D fPrimary;
  TH2D fRatio;
  ClassDef(AliFMDAnaBg,1) 
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
9b98d361	1	#include <AliFMDHit.h>
	2	#include <TH1D.h>
	3	#include <TH2D.h>
	4	#include <TParticle.h>
	5	#include <TBrowser.h>
	6	#include <TArrayF.h>
	7	#include <TCanvas.h>
	8	#include <TStyle.h>
	9	#include <TLegend.h>
	10	#include <TProfile.h>
	11	#include <TProfile2D.h>
	12	#include <TROOT.h>
	13	#include <AliFMDInput.h>
	14	#include <AliFMDAnaRing.h>
	15	#include <AliFMDGeometry.h>
	16	#include <AliCDBManager.h>
	17	#include <iostream>
	18
	19	//====================================================================
	20	class AliFMDAnaBgRing : public AliFMDAnaRing
	21	{
	22	public:
	23	AliFMDAnaBgRing(UShort_t det, Char_t ring,
	24	Float_t deta, Float_t etaMin, Float_t etaMax)
	25	: AliFMDAnaRing(det, ring, 0, 0, 0),
	26	fTotal(Form("total_%s", fName), "Total number of hits",
	27	Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax,
	28	fNSeq, 0, 2*TMath::Pi()),
	29	fPrimary(Form("primary_%s", fName), "Number of hits from primary",
	30	Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax,
	31	fNSeq, 0, 2*TMath::Pi()),
	32	fRatio(Form("ratio_%s", fName), "Correction factors",
	33	Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax,
	34	fNSeq, 0, 2*TMath::Pi()),
	35	fInverse(Form("inverse_%s", fName), "Inverse Correction factors",
	36	Int_t((etaMax-etaMin)/deta+.5), etaMin, etaMax,
	37	fNSeq, 0, 2*TMath::Pi())
	38	{
	39	fTotal.SetDirectory(0);
	40	fTotal.Sumw2();
	41	fTotal.SetXTitle("#eta");
	42	fTotal.SetYTitle("#varphi");
	43	fTotal.SetZTitle("# of hits");
	44	fPrimary.SetDirectory(0);
	45	fPrimary.Sumw2();
	46	fPrimary.SetXTitle("#eta");
	47	fPrimary.SetYTitle("#varphi");
	48	fPrimary.SetZTitle("# of primary");
	49	fRatio.SetDirectory(0);
	50	fRatio.Sumw2();
	51	fRatio.SetXTitle("#eta");
	52	fRatio.SetYTitle("#varphi");
	53	fRatio.SetZTitle("# of primary/# of hits");
	54	fInverse.SetDirectory(0);
	55	fInverse.Sumw2();
	56	fInverse.SetXTitle("#eta");
	57	fInverse.SetYTitle("#varphi");
	58	fInverse.SetZTitle("# of hits/# of primary");
	59	}
	60	void Fill(Float_t phi, Float_t eta, Float_t mult)
	61	{
	62	fTotal.Fill(eta, phi);
	63	if (mult >= 1) fPrimary.Fill(eta,phi);
	64	}
65	void Finish()
66	{
67	fRatio.Add(&fPrimary);
68	fRatio.Divide(&fTotal);
69	fInverse.Add(&fTotal);
70	fInverse.Divide(&fPrimary);
71	}
72	void Browse(TBrowser* b)
73	{
74	// AliFMDAnaRing::Browse(b);
75	b->Add(&fTotal);
76	b->Add(&fPrimary);
77	b->Add(&fRatio);
78	b->Add(&fInverse);
79	}
80	TH2D fTotal;
81	TH2D fPrimary;
82	TH2D fRatio;
83	TH2D fInverse;
84	};
85
86	//====================================================================
87	class AliFMDAnaBg : public AliFMDInput
88	{
89	public:
90	//__________________________________________________________________
91	/** Constructor */
92	AliFMDAnaBg(Float_t deta=.1)
93	: fEtaMin(-3.7),
94	fEtaMax(5.3),
95	fTotal("total", "Total number of hits",
96	Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax,
97	40, 0, 2*TMath::Pi()),
98	fPrimary("primary", "Number of hits from primary",
99	Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax,
100	40, 0, 2*TMath::Pi()),
101	fRatio("ratio", "Correction factors",
102	Int_t((fEtaMax-fEtaMin)/deta+.5), fEtaMin, fEtaMax,
103	40, 0, 2*TMath::Pi())
104	{
105	fTotal.SetDirectory(0);
106	fTotal.Sumw2();
107	fTotal.SetXTitle("#eta");
108	fTotal.SetYTitle("#varphi");
109	fTotal.SetZTitle("# of hits");
110	fPrimary.SetDirectory(0);
111	fPrimary.Sumw2();
112	fPrimary.SetXTitle("#eta");
113	fPrimary.SetYTitle("#varphi");
114	fPrimary.SetZTitle("# of primary");
115	fRatio.SetDirectory(0);
116	fRatio.Sumw2();
117	fRatio.SetXTitle("#eta");
118	fRatio.SetYTitle("#varphi");
119	fRatio.SetZTitle("# of primary/# of hits");
120
121	AliCDBManager::Instance()->SetRun(0);
122
123	AddLoad(kHits);
124	AddLoad(kKinematics);
125	AddLoad(kGeometry);
126	for (int i = 0; i < 5; i++) {
127	UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
128	Char_t rng = (i == 0 ? 'I' : (i % 2 == 1 ? 'I' : 'O'));
129	fRing[i] = new AliFMDAnaBgRing(det,rng,deta,EtaMin(i)-.2,EtaMax(i)+.12);
130	}
131	}
132	//__________________________________________________________________
133	Int_t FindRing(UShort_t det, Char_t ring) const
134	{
135	Int_t idx = -1;
136	switch (det) {
137	case 1: idx = 0; break;
138	case 2: idx = 1 + (ring == 'O' \|\| ring == 'o' ? 1 : 0); break;
139	case 3: idx = 3 + (ring == 'O' \|\| ring == 'o' ? 1 : 0); break;
140	}
141	return idx;
142	}
143	//__________________________________________________________________
144	/** Initialize */
145	Bool_t Init()
146	{
147	Bool_t ret = AliFMDInput::Init();
148	for (size_t i = 0; i < 5; i++) fRing[i]->Init();
149	AliFMDGeometry::Instance()->Init();
150	AliFMDGeometry::Instance()->InitTransformations();
151	return ret;
152	}
153	//__________________________________________________________________
154	/** Called at the beginning of an event */
155	Bool_t Begin(Int_t ev)
156	{
157	for (size_t i = 0; i < 5; i++) fRing[i]->Begin();
158	return AliFMDInput::Begin(ev);
159	}
160	//__________________________________________________________________
161	/** Called at the beginning of an event */
162	Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
163	{
164	Float_t mult = p->IsPrimary() ? 1 : 0;
165	UShort_t det = hit->Detector();
166	Char_t rng = hit->Ring();
167	UShort_t sec = hit->Sector();
168	UShort_t str = hit->Strip();
169	AliFMDAnaBgRing* frn = fRing[FindRing(det, rng)];
170	Double_t x, y, z;
171	AliFMDGeometry::Instance()->Detector2XYZ(det,rng,sec,str,x,y,z);
172	Double_t phi = TMath::ATan2(y,x);
173	Double_t r = TMath::Sqrt(x * x + y * y);
174	Double_t the = TMath::ATan2(r, z);
175	Double_t eta = -TMath::Log(TMath::Tan(the/2));
176	if (phi < 0) phi += 2*TMath::Pi();
177	frn->Fill(phi, eta, mult);
178	Fill(phi, eta, mult);
179	return kTRUE;
180	}
181
182	//__________________________________________________________________
183	/** Fill in one strip */
184	void Fill(Float_t phi, Float_t eta, Float_t mult)
185	{
186	fTotal.Fill(eta, phi);
187	if (mult >= 1) fPrimary.Fill(eta,phi);
188	}
189	//__________________________________________________________________
190	/** Called at the end of an event */
191	Bool_t End()
192	{
193	for (UShort_t i = 0; i < 5; i++) fRing[i]->End();
194	return AliFMDInput::End();
195	}
196	//__________________________________________________________________
197	/** called at the end of a run */
198	Bool_t Finish()
199	{
200	gStyle->SetPalette(1);
201	gStyle->SetOptTitle(0);
202	gStyle->SetOptStat(0);
203
204	fRatio.Add(&fPrimary);
205	fRatio.Divide(&fTotal);
206
207	for (UShort_t i = 0; i < 5; i++) fRing[i]->Finish();
208
209	TCanvas* c = new TCanvas("c1", "Result", 800, 600);
210	c->cd();
211	c->SetGridy();
212	c->SetTopMargin(0.02);
213	c->SetRightMargin(0.02);
214	c->SetFillColor(0);
215	c->SetBorderSize(0);
216	c->SetBorderMode(0);
217	c->SetHighLightColor(0);
218	c->Divide(5,1,0,0);
219	for (UShort_t i = 0; i < 5; i++) {
220	c->cd(i+1);
221	fRing[i]->fRatio.Draw("colz");
222	}
223
224	c = new TCanvas("c2", "All", 800, 600);
225	c->cd();
226	c->SetGridy();
227	c->SetTopMargin(0.02);
228	c->SetRightMargin(0.02);
229	c->SetFillColor(0);
230	c->SetBorderSize(0);
231	c->SetBorderMode(0);
232	c->SetHighLightColor(0);
233	c->Divide(5,1,0,0);
234	fRatio.Draw("colz");
235
236	// gROOT->GetListOfSpecials()->Add(this);
237	return kTRUE;
238	}
239	//__________________________________________________________________
240	Bool_t IsFolder() const { return kTRUE; }
241	//__________________________________________________________________
242	void Browse(TBrowser* b)
243	{
244	for (size_t i = 0; i < 5; i++)
245	b->Add(fRing[i], fRing[i]->Name());
246	}
247	//__________________________________________________________________
248	/** Get min pseudo rapidity of ring idx */
249	Float_t EtaMin(Int_t idx) const {
250	switch (idx) {
251	case 0: return 3.61728; break; // FMD1I
252	case 1: return 2.18122; break; // FMD2I
253	case 2: return 1.79821; break; // FMD2O
254	case 3: return -3.39913; break; // FMD3I
255	case 4: return -2.28925; break; // FMD3O
256	}
257	return 0;
258	}
259	//__________________________________________________________________
260	/** Get max pseudo rapidity of ring idx */
261	Float_t EtaMax(Int_t idx) const {
262	switch (idx) {
263	case 0: return 5.02643; break; // FMD1I
264	case 1: return 3.57899; break; // FMD2I
265	case 2: return 2.39085; break; // FMD2O
266	case 3: return -2.00644; break; // FMD3I
267	case 4: return -1.70080; break; // FMD3O
268	}
269	return 0;
270	}
271	/** Flow rings */
272	AliFMDAnaBgRing* fRing[5];
273	Float_t fEtaMin;
274	Float_t fEtaMax;
275	TH2D fTotal;
276	TH2D fPrimary;
277	TH2D fRatio;
278	ClassDef(AliFMDAnaBg,1)
279	};
280
281	//____________________________________________________________________
282	//
283	// EOF
284	//