Subversion Repositories f9daq

Rev

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(),
85 f9daq 93
          parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc*100);
73 f9daq 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;
84 f9daq 99
          cdata->cd(cpc++);
100
          TH2F* generated = detector->GetGenerated();
85 f9daq 101
          double nGenerated = generated->GetEntries();
102
          //int minimum = generated->GetBinContent(generated->GetMinimumBin());
84 f9daq 103
          int maximum = generated->GetBinContent(generated->GetMaximumBin());
104
          generated->GetZaxis()->SetRangeUser(0, 1.05*maximum);
85 f9daq 105
          //double variation = (maximum-minimum)/(double)nGenerated;
84 f9daq 106
          generated->SetTitle("Generated");
107
          generated->Draw("colz");
108
          cdata->cd(cpc++); //((detector->GetLG())->GetHIn())->Draw("COLZ");
109
          TH2F* histoPlate = detector->GetHPlate();
110
          histoPlate->Draw("COLZ");
73 f9daq 111
          cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
84 f9daq 112
          TH2F* histoGlass = (detector->GetHGlass());
113
          histoGlass->Draw("COLZ");
73 f9daq 114
 
84 f9daq 115
          cdata->cd(cpc++);
116
          TH2F* histoActive = (detector->GetHActive());
85 f9daq 117
          //double nAccepted = histoActive->GetEntries();
118
          //double variation = 1 / sqrt(nAccepted);
119
          double variation = sqrt(acc*(1-acc)/nGenerated);
120
          printf("Statistical error = %f perc.\n", variation*100);
84 f9daq 121
          histoActive->Draw("COLZ");
122
          cdata->cd(cpc++);
123
          TH2F* histoLaser = (detector->GetHLaser());
124
          histoLaser->Draw("COLZ");
125
          cdata->cd(cpc++);
126
          TH2F* histoDetector = (detector->GetHDetector());
127
          histoDetector->Draw("COLZ");
73 f9daq 128
 
84 f9daq 129
          cdata->cd(cpc++);
130
          TH1F* fate=((detector->GetLG())->GetHFate());
131
          fate->Draw();
132
          cdata->cd(cpc++);
133
          TH1F* reflections = ((detector->GetLG())->GetHNOdb_all());
134
          reflections->Draw();
135
          cdata->cd(cpc++);
136
          TH1F *reflectedExit = ((detector->GetLG())->GetHNOdb_exit());
137
          reflectedExit->Draw();
73 f9daq 138
      }else{
139
          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
140
          cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
141
          cdata->SaveAs("cdata.gif");
142
      }
143
  }
25 f9daq 144
}
145
//-----------------------------------------------------------------------------
146
void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
73 f9daq 147
    double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
