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

#include <TH1D.h>
#include <TH2D.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 "AliFMDAnaESD.h"
#include "AliFMDAnaFlowRing.h"
#include <TFile.h>
#if 1
# include "background.C"
// # include "background2.C"
#else
TH2* bgFMD1i() { return 0; }
TH2* bgFMD2i() { return 0; }
TH2* bgFMD2o() { return 0; }
TH2* bgFMD3o() { return 0; }
TH2* bgFMD3o() { return 0; }
#endif
#if 1
# include "reactionplane.C"
#else
Float_t reactionplane(int ev) { return 0; }
#endif
#include <iostream>


//====================================================================
class AliFMDAnaFlow : public AliFMDAnaESD
{
public:
  //__________________________________________________________________
  /** Constructor */
  AliFMDAnaFlow(Int_t n=2, Bool_t bg=false, Float_t cut0=.68) 
    : fAxis(n, -3.7, 5.3),
      fV2("all", "Full FMD", 2, 1, fAxis)
  {
    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'));
      AddRing(new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)+.1,EtaMax(i)-.1,
			   Bg(bg,i),cut0,Cut1(i)));
    }
  }
  //__________________________________________________________________
  /** Constructor */
  AliFMDAnaFlow(Float_t deta, Bool_t bg=false, Float_t cut0=.68) 
    : fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
      fV2("all", "Full FMD", 2, 1, fAxis)
  {
    for (int i = 0; i < 5; i++) { 
      UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
      Char_t   rng = (i % 2 == 0 ? 'I' : 'O');
      AddRing(new AliFMDAnaFlowRing(det,rng,fAxis,Bg(bg,i),cut0,Cut1(i)));
    }
  }
  //__________________________________________________________________
  /** Initialize */
  Bool_t Init()
  {
    fPhis.Set(51200);
    fEtas.Set(51200);
    fWeights.Set(51200);
    fNObserv = 0;

    Float_t pi2 = 2 * TMath::Pi();

    fPsi  = new TH2D("all_psi", "Event planes", 
		     fAxis.N(), fAxis.Bins(), 40, 0, pi2);
    fPsi->SetXTitle("#eta");
    fPsi->SetYTitle("#Psi");
    fPsi->SetZTitle("Frequency");
    fPsi->SetMarkerStyle(20);

    fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}", 
			 fAxis.N(), fAxis.Bins(), -pi2, pi2);
    fDPsi->SetXTitle("#eta");
    fDPsi->SetYTitle("#Delta#Psi=|#Psi-#Psi_{R}|");

    fDPhi = new TProfile2D("all_dphi","Distribution of phi signals", 
			  fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
    fDPhi->SetXTitle("#eta");
    fDPhi->SetYTitle("#varphi");
    fDPhi->SetZTitle("signal");
    fDPhi->SetMarkerStyle(20);

    return AliFMDAnaESD::Init();
  }
  //__________________________________________________________________
  /** Called at the beginning of an event */
  Bool_t Begin(Int_t ev) 
  {
    fNObserv = 0;
    Float_t psir = reactionplane(fNEvents);
    AliFMDAnaFlowRing* fring = 0;
    for (size_t i = 0; i < 5; i++) {
      if (!fRing[i]) continue;
      fring = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
      fring->fTrue.SetPsi(psir);
    }
    return AliFMDAnaESD::Begin(ev);
  }
  //__________________________________________________________________
  /** Fill in one strip */
  void Fill(Float_t phi, Float_t eta, Float_t mult)
  {
    Int_t imult = Int_t(mult+.33);
    if (fNObserv + imult >= fEtas.fN) { 
      ULong_t n = Int_t(fEtas.fN * 1.5)+imult;
      fEtas.Set(n);
      fPhis.Set(n);
      fWeights.Set(n);
    }
    for (Int_t i = 0; i < imult; i++) { 
      fEtas[fNObserv]    = eta;
      fPhis[fNObserv]    = phi;
      fWeights[fNObserv] = 1;
      fNObserv++;
    }
    fDPhi->Fill(eta, phi);
  }
  //__________________________________________________________________
  /** Called at the end of an event */
  Bool_t End() 
  {
    float       psir    = reactionplane(fNEvents-1);
    // std::cout << "PsiR=" << psir << std::endl;
    fV2.Event(fNObserv, fPhis.fArray, fEtas.fArray, 
	      fWeights.fArray, fWeights.fArray);
    for (UShort_t i = 0; i < fAxis.N(); i++) { 
      AliFMDFlowBin* bin = fV2.GetBin(i);
      if (bin && bin->Counts() > 0) {
	fPsi->Fill(fAxis.BinCenter(i),  bin->Psi() - psir);
	fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
      }
    }
    return AliFMDAnaESD::End();
  }
  //__________________________________________________________________
  /** Make a frame */
  TH1* MakeFrame(const char* xname, const char* yname, 
		 Float_t ymin, Float_t ymax)
  {
    static int id = 0;
    id++;
    TH1* frame = new TH1D(Form("frame%d", id), "Frame", 
			  fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
    frame->SetMinimum(ymin);
    frame->SetMaximum(ymax);
    frame->SetXTitle(xname);
    frame->SetYTitle(yname);
    frame->SetLabelFont(132, "XYZ");
    frame->SetTitleFont(132, "XYZ");
    frame->SetLabelSize(.08, "XYZ");
    frame->SetTitleSize(.08, "XYZ");
    frame->SetTitleOffset(.5, "XYZ");
    frame->SetNdivisions(8, "XYZ");
    frame->Draw();
    return frame;
  }
  //__________________________________________________________________
  /** called at the end of a run */
  Bool_t Finish()
  {
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetOptStat(0);
    AliFMDAnaFlowRing* fring[5];
    for (size_t i = 0; i < 5; i++) 
      fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
    
    TCanvas* c = new TCanvas("c1", "Result", 600, 750);
    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(1,3,0,0);

    TVirtualPad* p = 0;
    TH1*         f = 0;

    p = c->cd(1);
    p->SetGridx();
    p->SetGridy();
    p->SetFillColor(0);
    f = MakeFrame("#eta", "|#Psi_{2}-#Psi_{R}|", -TMath::Pi(), TMath::Pi());
    for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same");
    fDPsi->Draw("same");

    p = c->cd(2);
    p->SetGridx();
    p->SetGridy();
    p->SetFillColor(0);
    f = MakeFrame("#eta", "R_{k}", .7, 1);
    for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
    fV2.Draw("tr:same");

    p = c->cd(3);
    p->SetGridx();
    p->SetGridy();
    p->SetFillColor(0);
    f = MakeFrame("#eta", "v_{2}", 0, .1);
    for (int i = 0; i < 5; i++) { 
      fring[i]->fV2.Draw("th:same");
      fring[i]->fTrue.Draw("th:same");
    }
    fV2.Draw("th:same");


    for (int i = 0; i < 5; i++) { 
      fring[i]->fV2.Print("s");
      fring[i]->fTrue.Print("s");
    }
    fV2.Print("s");
    gROOT->GetListOfSpecials()->Add(this);
    return kTRUE;
  }
  void ToFile(const char* file) { 
    TFile* out = TFile::Open(file, "RECREATE");
    for (int i = 0; i < 5; i++) fRing[i]->Write(fRing[i]->Name());
    fV2.Write();
    if (fPsi) fPsi->Write();
    if (fDPsi) fDPsi->Write();
    if (fDPhi) fDPhi->Write();
    out->Close();
  }
    
    
  //__________________________________________________________________
  /** 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;
  }
  //__________________________________________________________________
  /** Get cut 1 of ring idx */
  Float_t Cut1(Int_t idx) const { 
    switch (idx) { 
    case 0: return 1.33857; break; // FMD1I
    case 1: return 1.30464; break; // FMD2I
    case 2: return 1.27263; break; // FMD2O
    case 3: return 1.28052; break; // FMD3I
    case 4: return 1.23041; break; // FMD3O
    }
    return 0;
  }
  //__________________________________________________________________
  /** Get background correction of ring idx */
  TH2* Bg(Bool_t use, Int_t idx) const { 
    if (!use) return 0;
    switch (idx) { 
    case 0: return bgFMD1i(); break;
    case 1: return bgFMD2i(); break;
    case 2: return bgFMD2o(); break;
    case 3: return bgFMD3i(); break;
    case 4: return bgFMD3o(); break;
    }
    return 0;
  }
  /** Flow axis */
  AliFMDFlowAxis fAxis;
  /** Flow object */
  AliFMDFlowBinned1D fV2;
  /** Cache of phi */
  TArrayD            fPhis;
  /** Cache of eta */
  TArrayD            fEtas;
  /** Cache of weights */
  TArrayD            fWeights;
  /** Number of observables */
  Int_t              fNObserv;
  /** Histogram of Psi */
  TH2D*              fPsi;       // Event planes
  /** Histogram of dPsi */
  TProfile*          fDPsi;   // Mean distance to Psi
  /** Histogram of dPhi */
  TProfile2D*        fDPhi;
  ClassDef(AliFMDAnaFlow,1) 
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
9b98d361	1	#include <TH1D.h>
	2	#include <TH2D.h>
	3	#include <TBrowser.h>
	4	#include <TArrayF.h>
	5	#include <TCanvas.h>
	6	#include <TStyle.h>
	7	#include <TLegend.h>
	8	#include <TProfile.h>
	9	#include <TProfile2D.h>
	10	#include <TROOT.h>
	11	#include "AliFMDAnaESD.h"
	12	#include "AliFMDAnaFlowRing.h"
	13	#include <TFile.h>
	14	#if 1
	15	# include "background.C"
	16	// # include "background2.C"
	17	#else
	18	TH2* bgFMD1i() { return 0; }
	19	TH2* bgFMD2i() { return 0; }
	20	TH2* bgFMD2o() { return 0; }
	21	TH2* bgFMD3o() { return 0; }
	22	TH2* bgFMD3o() { return 0; }
	23	#endif
	24	#if 1
	25	# include "reactionplane.C"
	26	#else
	27	Float_t reactionplane(int ev) { return 0; }
	28	#endif
	29	#include <iostream>
	30
	31
	32	//====================================================================
	33	class AliFMDAnaFlow : public AliFMDAnaESD
	34	{
	35	public:
	36	//__________________________________________________________________
	37	/** Constructor */
	38	AliFMDAnaFlow(Int_t n=2, Bool_t bg=false, Float_t cut0=.68)
	39	: fAxis(n, -3.7, 5.3),
	40	fV2("all", "Full FMD", 2, 1, fAxis)
	41	{
	42	for (int i = 0; i < 5; i++) {
	43	UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
	44	Char_t rng = (i == 0 ? 'I' : (i % 2 == 1 ? 'I' : 'O'));
	45	AddRing(new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)+.1,EtaMax(i)-.1,
	46	Bg(bg,i),cut0,Cut1(i)));
	47	}
	48	}
	49	//__________________________________________________________________
	50	/** Constructor */
	51	AliFMDAnaFlow(Float_t deta, Bool_t bg=false, Float_t cut0=.68)
	52	: fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
	53	fV2("all", "Full FMD", 2, 1, fAxis)
	54	{
	55	for (int i = 0; i < 5; i++) {
	56	UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
	57	Char_t rng = (i % 2 == 0 ? 'I' : 'O');
	58	AddRing(new AliFMDAnaFlowRing(det,rng,fAxis,Bg(bg,i),cut0,Cut1(i)));
	59	}
	60	}
	61	//__________________________________________________________________
	62	/** Initialize */
	63	Bool_t Init()
	64	{
65	fPhis.Set(51200);
66	fEtas.Set(51200);
67	fWeights.Set(51200);
68	fNObserv = 0;
69
70	Float_t pi2 = 2 * TMath::Pi();
71
72	fPsi = new TH2D("all_psi", "Event planes",
73	fAxis.N(), fAxis.Bins(), 40, 0, pi2);
74	fPsi->SetXTitle("#eta");
75	fPsi->SetYTitle("#Psi");
76	fPsi->SetZTitle("Frequency");
77	fPsi->SetMarkerStyle(20);
78
79	fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}",
80	fAxis.N(), fAxis.Bins(), -pi2, pi2);
81	fDPsi->SetXTitle("#eta");
82	fDPsi->SetYTitle("#Delta#Psi=\|#Psi-#Psi_{R}\|");
83
84	fDPhi = new TProfile2D("all_dphi","Distribution of phi signals",
85	fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
86	fDPhi->SetXTitle("#eta");
87	fDPhi->SetYTitle("#varphi");
88	fDPhi->SetZTitle("signal");
89	fDPhi->SetMarkerStyle(20);
90
91	return AliFMDAnaESD::Init();
92	}
93	//__________________________________________________________________
94	/** Called at the beginning of an event */
95	Bool_t Begin(Int_t ev)
96	{
97	fNObserv = 0;
98	Float_t psir = reactionplane(fNEvents);
99	AliFMDAnaFlowRing* fring = 0;
100	for (size_t i = 0; i < 5; i++) {
101	if (!fRing[i]) continue;
102	fring = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
103	fring->fTrue.SetPsi(psir);
104	}
105	return AliFMDAnaESD::Begin(ev);
106	}
107	//__________________________________________________________________
108	/** Fill in one strip */
109	void Fill(Float_t phi, Float_t eta, Float_t mult)
110	{
111	Int_t imult = Int_t(mult+.33);
112	if (fNObserv + imult >= fEtas.fN) {
113	ULong_t n = Int_t(fEtas.fN * 1.5)+imult;
114	fEtas.Set(n);
115	fPhis.Set(n);
116	fWeights.Set(n);
117	}
118	for (Int_t i = 0; i < imult; i++) {
119	fEtas[fNObserv] = eta;
120	fPhis[fNObserv] = phi;
121	fWeights[fNObserv] = 1;
122	fNObserv++;
123	}
124	fDPhi->Fill(eta, phi);
125	}
126	//__________________________________________________________________
127	/** Called at the end of an event */
128	Bool_t End()
129	{
130	float psir = reactionplane(fNEvents-1);
131	// std::cout << "PsiR=" << psir << std::endl;
132	fV2.Event(fNObserv, fPhis.fArray, fEtas.fArray,
133	fWeights.fArray, fWeights.fArray);
134	for (UShort_t i = 0; i < fAxis.N(); i++) {
135	AliFMDFlowBin* bin = fV2.GetBin(i);
136	if (bin && bin->Counts() > 0) {
137	fPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
138	fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
139	}
140	}
141	return AliFMDAnaESD::End();
142	}
143	//__________________________________________________________________
144	/** Make a frame */
145	TH1* MakeFrame(const char* xname, const char* yname,
146	Float_t ymin, Float_t ymax)
147	{
148	static int id = 0;
149	id++;
150	TH1* frame = new TH1D(Form("frame%d", id), "Frame",
151	fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
152	frame->SetMinimum(ymin);
153	frame->SetMaximum(ymax);
154	frame->SetXTitle(xname);
155	frame->SetYTitle(yname);
156	frame->SetLabelFont(132, "XYZ");
157	frame->SetTitleFont(132, "XYZ");
158	frame->SetLabelSize(.08, "XYZ");
159	frame->SetTitleSize(.08, "XYZ");
160	frame->SetTitleOffset(.5, "XYZ");
161	frame->SetNdivisions(8, "XYZ");
162	frame->Draw();
163	return frame;
164	}
165	//__________________________________________________________________
166	/** called at the end of a run */
167	Bool_t Finish()
168	{
169	gStyle->SetPalette(1);
170	gStyle->SetOptTitle(0);
171	gStyle->SetOptStat(0);
172	AliFMDAnaFlowRing* fring[5];
173	for (size_t i = 0; i < 5; i++)
174	fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
175
176	TCanvas* c = new TCanvas("c1", "Result", 600, 750);
177	c->cd();
178	c->SetGridy();
179	c->SetTopMargin(0.02);
180	c->SetRightMargin(0.02);
181	c->SetFillColor(0);
182	c->SetBorderSize(0);
183	c->SetBorderMode(0);
184	c->SetHighLightColor(0);
185	c->Divide(1,3,0,0);
186
187	TVirtualPad* p = 0;
188	TH1* f = 0;
189
190	p = c->cd(1);
191	p->SetGridx();
192	p->SetGridy();
193	p->SetFillColor(0);
194	f = MakeFrame("#eta", "\|#Psi_{2}-#Psi_{R}\|", -TMath::Pi(), TMath::Pi());
195	for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same");
196	fDPsi->Draw("same");
197
198	p = c->cd(2);
199	p->SetGridx();
200	p->SetGridy();
201	p->SetFillColor(0);
202	f = MakeFrame("#eta", "R_{k}", .7, 1);
203	for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
204	fV2.Draw("tr:same");
205
206	p = c->cd(3);
207	p->SetGridx();
208	p->SetGridy();
209	p->SetFillColor(0);
210	f = MakeFrame("#eta", "v_{2}", 0, .1);
211	for (int i = 0; i < 5; i++) {
212	fring[i]->fV2.Draw("th:same");
213	fring[i]->fTrue.Draw("th:same");
214	}
215	fV2.Draw("th:same");
216
217
218	for (int i = 0; i < 5; i++) {
219	fring[i]->fV2.Print("s");
220	fring[i]->fTrue.Print("s");
221	}
222	fV2.Print("s");
223	gROOT->GetListOfSpecials()->Add(this);
224	return kTRUE;
225	}
226	void ToFile(const char* file) {
227	TFile* out = TFile::Open(file, "RECREATE");
228	for (int i = 0; i < 5; i++) fRing[i]->Write(fRing[i]->Name());
229	fV2.Write();
230	if (fPsi) fPsi->Write();
231	if (fDPsi) fDPsi->Write();
232	if (fDPhi) fDPhi->Write();
233	out->Close();
234	}
235
236
237
238	//__________________________________________________________________
239	/** Get min pseudo rapidity of ring idx */
240	Float_t EtaMin(Int_t idx) const {
241	switch (idx) {
242	case 0: return 3.61728; break; // FMD1I
243	case 1: return 2.18122; break; // FMD2I
244	case 2: return 1.79821; break; // FMD2O
245	case 3: return -3.39913; break; // FMD3I
246	case 4: return -2.28925; break; // FMD3O
247	}
248	return 0;
249	}
250	//__________________________________________________________________
251	/** Get max pseudo rapidity of ring idx */
252	Float_t EtaMax(Int_t idx) const {
253	switch (idx) {
254	case 0: return 5.02643; break; // FMD1I
255	case 1: return 3.57899; break; // FMD2I
256	case 2: return 2.39085; break; // FMD2O
257	case 3: return -2.00644; break; // FMD3I
258	case 4: return -1.70080; break; // FMD3O
259	}
260	return 0;
261	}
262	//__________________________________________________________________
263	/** Get cut 1 of ring idx */
264	Float_t Cut1(Int_t idx) const {
265	switch (idx) {
266	case 0: return 1.33857; break; // FMD1I
267	case 1: return 1.30464; break; // FMD2I
268	case 2: return 1.27263; break; // FMD2O
269	case 3: return 1.28052; break; // FMD3I
270	case 4: return 1.23041; break; // FMD3O
271	}
272	return 0;
273	}
274	//__________________________________________________________________
275	/** Get background correction of ring idx */
276	TH2* Bg(Bool_t use, Int_t idx) const {
277	if (!use) return 0;
278	switch (idx) {
279	case 0: return bgFMD1i(); break;
280	case 1: return bgFMD2i(); break;
281	case 2: return bgFMD2o(); break;
282	case 3: return bgFMD3i(); break;
283	case 4: return bgFMD3o(); break;
284	}
285	return 0;
286	}
287	/** Flow axis */
288	AliFMDFlowAxis fAxis;
289	/** Flow object */
290	AliFMDFlowBinned1D fV2;
291	/** Cache of phi */
292	TArrayD fPhis;
293	/** Cache of eta */
294	TArrayD fEtas;
295	/** Cache of weights */
296	TArrayD fWeights;
297	/** Number of observables */
298	Int_t fNObserv;
299	/** Histogram of Psi */
300	TH2D* fPsi; // Event planes
301	/** Histogram of dPsi */
302	TProfile* fDPsi; // Mean distance to Psi
303	/** Histogram of dPhi */
304	TProfile2D* fDPhi;
305	ClassDef(AliFMDAnaFlow,1)
306	};
307
308	//____________________________________________________________________
309	//
310	// EOF
311	//