Subversion Repositories f9daq

Rev

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