Subversion Repositories f9daq

Rev

Rev 54 | Go to most recent revision | Details | 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
12
void Show3D(int b) {show_3d = b;}
13
void ShowData(int b) {show_data = b;}
14
 
15
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap)
16
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.52, TVector3(0.3, 0, 0));
17
 
18
//void SetLGType(int in = 1, int side = 1, int out = 0)
19
        //{detector->SetLGType(in, side, out);}
20
 
21
void SetLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n10 = 1.0, double n20 = 1.53, double n30 = 1.52)
22
        { parameters.setGuide(SiPM0, b0, d0, n10, n20, n30); }
23
 
24
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
25
  { parameters.setGap(x_gap0, y_gap0, z_gap0); }
26
 
27
void SetGlass(double glassOn, double glassD)
28
  { parameters.setGlass(glassOn, glassD); };
29
 
30
void SetPlate(int plateOn = 1, double plateWidth = 1)
31
  { parameters.setPlate(plateOn, plateWidth); }
32
 
33
void SetFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
34
 
35
int save_ary = 0;
36
//------------------------------------------------------------------------------------------
37
void Help()
38
{
39
        printf("void SetCenter(x, y, z)\n");
40
        printf("void SetLGType(in, side, out)\n");
41
        printf("SURF_DUMMY 0, ");
42
        printf("SURF_REFRA 1, ");
43
        printf("SURF_REFLE 2, ");
44
        printf("SURF_TOTAL 3, ");
45
        printf("SURF_IMPER 4\n");
46
        printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n");
47
        printf("void SetR(R0)\n");
48
        printf("void SetGlass(glass_on0, glass_d0)\n");
49
        printf("void SetGap(x_gap0, y_gap0, z_gap0)\n");
50
 
51
        printf("LGG(NN, theta)\n");
52
        printf("LGR(NN, theta, phi)\n");
53
        printf("LGI(NN, theta)\n");
54
 
55
        printf("LGI_gap(NN, min, max, steps, theta)\n");
56
        printf("LGR_gap(NN, min, max, steps, theta, phi)\n");
57
        printf("LGR_th(NN, min, max, steps, phi)\n");
58
        printf("LGG_th(NN, min, max, steps)\n");
59
 
60
        printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)");  
61
}
62
//------------------------------------------------------------------------------------------
63
int show_in_steps = 0;
64
void SetShowInSteps(int in = 1)
65
{                                                                                
66
        show_in_steps = in;
67
}
68
//------------------------------------------------------------------------------------------
69
int gs_set = 0;
70
void SetGS()
71
{                                                                                
72
        const UInt_t Number = 2;
73
        Double_t Red[Number]   = { 1.00, 0.00};
74
        Double_t Green[Number] = { 1.00, 0.00};
75
        Double_t Blue[Number]  = { 1.00, 0.00};
76
        Double_t Stops[Number] = { 0.00, 1.00};
77
 
78
        Int_t nb=50;
79
        if(!gs_set) {
80
                TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
81
                gs_set = 1;
82
        }
83
}
84
//------------------------------------------------------------------------------------------
85
void Draw1(CDetector *detector, int gs = 0, double range = 0.0)
86
{
87
        if(gs) SetGS();
88
        TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500);
89
        gStyle->SetOptTitle(0);
90
        gStyle->SetOptStat(0);
91
        cdata->cd(0);
92
                if(range > 0.1) {
93
                        ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range);
94
                        ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range);
95
                }
96
                (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]");
97
                (detector->GetHLaser())->Draw("COLZ");
98
        cdata->SaveAs("2dscan.gif");
99
        //gStyle->SetOptTitle(1);
100
        //gStyle->SetOptStat(1);
101
}
102
 
103
//------------------------------------------------------------------------------------------
104
 
105
 
106
TVector3 p_pol0(0.0, 0.0, 1.0);
107
void SetPol(double x, double y, double z)
108
{p_pol0.SetXYZ(x,y,z);}
109
 
110
void PolTest(double theta = 0.0)
111
{
112
  int p_type = 1;
113
  double p_n1 = 1.5;
114
  double p_n2 = 1.0;
115
        theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
116
        TVector3 p_pol;
117
 
118
        Init();
119
 
120
        double cx = 0.0;
121
        TVector3 vodnik_edge[4];
122
        double t = 3.0;
123
        vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t);
