Rev 73 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
25 | f9daq | 1 | #include "include/guide.h" |
2 | #include "src/raySimulator.cpp" |
||
3 | |||
4 | #include "TFile.h" |
||
5 | #include "TPolyMarker3D.h" |
||
6 | #include "TCanvas.h" |
||
7 | |||
8 | //extern int show_3d; |
||
9 | //extern int show_data; |
||
10 | |||
11 | // Set the simulation parameters |
||
70 | f9daq | 12 | void showVisual(int b) {show_3d = b;} |
13 | void showData(int b) {show_data = b;} |
||
25 | f9daq | 14 | |
72 | f9daq | 15 | // Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap, badCoupling) |
16 | DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false); |
||
70 | f9daq | 17 | // Print the detector parameters |
18 | void getParameters(); |
||
25 | f9daq | 19 | |
20 | //void SetLGType(int in = 1, int side = 1, int out = 0) |
||
73 | f9daq | 21 | //{detector->SetLGType(in, side, out);} |
25 | f9daq | 22 | |
70 | f9daq | 23 | void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46) |
73 | f9daq | 24 | { parameters.setGuide(SiPM0, b0, d0); |
25 | parameters.setIndices(n1, n2, n3); } |
||
26 | |||
70 | f9daq | 27 | void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
73 | f9daq | 28 | { parameters.setGap(x_gap0, y_gap0, z_gap0); } |
29 | |||
70 | f9daq | 30 | void setGlass(double glassOn, double glassD) |
73 | f9daq | 31 | { parameters.setGlass(glassOn, glassD); }; |
25 | f9daq | 32 | |
70 | f9daq | 33 | void setPlate(int plateOn = 1, double plateWidth = 1) |
73 | f9daq | 34 | { parameters.setPlate(plateOn, plateWidth); } |
35 | |||
70 | f9daq | 36 | void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); } |
25 | f9daq | 37 | |
70 | f9daq | 38 | // Refractive Indices |
39 | // n1 - around light guide - air |
||
40 | // n2 - light guide (plate) material - k9 glass |
||
41 | // n3 - the material at the exit - optical grease, epoxy, air, etc. |
||
42 | void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);} |
||
43 | |||
25 | f9daq | 44 | int save_ary = 0; |
45 | //------------------------------------------------------------------------------------------ |
||
46 | void Help() |
||
47 | { |
||
73 | f9daq | 48 | printf("void SetCenter(x, y, z)\n"); |
49 | printf("void SetLGType(in, side, out)\n"); |
||
50 | printf("SURF_DUMMY 0, "); |
||
51 | printf("SURF_REFRA 1, "); |
||
52 | printf("SURF_REFLE 2, "); |
||
53 | printf("SURF_TOTAL 3, "); |
||
54 | printf("SURF_IMPER 4\n"); |
||
55 | printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n"); |
||
56 | printf("void SetR(R0)\n"); |
||
57 | printf("void SetGlass(glass_on0, glass_d0)\n"); |
||
58 | printf("void SetGap(x_gap0, y_gap0, z_gap0)\n"); |
||
59 | |||
60 | printf("LGG(NN, theta)\n"); |
||
61 | printf("LGR(NN, theta, phi)\n"); |
||
62 | printf("LGI(NN, theta)\n"); |
||
63 | |||
64 | printf("LGI_gap(NN, min, max, steps, theta)\n"); |
||
65 | printf("LGR_gap(NN, min, max, steps, theta, phi)\n"); |
||
66 | printf("LGR_th(NN, min, max, steps, phi)\n"); |
||
67 | printf("LGG_th(NN, min, max, steps)\n"); |
||
68 | |||
69 | printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)"); |
||
25 | f9daq | 70 | } |
71 | //------------------------------------------------------------------------------------------ |
||
72 | int show_in_steps = 0; |
||
73 | void SetShowInSteps(int in = 1) |
||
74 | { |
||
73 | f9daq | 75 | show_in_steps = in; |
25 | f9daq | 76 | } |
77 | //------------------------------------------------------------------------------------------ |
||
78 | int gs_set = 0; |
||
79 | void SetGS() |
||
80 | { |
||
73 | f9daq | 81 | const UInt_t Number = 2; |
82 | Double_t Red[Number] = { 1.00, 0.00}; |
||
83 | Double_t Green[Number] = { 1.00, 0.00}; |
||
84 | Double_t Blue[Number] = { 1.00, 0.00}; |
||
85 | Double_t Stops[Number] = { 0.00, 1.00}; |
||
86 | |||
87 | Int_t nb=50; |
||
88 | if(!gs_set) { |
||
89 | TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb); |
||
90 | gs_set = 1; |
||
91 | } |
||
25 | f9daq | 92 | } |
93 | //------------------------------------------------------------------------------------------ |
||
94 | void Draw1(CDetector *detector, int gs = 0, double range = 0.0) |
||
95 | { |
||
73 | f9daq | 96 | if(gs) SetGS(); |
97 | TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500); |
||
98 | gStyle->SetOptTitle(0); |
||
99 | gStyle->SetOptStat(0); |
||
100 | cdata->cd(0); |
||
101 | if(range > 0.1) { |
||
102 | ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range); |
||
103 | ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range); |
||
104 | } |
||
105 | (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]"); |
||
106 | (detector->GetHLaser())->Draw("COLZ"); |
||
107 | cdata->SaveAs("2dscan.gif"); |
||
108 | //gStyle->SetOptTitle(1); |
||
109 | //gStyle->SetOptStat(1); |
||
25 | f9daq | 110 | } |
111 | |||
112 | //------------------------------------------------------------------------------------------ |
||
113 | |||
70 | f9daq | 114 | // Test function |
115 | // creates an optical boundary surface |
||
116 | // shows the propagation of light ray |
||
117 | // and polarization state of the incident ray (green) |
||
118 | // and surface normal (blue) |
||
25 | f9daq | 119 | void PolTest(double theta = 0.0) |
120 | { |
||
73 | f9daq | 121 | int p_type = SURF_REFRA; |
70 | f9daq | 122 | double p_n1 = parameters.getN1(); |
123 | double p_n2 = parameters.getN2(); |
||
73 | f9daq | 124 | theta = 3.141593*theta/180.0; |
125 | if(theta < 1e-6) theta = 1e-6; |
||
126 | |||
127 | Init(); |
||
128 | |||
129 | double cx = 0.0; |
||
130 | TVector3 vodnik_edge[4]; |
||
84 | f9daq | 131 | double t = 4.0; |
73 | f9daq | 132 | vodnik_edge[0].SetXYZ(cx, t,-t); |
133 | vodnik_edge[1].SetXYZ(cx, t, t); |
||
134 | vodnik_edge[2].SetXYZ(cx,-t, t); |
||
135 | vodnik_edge[3].SetXYZ(cx,-t,-t); |
||
136 | |||
84 | f9daq | 137 | CSurface surf; |
138 | //surf = CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); |
||
139 | surf.Set(p_type, vodnik_edge, p_n1, p_n2, 0.96); |
||
140 | surf.FlipN(); |
||
141 | surf.SetFresnel(1); |
||
142 | surf.Draw(); |
||
73 | f9daq | 143 | |
144 | CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta)); |
||
145 | ray->SetColor(kRed); |
||
146 | //p_pol = rotatey(p_pol0, -theta); |
||
147 | TVector3 sE(0,1,0); |
||
148 | TVector3 pE(0,0,1); |
||
149 | TVector3 p_pol = pE; |
||
150 | p_pol.RotateY(-theta); |
||
151 | //printf("p_pol = "); printv(p_pol); |
||
152 | ray->setPolarization(p_pol); |
||
153 | |||
154 | ray->DrawS(cx, -5.0); |
||
155 | |||
84 | f9daq | 156 | CRay out; |
157 | out.SetColor(kBlack); |
||
73 | f9daq | 158 | TVector3 *inters = new TVector3; |
84 | f9daq | 159 | surf.PropagateRay(*ray, out, inters); |
73 | f9daq | 160 | printf(" n1 = %f, n2 = %f\n", p_n1, p_n2); |
161 | //if(fate == 1) out->DrawS(cx, 5.0); |
||
84 | f9daq | 162 | out.DrawS(cx, 5.0); |
73 | f9daq | 163 | |
164 | CRay *incidentPolarization = new CRay; |
||
165 | incidentPolarization->SetColor(kGreen); |
||
166 | incidentPolarization->Set(ray->GetR(), p_pol); |
||
167 | incidentPolarization->DrawS(cx, 1.0); |
||
168 | printf(" GREEN: polarization vector\n"); |
||
169 | |||
170 | CRay *surfaceNormal = new CRay; |
||
171 | surfaceNormal->SetColor(kBlue); |
||
84 | f9daq | 172 | surfaceNormal->Set(ray->GetR(), surf.GetN()); |
73 | f9daq | 173 | surfaceNormal->DrawS(cx, 1.0); |
174 | printf(" BLUE: surface normal vector\n"); |
||
25 | f9daq | 175 | } |
176 | |||
177 | void ptt() |
||
178 | { |
||
73 | f9daq | 179 | for(double th = 0.0; th < 91.0; th += 5.0) { |
180 | printf("%lf ", th); |
||
181 | PolTest(th); |
||
182 | } |
||
25 | f9daq | 183 | } |
184 | //------------------------------------------------------------------------------------------ |
||
185 | // Propagate single ray and show the statistics |
||
186 | void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0) |
||
187 | { |
||
188 | CDetector *detector = new CDetector(CENTER, parameters); |
||
73 | f9daq | 189 | Init(); |
190 | double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi); |
||
191 | PrintGuideHead(); |
||
192 | PrintGuideStat(izkoristek); |
||
193 | DrawData(detector, parameters, theta, izkoristek); |
||
25 | f9daq | 194 | } |
195 | //------------------------------------------------------------------------------------------ |
||
196 | // Propagate NN rays generated as grid and show the statistics |
||
72 | f9daq | 197 | void LGG(int NN=10, double theta=0.0, bool coupling=false) |
25 | f9daq | 198 | { |
72 | f9daq | 199 | parameters.setCoupling(coupling); |
73 | f9daq | 200 | CDetector *detector = new CDetector(CENTER, parameters); |
201 | Init(); |
||
202 | double izkoristek = Grid(detector, parameters, NN, theta); |
||
203 | //printf("izkoristek = %.3lf\n", izkoristek); |
||
204 | PrintGuideHead(); |
||
205 | PrintGuideStat(izkoristek); |
||
206 | DrawData(detector, parameters, theta, izkoristek); |
||
25 | f9daq | 207 | } |
208 | //------------------------------------------------------------------------------------------ |
||
209 | // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics |
||
210 | void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0) |
||
211 | { |
||
72 | f9daq | 212 | CDetector *detector = new CDetector(CENTER, parameters); |
73 | f9daq | 213 | Init(); |
214 | double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30); |
||
215 | //printf("izkoristek = %.3lf\n", izkoristek); |
||
216 | PrintGuideHead(); PrintGuideStat(izkoristek); |
||
217 | DrawData(detector, parameters, theta, izkoristek); |
||
25 | f9daq | 218 | } |
219 | //------------------------------------------------------------------------------------------ |
||
220 | // Propagate NN rays isotropically generated in solid angle theta and show the statistics |
||
221 | void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0) |
||
222 | { |
||
223 | Init(); |
||
224 | CDetector *detector = new CDetector(CENTER, parameters); |
||
225 | //CDetector detector = new CDetector(); |
||
73 | f9daq | 226 | double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr); |
227 | //printf("izkoristek = %.3lf\n", izkoristek); |
||
228 | PrintGuideHead(); |
||
229 | PrintGuideStat(izkoristek); |
||
230 | DrawData(detector, parameters, theta, izkoristek); |
||
231 | //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
||
232 | //canvasDetector->cd(); |
||
233 | //detector->Draw(); |
||
25 | f9daq | 234 | } |
235 | |||
72 | f9daq | 236 | void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0) |
70 | f9daq | 237 | { |
238 | Init(); |
||
72 | f9daq | 239 | parameters.setCoupling(coupling); |
70 | f9daq | 240 | CDetector *detector = new CDetector(CENTER, parameters); |
73 | f9daq | 241 | |
242 | const double theta = 18.5; |
||
243 | double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr); |
||
244 | |||
245 | PrintGuideHead(); |
||
246 | PrintGuideStat(izkoristek); |
||
247 | DrawData(detector, parameters, 18.5, izkoristek); |
||
248 | |||
249 | //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
||
250 | //canvasDetector->cd(); |
||
251 | //detector->Draw(); |
||
70 | f9daq | 252 | } |
73 | f9daq | 253 | |
254 | |||
25 | f9daq | 255 | //------------------------------------------------------------------------------------------ |
256 | void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0) |
||
257 | { |
||
73 | f9daq | 258 | int tmp3d = show_3d, tmpdata = show_data; |
259 | double step[steps], acc[steps]; |
||
25 | f9daq | 260 | |
73 | f9daq | 261 | show_3d = 0; show_data = 0; |
262 | |||
263 | //printf(" x_gap | acceptance \n"); |
||
264 | PrintGuideHead(); |
||
265 | for(int i=0; i<=steps; i++) { |
||
266 | step[i] = min + i*(max - min)/(double)steps; |
||
267 | parameters.setGap(step[i], 0.0, 0.0); |
||
268 | CDetector *detector = new CDetector(CENTER, parameters); |
||
269 | Init(); |
||
270 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
||
271 | //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0); |
||
272 | PrintGuideStat(acc[i]); |
||
273 | } |
||
274 | |||
275 | //char sbuff[256]; |
||
276 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
277 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
278 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
279 | |||
280 | DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max); |
||
281 | |||
282 | show_3d = tmp3d; show_data = tmpdata; |
||
25 | f9daq | 283 | } |
284 | //------------------------------------------------------------------------------------------ |
||
285 | |||
70 | f9daq | 286 | void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
287 | { |
||
288 | //int tmp3d = show_3d, tmpdata = show_data; |
||
289 | double show_rays = 10; |
||
290 | double step[steps], acc[steps]; |
||
291 | |||
73 | f9daq | 292 | //show_3d = 0; show_data = 0; |
293 | |||
294 | //printf(" theta | acceptance \n"); |
||
295 | PrintGuideHead(); |
||
296 | for(int i=0; i<=steps; i++) { |
||
297 | Init(); |
||
298 | step[i] = min + i*(max - min)/(double)steps; |
||
299 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
300 | acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays); |
||
301 | //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
||
302 | PrintGuideStat(acc[i]); |
||
303 | delete detector; |
||
304 | } |
||
305 | |||
306 | //char sbuff[256]; |
||
307 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)", |
||
308 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
309 | // (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z()); |
||
310 | |||
311 | DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
312 | |||
313 | //show_3d = tmp3d; show_data = tmpdata; |
||
70 | f9daq | 314 | } |
315 | |||
316 | |||
25 | f9daq | 317 | //------------------------------------------------------------------------------------------ |
70 | f9daq | 318 | //void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
319 | void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15) |
||
25 | f9daq | 320 | { |
70 | f9daq | 321 | //int tmp3d = show_3d, tmpdata = show_data; |
322 | double step[steps], acc[steps]; |
||
25 | f9daq | 323 | |
73 | f9daq | 324 | //show_3d = 0; show_data = 0; |
325 | |||
326 | //printf(" theta | acceptance \n"); |
||
327 | PrintGuideHead(); |
||
328 | for(int i=0; i<=steps; i++) { |
||
329 | Init(); |
||
330 | step[i] = min + i*(max - min)/(double)steps; |
||
331 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
332 | //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30); |
||
333 | acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d); |
||
334 | //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
||
335 | PrintGuideStat(acc[i]); |
||
336 | delete detector; |
||
337 | } |
||
338 | |||
339 | //char sbuff[256]; |
||
340 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)", |
||
341 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
342 | // (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z()); |
||
343 | |||
344 | DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
345 | |||
346 | //show_3d = tmp3d; show_data = tmpdata; |
||
25 | f9daq | 347 | } |
348 | //------------------------------------------------------------------------------------------ |
||
349 | |||
70 | f9daq | 350 | void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10) |
351 | { |
||
352 | int tmp3d = show_3d, tmpdata = show_data; |
||
353 | show_3d = 1; show_data = 0; |
||
354 | //PrintGuideHead(); |
||
355 | TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi); |
||
356 | double min = 0.0; |
||
357 | printf("\nWait, this takes a while "); |
||
73 | f9daq | 358 | for(int i=0; i<=steps; i++) { |
359 | double theta = min + i*(maxTheta - min) / (double)steps; |
||
360 | for (int j=0; j<=steps; j++) { |
||
361 | Init(); |
||
362 | double phi = min + j*(maxPhi - min) / (double)steps; |
||
363 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
364 | double acc = RandYZ(detector, parameters, NN, theta, phi, 30); |
||
365 | //PrintGuideStat(acc); |
||
366 | hAcceptance->Fill(i, j, acc); |
||
367 | //printf("Acc: %f ", acc); |
||
368 | delete detector; |
||
70 | f9daq | 369 | } |
73 | f9daq | 370 | printf("."); |
70 | f9daq | 371 | } |
372 | printf("\n"); |
||
373 | //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
374 | show_3d = tmp3d; show_data = tmpdata; |
||
375 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000); |
||
376 | cp->cd(); |
||
377 | hAcceptance->Draw("COLZ"); |
||
378 | } |
||
379 | |||
25 | f9daq | 380 | // a vs. d |
381 | void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0) |
||
382 | { |
||
73 | f9daq | 383 | //char sbuff[256]; |
25 | f9daq | 384 | |
73 | f9daq | 385 | //show_3d = 0; // don't show simulations |
386 | //show_data = 1; |
||
25 | f9daq | 387 | |
73 | f9daq | 388 | //double d = detector->GetD(); |
389 | const double b = parameters.getB(); // upper side of LG |
||
390 | //const double SiPM = 3.0; // the length of the detector itself |
||
391 | //double reflectivity = detector->guide->getR(); |
||
392 | |||
393 | TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max); |
||
394 | |||
395 | // Use the Fresnel eq. instead of fixed reflectivity 96% |
||
396 | //detector->SetFresnel(1); |
||
397 | |||
398 | //printf(" d | a | Acceptance\n"); |
||
399 | getParameters(); |
||
400 | printf("Wait, this takes a while "); |
||
401 | for(int i=0; i<steps; i++) { |
||
402 | const double d = hAcceptance->GetXaxis()->GetBinCenter(i); |
||
403 | for(int j=0; j<steps; j++) { |
||
404 | Init(); |
||
405 | const double a = hAcceptance->GetYaxis()->GetBinCenter(j); |
||
406 | //const double M = b/a; |
||
407 | parameters.setGuide(a, b, d); |
||
408 | CDetector *detector = new CDetector(CENTER, parameters); |
||
409 | //detector->guide->setLG(y, M, x); |
||
410 | //Init(); exclude simulation |
||
411 | double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d); |
||
412 | //double acceptance = Grid(NN, theta); |
||
413 | //double acceptance = -1.0; |
||
414 | hAcceptance->Fill(d, a, acceptance); |
||
415 | //printf("%.2lf | %.2lf | ", d, a); |
||
416 | //PrintGuideStat(acceptance); |
||
417 | delete detector; //works fine, 50x50 grid takes ~4MB of RAM |
||
418 | } |
||
419 | printf("."); |
||
420 | } |
||
421 | printf("\n"); |
||
422 | |||
423 | |||
424 | |||
425 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200); |
||
426 | cp->cd(); |
||
427 | //TVirtualPad *pacc = cp->cd(0); |
||
428 | |||
429 | /* |
||
430 | pacc->SetRightMargin(0.10); |
||
431 | pacc->SetLeftMargin(0.10); |
||
432 | pacc->SetTopMargin(0.10); |
||
433 | pacc->SetBottomMargin(0.10); |
||
434 | */ |
||
435 | |||
436 | TFile *file = new TFile("acceptance.root","RECREATE"); |
||
437 | hAcceptance->Write(); |
||
438 | file->Close(); |
||
439 | //delete file; |
||
440 | |||
441 | //hAcceptance->SetContour(100); |
||
442 | //gStyle->SetPalette(1,0); |
||
443 | hAcceptance->SetTitle(";d [mm];a [mm]"); |
||
444 | hAcceptance->GetXaxis()->SetRangeUser(minD,maxD); |
||
445 | hAcceptance->Draw("COLZ"); |
||
446 | char filename[128]; |
||
447 | sprintf(filename,"LGI_ad%d.C", steps); |
||
448 | cp->SaveAs(filename); |
||
449 | |||
25 | f9daq | 450 | } |
451 | |||
452 | //------------------------------------------------------------------------------------------ |
||
453 | // collection efficiency vs length |
||
454 | // a changes, b rests the same, d changes |
||
455 | // still can't use this function, it has hardcoded tan 10 deg |
||
456 | // !very stupid! |
||
457 | // new 3D graph NEEDED |
||
458 | // d vs a vs acceptance |
||
459 | void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0) |
||
460 | { |
||
461 | TVector3 gapGuideSiPM(0.3, 0, 0); |
||
462 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
73 | f9daq | 463 | int tmp3d = show_3d, tmpdata = show_data; |
464 | double step[steps], acc[steps]; |
||
465 | double sipm_d, M0; |
||
466 | char sbuff[256]; |
||
25 | f9daq | 467 | |
73 | f9daq | 468 | show_3d = 1; show_data = 1; |
469 | |||
470 | M0 = parameters.getM(); |
||
471 | |||
472 | PrintGuideHead(); |
||
473 | for(int i=0; i<=steps; i++) { |
||
474 | step[i] = min + i*(max - min)/(double)steps; |
||
475 | sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg |
||
476 | parameters.setGuide(sipm_d, M0*sipm_d, step[i]); |
||
477 | //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
||
478 | printf("%.1lf | %.2lf | ", step[i], sipm_d); |
||
479 | Init(); |
||
480 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
||
481 | PrintGuideStat(acc[i]); |
||
482 | sprintf(sbuff, "d%2d.gif", i); |
||
483 | //c3dview->SaveAs(sbuff); |
||
484 | } |
||
485 | |||
486 | |||
487 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
488 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
489 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
490 | |||
491 | //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max); |
||
492 | |||
493 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480); |
||
494 | TVirtualPad *pacc = cp->cd(0); |
||
495 | |||
496 | pacc->SetRightMargin(0.05); |
||
497 | pacc->SetLeftMargin(0.13); |
||
498 | pacc->SetTopMargin(0.08); |
||
499 | |||
500 | TGraph *gacceptance = new TGraph(steps+1, step, acc); |
||
501 | gacceptance->SetTitle("; d [mm];Collection efficiency"); |
||
502 | gacceptance->SetMarkerStyle(8); |
||
503 | gacceptance->SetLineColor(12); |
||
504 | gacceptance->SetLineWidth(1); |
||
505 | (gacceptance->GetXaxis())->SetRangeUser(min, max); |
||
506 | //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
||
507 | |||
508 | gacceptance->Draw("ACP"); |
||
509 | |||
510 | |||
511 | show_3d = tmp3d; show_data = tmpdata; |
||
25 | f9daq | 512 | } |
513 | //------------------------------------------------------------------------------------------ |
||
514 | // Acceptance vs. the length |
||
515 | // The magnification ratio M rests the same all the time |
||
516 | void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0) |
||
517 | { |
||
73 | f9daq | 518 | int tmp3d = show_3d, tmpdata = show_data; |
519 | double step[steps], acc[steps]; |
||
520 | const double a = parameters.getA(); |
||
521 | const double magnif = parameters.getM(); |
||
25 | f9daq | 522 | |
73 | f9daq | 523 | //show_3d = show_in_steps; show_data = 1; |
524 | |||
525 | // Use Fresnel equations |
||
526 | //detector->SetFresnel(1); |
||
527 | |||
528 | // Set glass (n=1.5) at the exit |
||
529 | //detector->SetGlass(1,0); |
||
530 | |||
531 | PrintGuideHead(); |
||
532 | for(int i=0; i<=steps; i++) { |
||
533 | step[i] = min + i*(max - min)/(double)steps; |
||
534 | //parameters.setGuide(a, magnif, step[i]); |
||
535 | parameters.setGuide(a, magnif, step[i]); |
||
536 | CDetector *detector = new CDetector(CENTER, parameters); |
||
537 | //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
||
538 | printf("%.1lf | ", step[i]); |
||
539 | Init(); |
||
540 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
||
541 | PrintGuideStat(acc[i]); |
||
542 | delete detector; |
||
543 | } |
||
544 | |||
545 | //char sbuff[256]; |
||
546 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
547 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
548 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
549 | |||
550 | DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max); |
||
551 | |||
552 | show_3d = tmp3d; show_data = tmpdata; |
||
25 | f9daq | 553 | } |
554 | //------------------------------------------------------------------------------------------ |
||
555 | // Magnification optimization. a and d are fixed, M and b change in steps |
||
556 | void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0) |
||
557 | { |
||
73 | f9daq | 558 | int tmp3d = show_3d, tmpdata = show_data; |
559 | double step[steps], acc[steps]; |
||
25 | f9daq | 560 | |
73 | f9daq | 561 | const double a = parameters.getA(); |
562 | const double d = parameters.getD(); |
||
563 | |||
564 | show_3d = 0; show_data = 0; |
||
565 | |||
566 | PrintGuideHead(); |
||
567 | for(int i=0; i<=steps; i++) { |
||
568 | step[i] = min + i*(max - min)/(double)steps; |
||
569 | //parameters.setGuide(a, step[i], d); |
||
570 | parameters.setGuide(a, a*step[i], d); |
||
571 | CDetector *detector = new CDetector(CENTER, parameters); |
||
572 | Init(); |
||
573 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
||
574 | PrintGuideStat(acc[i]); |
||
575 | } |
||
576 | |||
577 | //char sbuff[256]; |
||
578 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
579 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
580 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
581 | |||
582 | DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max); |
||
583 | |||
584 | show_3d = tmp3d; show_data = tmpdata; |
||
25 | f9daq | 585 | } |
70 | f9daq | 586 | |
25 | f9daq | 587 | //------------------------------------------------------------------------------------------ |
588 | |||
70 | f9daq | 589 | void getParameters() |
590 | { |
||
591 | printf("LIGHT GUIDE\n" |
||
73 | f9daq | 592 | " b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA()); |
70 | f9daq | 593 | printf("MATERIAL REFRACITVE INDICES\n" |
73 | f9daq | 594 | " n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3()); |
70 | f9daq | 595 | printf("PLATE\n" |
73 | f9daq | 596 | " ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel()); |
70 | f9daq | 597 | return; |
598 | } |
||
73 | f9daq | 599 |