25 f9daq 148
{
73 f9daq 149
 
150
  cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
151
 
152
  (cacc->cd(0))->SetRightMargin(0.05);
153
  (cacc->cd(0))->SetLeftMargin(0.13);
154
  (cacc->cd(0))->SetTopMargin(0.08);
155
 
156
  TGraph *gacceptance = new TGraph(n, x, y);
157
  gacceptance->SetTitle(htitle);
158
  gacceptance->SetMarkerStyle(8);
159
  gacceptance->SetLineColor(12);
160
  gacceptance->SetLineWidth(1);
161
  if(xmin < xmax)
162
    (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
163
  (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
164
 
165
  gacceptance->Draw("ACP");
166
  //cacc->Print("acceptance.pdf");
167
  //delete gacceptance;
168
  //delete cacc;
25 f9daq 169
}
170
//-----------------------------------------------------------------------------
171
void PrintGlassStat(CDetector *detector)
172
{
73 f9daq 173
  int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
174
  int glass_out = (detector->GetHGlass())->GetEntries();
175
  printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
25 f9daq 176
}
177
//-----------------------------------------------------------------------------
178
void PrintGuideHead()
179
{
73 f9daq 180
  //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
181
  printf("Acceptance\n");
25 f9daq 182
}
183
//-----------------------------------------------------------------------------
54 f9daq 184
void PrintGuideStat(double acceptance)
25 f9daq 185
{
73 f9daq 186
  /*int n_active = (detector->GetHActive())->GetEntries();
25 f9daq 187
        int n_enter = (detector->GetLG())->GetEnteranceHits();
188
        double izkoristek = n_active/(double)n_enter;
189
        int fates[7]; (detector->GetLG())->GetVFate(fates);
73 f9daq 190
 
25 f9daq 191
//      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
73 f9daq 192
 
25 f9daq 193
        printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
194
                   (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
195
        for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       
73 f9daq 196
 
25 f9daq 197
        printf("%5.1lf\n", 100.0*izkoristek);*/
73 f9daq 198
 
199
  //double n_active = (detector->GetHActive())->GetEntries();
200
  //double r_acc = 100.0 * n_active / n_in;
201
  double r_acc = 100.0 * acceptance;
202
  printf("%7.3lf\n", r_acc);
25 f9daq 203
}
204
//-----------------------------------------------------------------------------
73 f9daq 205
 
25 f9daq 206
//-----------------------------------------------------------------------------
207
// en zarek
73 f9daq 208
double Single(CDetector *detector, DetectorParameters& parameters,
84 f9daq 209
    TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
25 f9daq 210
{
73 f9daq 211
  theta = Pi()*theta/180.0;
212
  if(theta < 1e-6) theta = 1e-6;
213
  phi = phi*Pi()/180.0;
214
  if(phi < 1e-6) phi = 1e-6;
215
 
216
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
217
 
218
  //detector->Init();
219
  if(show_3d) detector->Draw(draw_width);
220
 
221
  CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
222
      TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
223
      TMath::Sin(theta)*TMath::Cos(phi));
84 f9daq 224
  // Set y-polarization == vertical
225
  TVector3 k(ray0->GetK());
226
  TVector3 polarization(k.Orthogonal());
227
  polarization.Unit();
228
  polarization.Rotate(phi, k);
229
  if (polarization.Dot(ray0->GetK()) > 1e-2) printf("ERROR: pol not perep\n");
73 f9daq 230
  ray0->setPolarization(polarization);
84 f9daq 231
  CRay ray1;
73 f9daq 232
 
233
  detector->Propagate(*ray0, ray1, show_3d);
234
 
235
  CRay *incidentPolarization = new CRay;
236
  incidentPolarization->SetColor(kGreen);
237
  TVector3 drawPosition = ray0->GetR();
238
  drawPosition.SetX(drawPosition.X() - 5);
239
  incidentPolarization->Set(drawPosition, polarization);
240
  incidentPolarization->DrawS(drawPosition.X(), 1);
241
 
84 f9daq 242
  TVector3 outPolarization = ray1.GetP();
243
  drawPosition = ray1.GetR();
73 f9daq 244
  drawPosition.SetX(drawPosition.X() + 5);
245
  CRay* rayPol = new CRay(drawPosition, outPolarization);
246
  rayPol->SetColor(kBlack);
247
  rayPol->DrawS(drawPosition.X(), 1);
248
 
249
  delete ray0;
84 f9daq 250
  //delete ray1;
73 f9daq 251
 
252
  return (detector->GetHActive())->GetEntries() / (double)(1);
25 f9daq 253
}
254
//-----------------------------------------------------------------------------
255
// zarki, razporejeni v mrezi
73 f9daq 256
double Grid(CDetector *detector, DetectorParameters& parameters,
84 f9daq 257
    int NN = 10, double theta = 0.0)
25 f9daq 258
{
73 f9daq 259
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
260
  theta = Pi()*theta/180.0;
261
  if(theta < 1e-6) theta = 1e-6;
25 f9daq 262
 
73 f9daq 263
  //detector->Init();
264
  if(show_3d) detector->Draw(draw_width);
265
 
266
  const double b = parameters.getB();
84 f9daq 267
  double simulated = 0;
73 f9daq 268
 
269
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
84 f9daq 270
  const int GRID = 100;
271
  for(int i = 0; i < GRID; i++) {
272
      for(int j = 0; j < GRID; j++) {
273
          for(int jk = 0; jk<NN; jk++){
274
              CRay *ray0 = new CRay(CENTER.x() - offset,
275
                  (i*b/GRID + b/(2.0*GRID) - b/2.0),
276
                  (j*b/GRID + b/(2.0*GRID) - b/2.0),
277
                  TMath::Cos(theta),
278
                  0.0,
279
                  TMath::Sin(theta));
280
              TVector3 k(ray0->GetK());
281
              TVector3 polarization(k.Orthogonal());
282
              polarization.Unit();
283
              ray0->setPolarization(polarization);
284
              CRay ray1;
285
              if(i == (GRID/2) and jk == 0)
286
                detector->Propagate(*ray0, ray1, show_3d);
287
              else
288
                detector->Propagate(*ray0, ray1, 0);
289
              delete ray0;
290
              simulated++;
291
              //delete ray1;
292
          }
73 f9daq 293
      }
294
 
295
  }
296
  double acceptance = 0.0;
297
  /*
70 f9daq 298
        if( !(parameters.getPlateOn()) )
299
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
300
  else
301
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
73 f9daq 302
   */
84 f9daq 303
  acceptance = (detector->GetHActive())->GetEntries() / simulated;
73 f9daq 304
  return acceptance;
25 f9daq 305
}
306
//-----------------------------------------------------------------------------
307
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
308
// vsi pod kotom (theta, phi)
73 f9daq 309
double RandYZ(CDetector *detector, DetectorParameters& parameters,
84 f9daq 310
    int NN, double theta, double phi, int show_rays)
25 f9daq 311
{
73 f9daq 312
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
313
  theta = theta*3.14159265358979312/180.0;
314
  if(theta < MARGIN) theta = MARGIN;
315
  phi = phi*3.14159265358979312/180.0;
316
  if(phi < MARGIN) phi = MARGIN;
317
 
318
  TDatime now;
84 f9daq 319
  TRandom2 rand;
73 f9daq 320
  rand.SetSeed(now.Get());
321
 
322
  //detector->Init();
323
  if(show_3d) detector->Draw(draw_width);
324
 
84 f9daq 325
  double b = parameters.getB();
73 f9daq 326
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
327
 
328
  for(int i = 0; i < NN; i++) {
329
 
330
      double rx = CENTER.x() - offset;
84 f9daq 331
      double ry = rand.Uniform(-b/2.0, b/2.0);
332
      double rz = rand.Uniform(-b/2.0, b/2.0);
73 f9daq 333
 
334
      double rl = TMath::Cos(theta);
335
      double rm = TMath::Sin(theta)*TMath::Sin(phi);
336
      double rn = TMath::Sin(theta)*TMath::Cos(phi);
337
 
338
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
84 f9daq 339
      TVector3 k(ray0->GetK());
340
      TVector3 polarization(k.Orthogonal());
341
      polarization.Unit();
342
      polarization.Rotate(phi, k);
73 f9daq 343
      ray0->setPolarization(polarization);
84 f9daq 344
      CRay ray1;
73 f9daq 345
 
346
      if(i < show_rays)
347
        detector->Propagate(*ray0, ray1, show_3d);
348
      else
349
        detector->Propagate(*ray0, ray1, 0);
84 f9daq 350
      //delete ray1;
73 f9daq 351
      delete ray0;
352
  }
353
 
354
  double acceptance = 0.0;
355
  /*
70 f9daq 356
        if( !(parameters.getPlateOn()) )
357
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
358
  else
359
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
360
        //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/
73 f9daq 361
  acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
362
  return acceptance;
25 f9daq 363
}
364
//-----------------------------------------------------------------------------
365
// zarki, izotropno porazdeljeni znotraj kota theta
366
// = nakljucni vstopni polozaj in kot
70 f9daq 367
// NN - number of rays to be simulated with angles theta distributed uniformly
368
//      in cos theta and phi uniformly:
369
//      theta [0, theta]
370
//      phi [0,360]
73 f9daq 371
double RandIso(CDetector *detector, DetectorParameters& parameters,
84 f9daq 372
    int NN = 1e3, double thetaMin=0.0, double thetaMax = 30.0, int show_rays = 30, int show_rand = 0)
25 f9daq 373
{
73 f9daq 374
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
375
  double pi = 3.14159265358979312;
84 f9daq 376
  thetaMax = thetaMax*pi/180.0;
377
  thetaMin = thetaMin*pi/180.0;
378
  if(thetaMin < MARGIN) thetaMin = MARGIN;
379
  if(thetaMax < MARGIN) thetaMax = MARGIN;
25 f9daq 380
 
73 f9daq 381
  TDatime now;
84 f9daq 382
  TRandom2 rand;
73 f9daq 383
  rand.SetSeed(now.Get());
384
  //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
84 f9daq 385
  double theta_min_rad = TMath::Cos(thetaMax);
73 f9daq 386
  double theta_max_rad = TMath::Cos(thetaMin);
25 f9daq 387
 
73 f9daq 388
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
389
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
390
  hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
391
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
392
  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
393
  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
394
  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
395
  hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
396
  hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
397
  hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
398
  hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
399
  hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
400
  hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
25 f9daq 401
 
73 f9daq 402
 
403
  //detector->Init();
404
  if (show_3d) detector->Draw(draw_width);
405
 
84 f9daq 406
  double b = parameters.getB();
73 f9daq 407
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
408
 
409
  for(int i = 0; i < NN; i++) {
410
 
84 f9daq 411
      double rx = CENTER.x() - offset;
412
      double ry = rand.Uniform(-b/2.0, b/2.0);
413
      double rz = rand.Uniform(-b/2.0, b/2.0);
73 f9daq 414
 
84 f9daq 415
      double phi = rand.Uniform(0.0, 2.0*pi);
73 f9daq 416
      //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
417
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
84 f9daq 418
      double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
73 f9daq 419
 
84 f9daq 420
      double rl = TMath::Cos(theta);
421
      double rm = TMath::Sin(theta)*TMath::Sin(phi);
422
      double rn = TMath::Sin(theta)*TMath::Cos(phi);
73 f9daq 423
 
424
      if(show_rand) {
84 f9daq 425
          htheta->Fill(theta);
426
          hcostheta->Fill( TMath::Cos(theta) );
427
          hphi->Fill(phi);
428
          hl->Fill(rl);
429
          hm->Fill(rm);
430
          hn->Fill(rn);
73 f9daq 431
      }
432
 
433
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
84 f9daq 434
      TVector3 k(ray0->GetK());
435
      TVector3 polarization(k.Orthogonal());
436
      polarization.Unit();
437
      polarization.Rotate(phi, k);
438
      ray0->setPolarization(polarization);
439
      CRay ray1;
73 f9daq 440
 
441
      if(i < show_rays) {
442
          detector->Propagate(*ray0, ray1, show_3d);
443
      }
444
      else {
445
          detector->Propagate(*ray0, ray1, 0);
446
      }
447
 
84 f9daq 448
      //delete ray1;
73 f9daq 449
      delete ray0;
450
  }
451
 
452
  if(show_rand) {
453
      TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
454
      if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
455
      else c2rand->Clear();
456
      c2rand->Divide(2,3);
457
 
458
      c2rand->cd(1); hphi->Draw();
459
      c2rand->cd(3); htheta->Draw();
460
      c2rand->cd(5); hcostheta->Draw();
461
      c2rand->cd(2); hl->Draw();
462
      c2rand->cd(4); hm->Draw();
463
      c2rand->cd(6); hn->Draw();
464
  }
465
 
466
  double acceptance = 0.0;
467
  /*
70 f9daq 468
        if( !(parameters.getPlateOn()) )
469
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
470
  else
471
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
73 f9daq 472
   */
70 f9daq 473
  acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
73 f9daq 474
  return acceptance;
70 f9daq 475
}
476
 
477
// Beamtest distribution
478
// with fixed theta and phi in interval phiMin, phiMax
73 f9daq 479
double beamtest(CDetector *detector, DetectorParameters& parameters,
84 f9daq 480
    int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
70 f9daq 481
{
73 f9daq 482
  double pi = 3.14159265358979312;
483
  theta *= pi/180.0;
484
  if(theta < MARGIN) theta = MARGIN;
485
  phiMin *= pi/180.0;
486
  phiMax *= pi/180.0;
70 f9daq 487
 
73 f9daq 488
  TDatime now;
84 f9daq 489
  TRandom2 rand;
73 f9daq 490
  rand.SetSeed(now.Get());
491
  double rx, ry, rz, rl, rm, rn;
84 f9daq 492
  double phi;
70 f9daq 493
 
73 f9daq 494
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
495
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
496
  hphi = new TH1F("hphi", "hphi", 100, -pi, pi);
497
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
498
  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
499
  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
500
  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
501
  hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
502
  hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
503
  hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
504
  hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
505
  hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
506
  hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
507
 
508
  //detector->Init();
509
  if (show_3d) detector->Draw(draw_width);
510
 
511
  double b = parameters.getB();
512
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
513
 
514
  for(int i = 0; i < NN; i++) {
515
 
516
      rx = CENTER.x() - offset;
517
      ry = rand.Uniform(-b/2.0, b/2.0);
518
      rz = rand.Uniform(-b/2.0, b/2.0);
519
 
84 f9daq 520
      phi = rand.Uniform(phiMin, phiMax);
73 f9daq 521
      //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
522
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
523
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
524
 
525
      rl = TMath::Cos(theta);
84 f9daq 526
      rm = TMath::Sin(theta)*TMath::Sin(phi);
527
      rn = TMath::Sin(theta)*TMath::Cos(phi);
73 f9daq 528
 
529
      if(show_rand) {
530
          //htheta->Fill(rtheta);
531
          //hcostheta->Fill( TMath::Cos(rtheta) );
84 f9daq 532
          hphi->Fill(phi);
73 f9daq 533
          hl->Fill(rl);
534
          hm->Fill(rm);
535
          hn->Fill(rn);
536
      }
537
 
538
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
539
      // inital polarizaton
84 f9daq 540
      TVector3 k(ray0->GetK());
541
      TVector3 polarization(k.Orthogonal());
542
      polarization.Unit();
543
      polarization.Rotate(phi, k);
73 f9daq 544
      ray0->setPolarization(polarization);
84 f9daq 545
      CRay ray1;
73 f9daq 546
 
547
      if(i < show_rays) {
548
          detector->Propagate(*ray0, ray1, show_3d);
549
      }
550
      else {
551
          detector->Propagate(*ray0, ray1, 0);
552
      }
553
 
84 f9daq 554
      //delete ray1;
73 f9daq 555
      delete ray0;
556
  }
557
 
558
  if(show_rand) {
559
      TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
560
      if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
561
      else c2rand->Clear();
562
      c2rand->Divide(2,3);
563
 
564
      c2rand->cd(1); hphi->Draw();
565
      c2rand->cd(3); htheta->Draw();
566
      c2rand->cd(5); hcostheta->Draw();
567
      c2rand->cd(2); hl->Draw();
568
      c2rand->cd(4); hm->Draw();
569
      c2rand->cd(6); hn->Draw();
570
  }
571
 
572
  double acceptance = 0.0;
573
  /*
70 f9daq 574
        if( !(parameters.getPlateOn()) )
575
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
576
  else
577
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
73 f9daq 578
   */
70 f9daq 579
  acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
73 f9daq 580
  return acceptance;
25 f9daq 581
}