Subversion Repositories f9daq

Rev

Rev 71 | Rev 84 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
25 f9daq 1
// own include
2
//#include "include/raySimulator.h"
3
#include "include/RTUtil.h"
4
#include "include/guide.h"
5
 
6
// general include
7
#include "TROOT.h"
8
#include "TSystem.h"
9
#include "TStyle.h"
10
#include "TCanvas.h"
11
#include "TMath.h"
12
#include "TView3D.h"
13
#include "TH1F.h"
14
#include "TH2F.h"
15
#include "TBenchmark.h"
16
#include "TPolyMarker.h"
17
#include "TGraph.h"
18
#include "TF1.h"
19
 
20
//=================================================================================
21
TCanvas *c3dview;
22
TCanvas *cacc;
23
int show_3d = 0, show_data = 0;
24
int draw_width = 2;
25
 
73 f9daq 26
// calls default constructor CVodnik.cpp
25 f9daq 27
 
28
//void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
73 f9daq 29
//{center.SetXYZ(x,y,z); detector = new CDetector(center);}
25 f9daq 30
 
31
 
32
//void SetLGType(int in = 1, int side = 1, int out = 0)
73 f9daq 33
//{detector->SetLGType(in, side, out);}
25 f9daq 34
 
35
//void SetLG(double SiPM0 = 3.0, double M0 = 1.666, double d0 = 3.0, double n10 = 1.0, double n20 = 1.48, double n30 = 1.48, double R0 = 0.96)
73 f9daq 36
//{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
37
 
25 f9daq 38
//void SetR(double R0) {detector->SetR(R0);}
73 f9daq 39
/*
25 f9daq 40
void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)
41
        {detector->SetGlass(glass_on0, glass_d0);}
42
 
43
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
44
        {detector->SetGap(x_gap0 , y_gap0, z_gap0);}
73 f9daq 45
 
25 f9daq 46
void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)
47
        {detector->SetRCol(in, lg, out, gla);};
48
void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)
49
        {detector->SetDCol(LG0, glass0, active0);};
50
 
73 f9daq 51
 
25 f9daq 52
void SetDWidth(int w = 1) {draw_width = w;}
53
 
54
void SetFresnel(int b = 1) {detector->SetFresnel(b);}
55
void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}
56
 
57
void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}
58
 
59
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
60
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
73 f9daq 61
 */
25 f9daq 62
//-----------------------------------------------------------------------------
73 f9daq 63
void Init(double rsys_scale = 10.0)
25 f9daq 64
{      
73 f9daq 65
  RTSetStyle(gStyle);
66
  gStyle->SetOptStat(10);
67
 
68
  if(show_3d) {
69
      double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
70
      double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
71
 
72
      c3dview = (TCanvas*)gROOT->FindObject("c3dview");
73
      if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
74
      c3dview->SetFillColor(0);}
75
      else c3dview->Clear();
76
 
77
      TView3D *view = new TView3D(1, rsys0, rsys1);
78
      view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
79
      //view->SetPerspective();
80
      view->SetParallel();
81
      view->FrontView();
82
      view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
83
      //view->ToggleRulers();
84
  }
25 f9daq 85
}
86
//-----------------------------------------------------------------------------
87
void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
88
{
73 f9daq 89
  if(show_data) {
90
      char sbuff[256];
91
      sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
92
          parameters.getA(), parameters.getLightYield()*acc, parameters.getD(),
93
          parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
94
 
95
      if(!only2d) {
96
          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
97
          cdata->Divide(3,3);
98
          int cpc = 1;
99
          cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
100
          cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
101
          cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
102
          (detector->GetHGlass())->Draw("COLZ");
103
 
104
          cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
105
          cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
106
          cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
107
 
108
          cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
109
          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
110
          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
111
      }else{
112
          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
113
          cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
114
          cdata->SaveAs("cdata.gif");
115
      }
116
  }
25 f9daq 117
}
118
//-----------------------------------------------------------------------------
119
void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
73 f9daq 120
    double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
