Some extra tests
[u/mrichter/AliRoot.git] / FMD / scripts / TestAcc.C
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(dx*dx+dy*dy);
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(dx*dx+dy*dy);
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 }