124
        vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t);
125
        /*
126
#define SURF_DUMMY 0
127
#define SURF_REFRA 1
128
#define SURF_REFLE 2
129
#define SURF_TOTAL 3
130
#define SURF_IMPER 4
131
        */
132
        CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
133
        surf->SetFresnel(0);
134
        surf->Draw();
135
 
136
        CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
137
        ray->SetColor(1);
138
        //p_pol = rotatey(p_pol0, -theta);
139
        p_pol = p_pol0; p_pol.RotateY(-theta);
140
        //printf("p_pol = "); printv(p_pol);
141
        ray->SetPolarization(p_pol);
142
 
143
        ray->DrawS(cx, -5.0);
144
 
145
        CRay *out = new CRay; out->SetColor(2);
146
        TVector3 *inters = new TVector3;
147
        int fate = surf->PropagateRay(*ray, out, inters);
148
        if(fate == 1) out->DrawS(cx, 5.0);
149
 
150
        CRay *pp = new CRay; pp->SetColor(3);
151
        pp->Set(ray->GetR(), p_pol);
152
        pp->DrawS(cx, 2.0);
153
 
154
        CRay *pn = new CRay; pn->SetColor(4);
155
        pn->Set(ray->GetR(), surf->GetN());
156
        pn->DrawS(cx, 3.0);
157
}
158
 
159
void ptt()
160
{
161
        for(double th = 0.0; th < 91.0; th += 5.0) {
162
                printf("%lf ", th);
163
                PolTest(th);
164
        }
165
}
166
//------------------------------------------------------------------------------------------
167
// Propagate single ray and show the statistics
168
void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
169
{
170
  CDetector *detector = new CDetector(CENTER, parameters);
171
        Init();
172
        double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
173
        PrintGuideHead();
174
        PrintGuideStat(detector, 1);
175
        DrawData(detector, parameters, theta, izkoristek);
176
}
177
//------------------------------------------------------------------------------------------
178
// Propagate NN rays generated as grid and show the statistics
179
void LGG(int NN = 10, double theta = 0.0)
180
{
181
  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
182
        Init();
183
        double izkoristek = Grid(detector, parameters, NN, theta);
184
        //printf("izkoristek = %.3lf\n", izkoristek);
185
        PrintGuideHead();
186
        PrintGuideStat(detector, (NN+1)*(NN+1));
187
        DrawData(detector, parameters, theta, izkoristek);
188
}
189
//------------------------------------------------------------------------------------------
190
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
191
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
192
{
193
  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
194
        Init();
195
        double izkoristek = RandYZ(detector, parameters, NN, theta, phi);
196
        //printf("izkoristek = %.3lf\n", izkoristek);
197
        PrintGuideHead(); PrintGuideStat(detector, NN);
198
        DrawData(detector, parameters, theta, izkoristek);
199
}
200
//------------------------------------------------------------------------------------------
201
// Propagate NN rays isotropically generated in solid angle theta and show the statistics
202
void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
203
{
204
  Init();
205
  CDetector *detector = new CDetector(CENTER, parameters);
206
  //CDetector detector = new CDetector();
207
        double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr);
208
        //printf("izkoristek = %.3lf\n", izkoristek);
209
        PrintGuideHead(); PrintGuideStat(detector, NN);
210
        DrawData(detector, parameters, theta, izkoristek);
211
        //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
212
        //canvasDetector->cd();
213
        //detector->Draw();
214
        //delete detector;
215
}
216
 
217
 
218
 