25 f9daq 121
{
73 f9daq 122
 
123
  cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
124
 
125
  (cacc->cd(0))->SetRightMargin(0.05);
126
  (cacc->cd(0))->SetLeftMargin(0.13);
127
  (cacc->cd(0))->SetTopMargin(0.08);
128
 
129
  TGraph *gacceptance = new TGraph(n, x, y);
130
  gacceptance->SetTitle(htitle);
131
  gacceptance->SetMarkerStyle(8);
132
  gacceptance->SetLineColor(12);
133
  gacceptance->SetLineWidth(1);
134
  if(xmin < xmax)
135
    (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
136
  (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
137
 
138
  gacceptance->Draw("ACP");
139
  //cacc->Print("acceptance.pdf");
140
  //delete gacceptance;
141
  //delete cacc;
25 f9daq 142
}
143
//-----------------------------------------------------------------------------
144
void PrintGlassStat(CDetector *detector)
145
{
73 f9daq 146
  int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
147
  int glass_out = (detector->GetHGlass())->GetEntries();
148
  printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
25 f9daq 149
}
150
//-----------------------------------------------------------------------------
151
void PrintGuideHead()
152
{
73 f9daq 153
  //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
154
  printf("Acceptance\n");
25 f9daq 155
}
156
//-----------------------------------------------------------------------------
54 f9daq 157
void PrintGuideStat(double acceptance)
25 f9daq 158
{
73 f9daq 159
  /*int n_active = (detector->GetHActive())->GetEntries();
25 f9daq 160
        int n_enter = (detector->GetLG())->GetEnteranceHits();
161
        double izkoristek = n_active/(double)n_enter;
162
        int fates[7]; (detector->GetLG())->GetVFate(fates);
73 f9daq 163
 
25 f9daq 164
//      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
73 f9daq 165
 
25 f9daq 166
        printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
167
                   (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
168
        for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       
73 f9daq 169
 
25 f9daq 170
        printf("%5.1lf\n", 100.0*izkoristek);*/
73 f9daq 171
 
172
  //double n_active = (detector->GetHActive())->GetEntries();
173
  //double r_acc = 100.0 * n_active / n_in;
174
  double r_acc = 100.0 * acceptance;
175
  printf("%7.3lf\n", r_acc);
25 f9daq 176
}
177
//-----------------------------------------------------------------------------
73 f9daq 178
 
25 f9daq 179
//-----------------------------------------------------------------------------
180
// en zarek
73 f9daq 181
double Single(CDetector *detector, DetectorParameters& parameters,
182
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
25 f9daq 183
{
73 f9daq 184
  theta = Pi()*theta/180.0;
185
  if(theta < 1e-6) theta = 1e-6;
186
  phi = phi*Pi()/180.0;
187
  if(phi < 1e-6) phi = 1e-6;
188
 
189
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
190
 
191
  //detector->Init();
192
  if(show_3d) detector->Draw(draw_width);
193
 
194
  CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
195
      TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
196
      TMath::Sin(theta)*TMath::Cos(phi));
197
  // Set z-polarization == vertical
198
  TVector3 polarization(0,0,1);
199
  polarization.RotateY(-theta);
200
  polarization.RotateZ(phi);
201
  if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n");
202
  ray0->setPolarization(polarization);
203
  CRay *ray1 = new CRay;
204
 
205
  detector->Propagate(*ray0, ray1, show_3d);
206
 
207
  CRay *incidentPolarization = new CRay;
208
  incidentPolarization->SetColor(kGreen);
209
  TVector3 drawPosition = ray0->GetR();
210
  drawPosition.SetX(drawPosition.X() - 5);
211
  incidentPolarization->Set(drawPosition, polarization);
212
  incidentPolarization->DrawS(drawPosition.X(), 1);
213
 
214
  TVector3 outPolarization = ray1->GetP();
215
  drawPosition = ray1->GetR();
216
  drawPosition.SetX(drawPosition.X() + 5);
217
  CRay* rayPol = new CRay(drawPosition, outPolarization);
218
  rayPol->SetColor(kBlack);
219
  rayPol->DrawS(drawPosition.X(), 1);
220
 
221
  delete ray0;
222
  delete ray1;
223
 
224
  return (detector->GetHActive())->GetEntries() / (double)(1);
25 f9daq 225
}
226
//-----------------------------------------------------------------------------
227
// zarki, razporejeni v mrezi
73 f9daq 228
double Grid(CDetector *detector, DetectorParameters& parameters,
229
int NN = 10, double theta = 0.0)
25 f9daq 230
{
73 f9daq 231
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
232
  theta = Pi()*theta/180.0;
233
  if(theta < 1e-6) theta = 1e-6;
25 f9daq 234
 
73 f9daq 235
  //detector->Init();
236
  if(show_3d) detector->Draw(draw_width);
237
 
238
  const double b = parameters.getB();
239
 
240
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
241
 
242
  for(int i = 0; i < NN; i++) {
243
      for(int j = 0; j < NN; j++) {
244
          CRay *ray0 = new CRay(CENTER.x() - offset,
245
              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
246
              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
247
              TMath::Cos(theta),
248
              0.0,
249
              TMath::Sin(theta));
250
          TVector3 polarization(0, 1, 0);
251
          polarization.RotateY(-theta);
252
          ray0->setPolarization(polarization);
253
          CRay *ray1 = new CRay;
254
          if(i == (NN/2))
255
            detector->Propagate(*ray0, ray1, show_3d);
256
          else
257
            detector->Propagate(*ray0, ray1, 0);
258
          delete ray0;
259
          delete ray1;
260
      }
261
 
262
  }
263
  double acceptance = 0.0;
264
  /*
70 f9daq 265
        if( !(parameters.getPlateOn()) )
266
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
267
  else
268
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
73 f9daq 269
   */
270
  acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
271
  return acceptance;
25 f9daq 272
}
273
//-----------------------------------------------------------------------------
274
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
275
// vsi pod kotom (theta, phi)
73 f9daq 276
double RandYZ(CDetector *detector, DetectorParameters& parameters,
277
int NN, double theta, double phi, int show_rays)
25 f9daq 278
{
73 f9daq 279
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
280
  theta = theta*3.14159265358979312/180.0;
281
  if(theta < MARGIN) theta = MARGIN;
282
  phi = phi*3.14159265358979312/180.0;
283
  if(phi < MARGIN) phi = MARGIN;
284
 
285
  TDatime now;
286
  TRandom rand;
287
  rand.SetSeed(now.Get());
288
 
289
  //detector->Init();
290
  if(show_3d) detector->Draw(draw_width);
291
 
292
  double SiPM = parameters.getA();
293
  double    M = parameters.getM();
294
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
295
 
296
  for(int i = 0; i < NN; i++) {
297
 
298
      double rx = CENTER.x() - offset;
299
      double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
300
      double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
301
 
302
      double rl = TMath::Cos(theta);
303
      double rm = TMath::Sin(theta)*TMath::Sin(phi);
304
      double rn = TMath::Sin(theta)*TMath::Cos(phi);
305
 
306
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
307
      TVector3 polarization(0, 0, 1);
308
      polarization.RotateY(-theta);
309
      polarization.RotateZ(phi);
310
      ray0->setPolarization(polarization);
311
      CRay *ray1 = new CRay;
312
 
313
      if(i < show_rays)
314
        detector->Propagate(*ray0, ray1, show_3d);
315
      else
316
        detector->Propagate(*ray0, ray1, 0);
317
      delete ray1;
318
      delete ray0;
319
  }
320
 
321
  double acceptance = 0.0;
322
  /*
70 f9daq 323
        if( !(parameters.getPlateOn()) )
324
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
325
  else
326
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
327
        //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/
73 f9daq 328
  acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
329
  return acceptance;
25 f9daq 330
}
331
//-----------------------------------------------------------------------------
332
// zarki, izotropno porazdeljeni znotraj kota theta
333
// = nakljucni vstopni polozaj in kot
70 f9daq 334
// NN - number of rays to be simulated with angles theta distributed uniformly
335
//      in cos theta and phi uniformly:
336
//      theta [0, theta]
337
//      phi [0,360]
73 f9daq 338
double RandIso(CDetector *detector, DetectorParameters& parameters,
339
int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
25 f9daq 340
{
73 f9daq 341
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
342
  double pi = 3.14159265358979312;
343
  theta = theta*3.14159265358979312/180.0;
344
  if(theta < MARGIN) theta = MARGIN;
25 f9daq 345
 
73 f9daq 346
  TDatime now;
347
  TRandom rand;
348
  rand.SetSeed(now.Get());
349
  double rx, ry, rz, rl, rm, rn;
350
  double rphi, rtheta;
351
  //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
352
  double theta_min_rad = TMath::Cos(theta);
353
  double theta_max_rad = TMath::Cos(thetaMin);
25 f9daq 354
 
73 f9daq 355
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
356
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
357
  hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
358
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
359
  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
360
  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
361
  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
362
  hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
363
  hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
364
  hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
365
  hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
366
  hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
367
  hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
25 f9daq 368
 
73 f9daq 369
 
370
  //detector->Init();
371
  if (show_3d) detector->Draw(draw_width);
372
 
373
  double SiPM = parameters.getA();
374
  double M = parameters.getM();
375
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
376
 
377
  for(int i = 0; i < NN; i++) {
378
 
379
      rx = CENTER.x() - offset;
380
      ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
381
      rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
382
 
383
      rphi = rand.Uniform(0.0, 2.0*pi);
384
      //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
385
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
386
      rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
387
 
388
      rl = TMath::Cos(rtheta);
389
      rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
390
      rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
391
 
392
      if(show_rand) {
393
          htheta->Fill(rtheta);
394
          hcostheta->Fill( TMath::Cos(rtheta) );
395
          hphi->Fill(rphi);
396
          hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
397
      }
398
 
399
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
400
      CRay *ray1 = new CRay;
401
 
402
      if(i < show_rays) {
403
          detector->Propagate(*ray0, ray1, show_3d);
404
      }
405
      else {
406
          detector->Propagate(*ray0, ray1, 0);
407
      }
408
 
409
      delete ray0;
410
      delete ray1;
411
  }
412
 
413
  if(show_rand) {
414
      TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
415
      if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
416
      else c2rand->Clear();
417
      c2rand->Divide(2,3);
418
 
419
      c2rand->cd(1); hphi->Draw();
420
      c2rand->cd(3); htheta->Draw();
421
      c2rand->cd(5); hcostheta->Draw();
422
      c2rand->cd(2); hl->Draw();
423
      c2rand->cd(4); hm->Draw();
424
      c2rand->cd(6); hn->Draw();
425
  }
426
 
427
  double acceptance = 0.0;
428
  /*
70 f9daq 429
        if( !(parameters.getPlateOn()) )
430
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
431
  else
432
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
73 f9daq 433
   */
70 f9daq 434
  acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
73 f9daq 435
  return acceptance;
70 f9daq 436
}
437
 
438
// Beamtest distribution
439
// with fixed theta and phi in interval phiMin, phiMax
73 f9daq 440
double beamtest(CDetector *detector, DetectorParameters& parameters,
441
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
70 f9daq 442
{
73 f9daq 443
  double pi = 3.14159265358979312;
444
  theta *= pi/180.0;
445
  if(theta < MARGIN) theta = MARGIN;
446
  phiMin *= pi/180.0;
447
  phiMax *= pi/180.0;
70 f9daq 448
 
73 f9daq 449
  TDatime now;
450
  TRandom rand;
451
  rand.SetSeed(now.Get());
452
  double rx, ry, rz, rl, rm, rn;
453
  double rphi;
70 f9daq 454
 
73 f9daq 455
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
456
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
457
  hphi = new TH1F("hphi", "hphi", 100, -pi, pi);
458
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
459
  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
460
  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
461
  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
462
  hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
463
  hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
464
  hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
465
  hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
466
  hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
467
  hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
468
 
469
  //detector->Init();
470
  if (show_3d) detector->Draw(draw_width);
471
 
472
  double b = parameters.getB();
473
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
474
 
475
  for(int i = 0; i < NN; i++) {
476
 
477
      rx = CENTER.x() - offset;
478
      ry = rand.Uniform(-b/2.0, b/2.0);
479
      rz = rand.Uniform(-b/2.0, b/2.0);
480
 
481
      rphi = rand.Uniform(phiMin, phiMax);
482
      //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
483
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
484
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
485
 
486
      rl = TMath::Cos(theta);
487
      rm = TMath::Sin(theta)*TMath::Sin(rphi);
488
      rn = TMath::Sin(theta)*TMath::Cos(rphi);
489
 
490
      if(show_rand) {
491
          //htheta->Fill(rtheta);
492
          //hcostheta->Fill( TMath::Cos(rtheta) );
493
          hphi->Fill(rphi);
494
          hl->Fill(rl);
495
          hm->Fill(rm);
496
          hn->Fill(rn);
497
      }
498
 
499
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
500
      // inital polarizaton
501
      TVector3 polarization(0, 0, 1);
502
      polarization.RotateY(-theta);
503
      polarization.RotateX(rphi);
504
      ray0->setPolarization(polarization);
505
      CRay *ray1 = new CRay;
506
 
507
      if(i < show_rays) {
508
          detector->Propagate(*ray0, ray1, show_3d);
509
      }
510
      else {
511
          detector->Propagate(*ray0, ray1, 0);
512
      }
513
 
514
      delete ray0;
515
      delete ray1;
516
  }
517
 
518
  if(show_rand) {
519
      TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
520
      if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
521
      else c2rand->Clear();
522
      c2rand->Divide(2,3);
523
 
524
      c2rand->cd(1); hphi->Draw();
525
      c2rand->cd(3); htheta->Draw();
526
      c2rand->cd(5); hcostheta->Draw();
527
      c2rand->cd(2); hl->Draw();
528
      c2rand->cd(4); hm->Draw();
529
      c2rand->cd(6); hn->Draw();
530
  }
531
 
532
  double acceptance = 0.0;
533
  /*
70 f9daq 534
        if( !(parameters.getPlateOn()) )
535
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
536
  else
537
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
73 f9daq 538
   */
70 f9daq 539
  acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
73 f9daq 540
  return acceptance;
25 f9daq 541
}