Subversion Repositories f9daq

Rev

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

  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.  
  26. // calls default constructor CVodnik.cpp
  27.  
  28. //void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
  29. //{center.SetXYZ(x,y,z); detector = new CDetector(center);}
  30.  
  31.  
  32. //void SetLGType(int in = 1, int side = 1, int out = 0)
  33. //{detector->SetLGType(in, side, out);}
  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)
  36. //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
  37.  
  38. //void SetR(double R0) {detector->SetR(R0);}
  39. /*
  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);}
  45.  
  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.  
  51.  
  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);}
  61.  */
  62. //-----------------------------------------------------------------------------
  63. void Init(double rsys_scale = 10.0)
  64. {      
  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.   }
  85. }
  86. //-----------------------------------------------------------------------------
  87. void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
  88. {
  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.   }
  117. }
  118. //-----------------------------------------------------------------------------
  119. void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
  120.     double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
  121. {
  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;
  142. }
  143. //-----------------------------------------------------------------------------
  144. void PrintGlassStat(CDetector *detector)
  145. {
  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);
  149. }
  150. //-----------------------------------------------------------------------------
  151. void PrintGuideHead()
  152. {
  153.   //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  154.   printf("Acceptance\n");
  155. }
  156. //-----------------------------------------------------------------------------
  157. void PrintGuideStat(double acceptance)
  158. {
  159.   /*int n_active = (detector->GetHActive())->GetEntries();
  160.         int n_enter = (detector->GetLG())->GetEnteranceHits();
  161.         double izkoristek = n_active/(double)n_enter;
  162.         int fates[7]; (detector->GetLG())->GetVFate(fates);
  163.  
  164. //      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  165.  
  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);       
  169.  
  170.         printf("%5.1lf\n", 100.0*izkoristek);*/
  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);
  176. }
  177. //-----------------------------------------------------------------------------
  178.  
  179. //-----------------------------------------------------------------------------
  180. // en zarek
  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)
  183. {
  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);
  225. }
  226. //-----------------------------------------------------------------------------
  227. // zarki, razporejeni v mrezi
  228. double Grid(CDetector *detector, DetectorParameters& parameters,
  229. int NN = 10, double theta = 0.0)
  230. {
  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;
  234.  
  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.   /*
  265.         if( !(parameters.getPlateOn()) )
  266.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
  267.   else
  268.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
  269.    */
  270.   acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
  271.   return acceptance;
  272. }
  273. //-----------------------------------------------------------------------------
  274. // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
  275. // vsi pod kotom (theta, phi)
  276. double RandYZ(CDetector *detector, DetectorParameters& parameters,
  277. int NN, double theta, double phi, int show_rays)
  278. {
  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.   /*
  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));*/
  328.   acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
  329.   return acceptance;
  330. }
  331. //-----------------------------------------------------------------------------
  332. // zarki, izotropno porazdeljeni znotraj kota theta
  333. // = nakljucni vstopni polozaj in kot
  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]
  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)
  340. {
  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;
  345.  
  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);
  354.  
  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);
  368.  
  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.   /*
  429.         if( !(parameters.getPlateOn()) )
  430.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  431.   else
  432.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  433.    */
  434.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  435.   return acceptance;
  436. }
  437.  
  438. // Beamtest distribution
  439. // with fixed theta and phi in interval phiMin, phiMax
  440. double beamtest(CDetector *detector, DetectorParameters& parameters,
  441. int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
  442. {
  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;
  448.  
  449.   TDatime now;
  450.   TRandom rand;
  451.   rand.SetSeed(now.Get());
  452.   double rx, ry, rz, rl, rm, rn;
  453.   double rphi;
  454.  
  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.   /*
  534.         if( !(parameters.getPlateOn()) )
  535.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  536.   else
  537.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  538.    */
  539.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  540.   return acceptance;
  541. }
  542.