Subversion Repositories f9daq

Rev

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