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