219
//------------------------------------------------------------------------------------------
220
void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
221
{
222
int tmp3d = show_3d, tmpdata = show_data;
223
double step[steps], acc[steps];
224
 
225
        show_3d = 0; show_data = 0;
226
 
227
        //printf(" x_gap | acceptance \n");
228
        PrintGuideHead();
229
        for(int i=0; i<=steps; i++) {
230
                step[i] = min + i*(max - min)/(double)steps;
231
                parameters.setGap(step[i], 0.0, 0.0);
232
                CDetector *detector = new CDetector(CENTER, parameters);
233
                Init();
234
                acc[i] = RandIso(detector, parameters, NN, theta);
235
                //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
236
                PrintGuideStat(detector, NN);
237
        }
238
 
239
        //char sbuff[256];
240
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
241
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(), 
242
                //      (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
243
 
244
        DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
245
 
246
        show_3d = tmp3d; show_data = tmpdata;
247
}
248
//------------------------------------------------------------------------------------------
249
 
250
//------------------------------------------------------------------------------------------
251
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 5.0)
252
{
253
int tmp3d = show_3d, tmpdata = show_data;
254
double step[steps], acc[steps];
255
 
256
        show_3d = 0; show_data = 0;
257
 
258
        //printf(" theta | acceptance \n");
259
        PrintGuideHead();
260
        for(int i=0; i<=steps; i++) {
261
                Init();
262
                step[i] = min + i*(max - min)/(double)steps;
263
                CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
264
                acc[i] = RandYZ(detector, parameters, NN, step[i], phi);
265
                //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
266
                PrintGuideStat(detector, NN);
267
                delete detector;
268
        }
269
 
270
        //char sbuff[256];
271
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
272
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(), 
273
                //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
274
 
275
        DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
276
 
277
        show_3d = tmp3d; show_data = tmpdata;
278
}
279
//------------------------------------------------------------------------------------------
280
 
281
// a vs. d
282
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)
283
{
284
  int tmp3d = show_3d;
285
        int tmpdata = show_data;
286
        //char sbuff[256];
287
 
288
        show_3d = 0; // don't show simulations 
289
        show_data = 1;
290
 
291
        //double d  = detector->GetD();
292
        const double b = parameters.getB(); // upper side of LG
293
        //const double SiPM = 3.0; // the length of the detector itself
294
        //double reflectivity = detector->guide->getR();
295
 
296
        TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
297
 
298
        // Use the Fresnel eq. instead of fixed reflectivity 96%
299
        //detector->SetFresnel(1);
300
 
301
        printf("   d |   a  | Acceptance\n");
302
        for(int i=0; i<steps; i++) {
303
                const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
304
                for(int j=0; j<steps; j++) {
305
                  const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
306
                  const double M = b/a;
307
      parameters.setGuide(a, M, d);
308
                  CDetector *detector = new CDetector(CENTER, parameters);
309
                  //detector->guide->setLG(y, M, x);
310
                  //Init(); exclude simulation
311
                  double acceptance = RandIso(detector, parameters, NN, theta, 0, 0);
312
                  //double acceptance = Grid(NN, theta);
313
                  //double acceptance = -1.0;
314
                  hAcceptance->Fill(d, a, acceptance);
315
                  //printf("%.2lf | %.2lf | ", x,y);
316
                  //PrintGuideStat(detector, NN);
317
                  //delete detector; //works fine, 50x50 grid takes ~4MB of RAM
318
                  }
319
        }
320
 
321
 
322
 
323
        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
324
        TVirtualPad *pacc = cp->cd(0);
325
 
326
        pacc->SetRightMargin(0.10);
327
        pacc->SetLeftMargin(0.10);
328
        pacc->SetTopMargin(0.10);
329
 
330
        TFile *file = new TFile("acceptance.root","RECREATE");
331
        hAcceptance->Write();
332
        file->Close();
333
        delete file;
334
 
335
        //hAcceptance->SetContour(100);
336
        //gStyle->SetPalette(1,0);
337
        hAcceptance->SetTitle(";d [mm];a [mm]");
338
        hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
339
        hAcceptance->Draw("COLZ");
340
 
341
        show_3d = tmp3d; show_data = tmpdata;
342
}
343
 
