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

#include <TMath.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TMarker.h>
#include <TLine.h>
#include <TArc.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <THStack.h>

//_____________________________________________________________________
void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
{
  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  //AliFMDRing fmdring(ring);
  //fmdring.Init();
  // Float_t   rad       = pars->GetMaxR(r)-pars->GetMinR(r);
  // Float_t   slen      = pars->GetStripLength(r,t);
  // Float_t   sblen     = pars->GetBaseStripLength(r,t);
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };

  Double_t  minR      = (r == 'I' ?  4.5213 : 15.4);
  Double_t  maxR      = (r == 'I' ? 17.2    : 28.0);
  Double_t  rad       = maxR - minR;
  
  Int_t     nStrips   = (r == 'I' ? 512 : 256);
  Int_t     nSec      = (r == 'I' ?  20 :  40);
  Float_t   segment   = rad / nStrips;
  Float_t   basearc   = 2 * TMath::Pi() / (0.5 * nSec);
  Float_t   radius    = minR + t * segment;
  Float_t   baselen   = 0.5 * basearc * radius;

  const Double_t* c1  = (r == 'I' ? ic1 : oc1);
  const Double_t* c2  = (r == 'I' ? ic2 : oc2);

  // If the radius of the strip is smaller than the radius corresponding 
  // to the first corner we have a full strip length 
  Float_t   cr        = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
  if (radius <= cr) newm = 1;
  else {
    // Next, we should find the end-point of the strip - that is, 
    // the coordinates where the line from c1 to c2 intersects a circle 
    // with radius given by the strip. 
    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
    Float_t dx  = c2[0] - c1[0];
    Float_t dy  = c2[1] - c1[1];
    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
    Float_t det = radius * radius * dr * dr - D*D;
    
    if      (det <  0) newm = 1; // No intersection
    else if (det == 0) newm = 1; // Exactly tangent
    else {
      Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
      Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
      Float_t th  = TMath::ATan2(x, y);

      newm = th / (basearc/2);
    }
  }

  Float_t   slope     = (c1[1] - c2[1]) / (c1[0] - c2[0]);
  Float_t   constant  = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
  
  Float_t   d         = (TMath::Power(TMath::Abs(radius*slope),2) + 
			 TMath::Power(radius,2) - TMath::Power(constant,2));
  
  // If below corners return 1
  Float_t arclen = baselen;
  if (d > 0) {
    Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
			   (1+TMath::Power(slope, 2)));
    Float_t   y         = slope*x + constant;
    
    // If x is larger than corner x or y less than corner y, we have full
    // length strip
    if(x < c1[0] && y> c1[1]) {
      //One sector since theta is by definition half-hybrid
      Float_t   theta     = TMath::ATan2(x,y);
      arclen              = radius * theta;
    }
  }
  // Calculate the area of a strip with no cut
  Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
  Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
  Float_t   basearea  = basearea1 - basearea2;
  
  // Calculate the area of a strip with cut
  Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
  Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
  Float_t   area      = area1 - area2;
  
  // Acceptance is ratio 
  oldm = area/basearea;
}

