Subversion Repositories f9daq

Rev

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