344
//------------------------------------------------------------------------------------------
345
// collection efficiency vs length
346
// a changes, b rests the same, d changes
347
// still can't use this function, it has hardcoded tan 10 deg
348
// !very stupid!
349
// new 3D graph NEEDED
350
// d vs a vs acceptance
351
void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
352
{
353
  TVector3 gapGuideSiPM(0.3, 0, 0);
354
  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
355
        int tmp3d = show_3d, tmpdata = show_data;
356
        double step[steps], acc[steps];
357
        double sipm_d, M0;
358
        char sbuff[256];
359
 
360
        show_3d = 1; show_data = 1;
361
 
362
        M0 = parameters.getM();
363
 
364
        PrintGuideHead();
365
        for(int i=0; i<=steps; i++) {
366
                step[i] = min + i*(max - min)/(double)steps;           
367
                sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
368
                parameters.setGuide(sipm_d, M0/sipm_d, step[i]);
369
                //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
370
                printf("%.1lf | %.2lf | ", step[i], sipm_d);
371
                Init();
372
                acc[i] = RandIso(detector, parameters, NN, theta);
373
                PrintGuideStat(detector, NN);
374
                sprintf(sbuff, "d%2d.gif", i);
375
                //c3dview->SaveAs(sbuff);
376
        }
377
 
378
 
379
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
380
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(), 
381
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
382
 
383
        //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
384
 
385
        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
386
        TVirtualPad *pacc = cp->cd(0);
387
 
388
        pacc->SetRightMargin(0.05);
389
        pacc->SetLeftMargin(0.13);
390
        pacc->SetTopMargin(0.08);
391
 
392
        TGraph *gacceptance = new TGraph(steps+1, step, acc);
393
        gacceptance->SetTitle("; d [mm];Collection efficiency");
394
        gacceptance->SetMarkerStyle(8);
395
        gacceptance->SetLineColor(12);
396
        gacceptance->SetLineWidth(1);
397
        (gacceptance->GetXaxis())->SetRangeUser(min, max);
398
        //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
399
 
400
        gacceptance->Draw("ACP");
401
 
402
 
403
        show_3d = tmp3d; show_data = tmpdata;
404
}
405
//------------------------------------------------------------------------------------------
406
// Acceptance vs. the length
407
// The magnification ratio M rests the same all the time
408
void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
409
{
410
        int tmp3d = show_3d, tmpdata = show_data;
411
        double step[steps], acc[steps];
412
        const double a = parameters.getA();
413
        const double magnif = parameters.getM();
414
 
415
        //show_3d = show_in_steps; show_data = 1;
416
 
417
        // Use Fresnel equations
418
        //detector->SetFresnel(1);
419
 
420
        // Set glass (n=1.5) at the exit
421
        //detector->SetGlass(1,0);
422
 
423
        PrintGuideHead();
424
        for(int i=0; i<=steps; i++) {
425
                step[i] = min + i*(max - min)/(double)steps;
426
                parameters.setGuide(a, magnif, step[i]);
427
                CDetector *detector = new CDetector(CENTER, parameters);
428
                //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
429
                printf("%.1lf | ", step[i]);
430
                Init();
431
                acc[i] = RandIso(detector, parameters, NN, theta);
432
                PrintGuideStat(detector, NN);
433
                delete detector;
434
        }
435
 
436
        //char sbuff[256];
437
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
438
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(), 
439
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
440
 
441
        DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
442
 
443
        show_3d = tmp3d; show_data = tmpdata;
444
}
445
//------------------------------------------------------------------------------------------
446
// Magnification optimization. a and d are fixed, M and b change in steps
447
void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
448
{
449
        int tmp3d = show_3d, tmpdata = show_data;
450
        double step[steps], acc[steps];
451
 
452
        const double a = parameters.getA();
453
        const double d = parameters.getD();
454
 
455
        show_3d = 0; show_data = 0;
456
 
457
        PrintGuideHead();
458
        for(int i=0; i<=steps; i++) {
459
                step[i] = min + i*(max - min)/(double)steps;
460
                parameters.setGuide(a, step[i], d);
461
                CDetector *detector = new CDetector(CENTER, parameters);
462
                Init();
463
                acc[i] = RandIso(detector, parameters, NN, theta);
464
                PrintGuideStat(detector, NN);
465
        }
466
 
467
        //char sbuff[256];
468
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
469
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(), 
470
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
471
 
472
        DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
473
 
474
        show_3d = tmp3d; show_data = tmpdata;
475
}
476
//------------------------------------------------------------------------------------------
477