void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
{
  TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ? 
			   "Inner" : "Outer", 800, 500);
  
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };

  const Double_t* c1   = (r == 'I' ? ic1 : oc1);
  const Double_t* c2   = (r == 'I' ? ic2 : oc2);
  Double_t  minR       = (r == 'I' ?  4.5213 : 15.4);
  Double_t  maxR       = (r == 'I' ? 17.2    : 28.0);
  Double_t  rad        = maxR - minR;
  Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
  Double_t  theta      = (r == 'I' ? 18 : 9);
  TH2*   frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
  frame->Draw();

  TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
  line->SetLineColor(kBlue+1);
  line->Draw();

  UShort_t nT = (r == 'I' ? 512 : 256);
  for (Int_t t = offT; t < nT; t += dt) {
    Double_t radius = minR + t * rad / nT;

    TArc*    circle = new TArc(0, 0, radius, 90-theta, 90+theta);
    circle->SetFillColor(0);
    circle->SetFillStyle(0);
    circle->Draw();

    // Now find the intersection 
    if (radius <= cr) continue;

    // Next, we should find the end-point of the strip - that is, 
    // the coordinates where the line from c1 to c2 intersects a circle 
    // with radius given by the strip. 
    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
    Float_t dx  = c2[0] - c1[0];
    Float_t dy  = c2[1] - c1[1];
    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
    Float_t det = radius * radius * dr * dr - D*D;
    
    if (det <  0) continue;
    if (det == 0) continue;


    Float_t x   = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
    Float_t y   = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
    
    TMarker* m = new TMarker(x, y, 20);
    m->SetMarkerColor(kRed+1);
    m->Draw();

    x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
    y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;

    // Float_t th = TMath::ATan2(x,y);
    // Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);

    m = new TMarker(x, y, 20);
    m->SetMarkerColor(kGreen+1);
    m->Draw();
  }
  c->cd();
}
    
  
void TestAcc()
{
  TCanvas* c =  new TCanvas("c", "C");
  TGraph*      innerOld = new TGraph(512);
  TGraph*      outerOld = new TGraph(256);
  TGraph*      innerNew = new TGraph(512);
  TGraph*      outerNew = new TGraph(256);
  innerOld->SetName("innerOld");
  innerOld->SetName("Inner type, old");
  innerOld->SetMarkerStyle(21);
  innerOld->SetMarkerColor(kRed+1);
  outerOld->SetName("outerOld");
  outerOld->SetName("Outer type, old");
  outerOld->SetMarkerStyle(21);
  outerOld->SetMarkerColor(kBlue+1);
  innerNew->SetName("innerNew");
  innerNew->SetName("Inner type, new");
  innerNew->SetMarkerStyle(20);
  innerNew->SetMarkerColor(kGreen+1);
  outerNew->SetName("outerNew");
  outerNew->SetName("Outer type, new");
  outerNew->SetMarkerStyle(20);
  outerNew->SetMarkerColor(kMagenta+1);

  TMultiGraph* all      = new TMultiGraph("all", "Acceptances");
  all->Add(innerOld);
  all->Add(outerOld);
  all->Add(innerNew);
  all->Add(outerNew);

  TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
  TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
  innerCor->SetMarkerStyle(20);
  outerCor->SetMarkerStyle(20);
  innerCor->SetMarkerColor(kRed+1);
  outerCor->SetMarkerColor(kBlue+1);
  // innerCor->SetMarkerSize(1.2);
  // outerCor->SetMarkerSize(1.2);
  THStack* stack = new THStack("corr", "Correlations");
  stack->Add(innerCor);
  stack->Add(outerCor);

  for (Int_t i = 0; i < 512; i++) { 
    Float_t oldAcc, newAcc;
    
    AcceptanceCorrection('I', i, oldAcc, newAcc);
    innerOld->SetPoint(i, i, oldAcc);
    innerNew->SetPoint(i, i, newAcc);
    // Printf("Inner strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);

    innerCor->Fill(oldAcc, newAcc);
    if (i >= 256) continue;
    
    AcceptanceCorrection('O', i, oldAcc, newAcc);
    outerOld->SetPoint(i, i, oldAcc);
    outerNew->SetPoint(i, i, newAcc);
    // Printf("Outer strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
    outerCor->Fill(oldAcc, newAcc);
  }

  all->Draw("ap");

  DrawSolution('I',4,256);
  DrawSolution('O',4,128);

  TCanvas* c2 = new TCanvas("cc", "cc");
  stack->Draw("nostack p");
  TLine* l = new TLine(0,0,1,1);
  l->SetLineStyle(2);
  l->Draw();

  c2->cd();
}
Commit	Line	Data
5cbf9b21	1	#include <TMath.h>
	2	#include <TGraph.h>
	3	#include <TMultiGraph.h>
	4	#include <TMarker.h>
	5	#include <TLine.h>
	6	#include <TArc.h>
	7	#include <TCanvas.h>
	8	#include <TH2F.h>
	9	#include <THStack.h>
	10
	11	//_____________________________________________________________________
	12	void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
	13	{
	14	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	15
	16	//AliFMDRing fmdring(ring);
	17	//fmdring.Init();
	18	// Float_t rad = pars->GetMaxR(r)-pars->GetMinR(r);
	19	// Float_t slen = pars->GetStripLength(r,t);
	20	// Float_t sblen = pars->GetBaseStripLength(r,t);
	21	const Double_t ic1[] = { 4.9895, 15.3560 };
	22	const Double_t ic2[] = { 1.8007, 17.2000 };
	23	const Double_t oc1[] = { 4.2231, 26.6638 };
	24	const Double_t oc2[] = { 1.8357, 27.9500 };
	25
	26	Double_t minR = (r == 'I' ? 4.5213 : 15.4);
	27	Double_t maxR = (r == 'I' ? 17.2 : 28.0);
	28	Double_t rad = maxR - minR;
	29
	30	Int_t nStrips = (r == 'I' ? 512 : 256);
	31	Int_t nSec = (r == 'I' ? 20 : 40);
	32	Float_t segment = rad / nStrips;
	33	Float_t basearc = 2 * TMath::Pi() / (0.5 * nSec);
	34	Float_t radius = minR + t * segment;
	35	Float_t baselen = 0.5 * basearc * radius;
	36
	37	const Double_t* c1 = (r == 'I' ? ic1 : oc1);
	38	const Double_t* c2 = (r == 'I' ? ic2 : oc2);
	39
	40	// If the radius of the strip is smaller than the radius corresponding
	41	// to the first corner we have a full strip length
	42	Float_t cr = TMath::Sqrt(c1[0]c1[0]+c1[1]c1[1]);
	43	if (radius <= cr) newm = 1;
	44	else {
	45	// Next, we should find the end-point of the strip - that is,
	46	// the coordinates where the line from c1 to c2 intersects a circle
	47	// with radius given by the strip.
	48	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
	49	Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
	50	Float_t dx = c2[0] - c1[0];
	51	Float_t dy = c2[1] - c1[1];
	52	Float_t dr = TMath::Sqrt(dxdx+dydy);
	53	Float_t det = radius * radius * dr * dr - D*D;
	54
	55	if (det < 0) newm = 1; // No intersection
	56	else if (det == 0) newm = 1; // Exactly tangent
	57	else {
	58	Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
	59	Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
	60	Float_t th = TMath::ATan2(x, y);
	61
	62	newm = th / (basearc/2);
	63	}
	64	}
65
66	Float_t slope = (c1[1] - c2[1]) / (c1[0] - c2[0]);
67	Float_t constant = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
68
69	Float_t d = (TMath::Power(TMath::Abs(radius*slope),2) +
70	TMath::Power(radius,2) - TMath::Power(constant,2));
71
72	// If below corners return 1
73	Float_t arclen = baselen;
74	if (d > 0) {
75	Float_t x = ((-TMath::Sqrt(d) - slope * constant) /
76	(1+TMath::Power(slope, 2)));
77	Float_t y = slope*x + constant;
78
79	// If x is larger than corner x or y less than corner y, we have full
80	// length strip
81	if(x < c1[0] && y> c1[1]) {
82	//One sector since theta is by definition half-hybrid
83	Float_t theta = TMath::ATan2(x,y);
84	arclen = radius * theta;
85	}
86	}
87	// Calculate the area of a strip with no cut
88	Float_t basearea1 = 0.5 * baselen * TMath::Power(radius,2);
89	Float_t basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
90	Float_t basearea = basearea1 - basearea2;
91
92	// Calculate the area of a strip with cut
93	Float_t area1 = 0.5 * arclen * TMath::Power(radius,2);
94	Float_t area2 = 0.5 * arclen * TMath::Power((radius-segment),2);
95	Float_t area = area1 - area2;
96
97	// Acceptance is ratio
98	oldm = area/basearea;
99	}
100
101	void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
102	{
103	TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ?
104	"Inner" : "Outer", 800, 500);
105
106	const Double_t ic1[] = { 4.9895, 15.3560 };
107	const Double_t ic2[] = { 1.8007, 17.2000 };
108	const Double_t oc1[] = { 4.2231, 26.6638 };
109	const Double_t oc2[] = { 1.8357, 27.9500 };
110
111	const Double_t* c1 = (r == 'I' ? ic1 : oc1);
112	const Double_t* c2 = (r == 'I' ? ic2 : oc2);
113	Double_t minR = (r == 'I' ? 4.5213 : 15.4);
114	Double_t maxR = (r == 'I' ? 17.2 : 28.0);
115	Double_t rad = maxR - minR;
116	Float_t cr = TMath::Sqrt(c1[0]c1[0]+c1[1]c1[1]);
117	Double_t theta = (r == 'I' ? 18 : 9);
118	TH2* frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
119	frame->Draw();
120
121	TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
122	line->SetLineColor(kBlue+1);
123	line->Draw();
124
125	UShort_t nT = (r == 'I' ? 512 : 256);
126	for (Int_t t = offT; t < nT; t += dt) {
127	Double_t radius = minR + t * rad / nT;
128
129	TArc* circle = new TArc(0, 0, radius, 90-theta, 90+theta);
130	circle->SetFillColor(0);
131	circle->SetFillStyle(0);
132	circle->Draw();
133
134	// Now find the intersection
135	if (radius <= cr) continue;
136
137	// Next, we should find the end-point of the strip - that is,
138	// the coordinates where the line from c1 to c2 intersects a circle
139	// with radius given by the strip.
140	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
141	Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
142	Float_t dx = c2[0] - c1[0];
143	Float_t dy = c2[1] - c1[1];
144	Float_t dr = TMath::Sqrt(dxdx+dydy);
145	Float_t det = radius * radius * dr * dr - D*D;
146
147	if (det < 0) continue;
148	if (det == 0) continue;
149
150
151	Float_t x = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
152	Float_t y = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
153
154	TMarker* m = new TMarker(x, y, 20);
155	m->SetMarkerColor(kRed+1);
156	m->Draw();
157
158	x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
159	y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
160
161	// Float_t th = TMath::ATan2(x,y);
162	// Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);
163
164	m = new TMarker(x, y, 20);
165	m->SetMarkerColor(kGreen+1);
166	m->Draw();
167	}
168	c->cd();
169	}
170
171
172	void TestAcc()
173	{
174	TCanvas* c = new TCanvas("c", "C");
175	TGraph* innerOld = new TGraph(512);
176	TGraph* outerOld = new TGraph(256);
177	TGraph* innerNew = new TGraph(512);
178	TGraph* outerNew = new TGraph(256);
179	innerOld->SetName("innerOld");
180	innerOld->SetName("Inner type, old");
181	innerOld->SetMarkerStyle(21);
182	innerOld->SetMarkerColor(kRed+1);
183	outerOld->SetName("outerOld");
184	outerOld->SetName("Outer type, old");
185	outerOld->SetMarkerStyle(21);
186	outerOld->SetMarkerColor(kBlue+1);
187	innerNew->SetName("innerNew");
188	innerNew->SetName("Inner type, new");
189	innerNew->SetMarkerStyle(20);
190	innerNew->SetMarkerColor(kGreen+1);
191	outerNew->SetName("outerNew");
192	outerNew->SetName("Outer type, new");
193	outerNew->SetMarkerStyle(20);
194	outerNew->SetMarkerColor(kMagenta+1);
195
196	TMultiGraph* all = new TMultiGraph("all", "Acceptances");
197	all->Add(innerOld);
198	all->Add(outerOld);
199	all->Add(innerNew);
200	all->Add(outerNew);
201
202	TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
203	TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
204	innerCor->SetMarkerStyle(20);
205	outerCor->SetMarkerStyle(20);
206	innerCor->SetMarkerColor(kRed+1);
207	outerCor->SetMarkerColor(kBlue+1);
208	// innerCor->SetMarkerSize(1.2);
209	// outerCor->SetMarkerSize(1.2);
210	THStack* stack = new THStack("corr", "Correlations");
211	stack->Add(innerCor);
212	stack->Add(outerCor);
213
214	for (Int_t i = 0; i < 512; i++) {
215	Float_t oldAcc, newAcc;
216
217	AcceptanceCorrection('I', i, oldAcc, newAcc);
218	innerOld->SetPoint(i, i, oldAcc);
219	innerNew->SetPoint(i, i, newAcc);
220	// Printf("Inner strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
221
222	innerCor->Fill(oldAcc, newAcc);
223	if (i >= 256) continue;
224
225	AcceptanceCorrection('O', i, oldAcc, newAcc);
226	outerOld->SetPoint(i, i, oldAcc);
227	outerNew->SetPoint(i, i, newAcc);
228	// Printf("Outer strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
229	outerCor->Fill(oldAcc, newAcc);
230	}
231
232	all->Draw("ap");
233
234	DrawSolution('I',4,256);
235	DrawSolution('O',4,128);
236
237	TCanvas* c2 = new TCanvas("cc", "cc");
238	stack->Draw("nostack p");
239	TLine* l = new TLine(0,0,1,1);
240	l->SetLineStyle(2);
241	l->Draw();
242
243	c2->cd();
244	}