Subversion Repositories f9daq

Rev

Rev 84 | 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*100);
  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.           double 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.           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");
  111.           cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
  112.           TH2F* histoGlass = (detector->GetHGlass());
  113.           histoGlass->Draw("COLZ");
  114.  
  115.           cdata->cd(cpc++);
  116.           TH2F* histoActive = (detector->GetHActive());
  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);
  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");
  128.  
  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();
  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.   }
  144. }
  145. //-----------------------------------------------------------------------------
  146. void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
  147.     double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
  148. {
  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;
  169. }
  170. //-----------------------------------------------------------------------------
  171. void PrintGlassStat(CDetector *detector)
  172. {
  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);
  176. }
  177. //-----------------------------------------------------------------------------
  178. void PrintGuideHead()
  179. {
  180.   //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  181.   printf("Acceptance\n");
  182. }
  183. //-----------------------------------------------------------------------------
  184. void PrintGuideStat(double acceptance)
  185. {
  186.   /*int n_active = (detector->GetHActive())->GetEntries();
  187.         int n_enter = (detector->GetLG())->GetEnteranceHits();
  188.         double izkoristek = n_active/(double)n_enter;
  189.         int fates[7]; (detector->GetLG())->GetVFate(fates);
  190.  
  191. //      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  192.  
  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);       
  196.  
  197.         printf("%5.1lf\n", 100.0*izkoristek);*/
  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);
  203. }
  204. //-----------------------------------------------------------------------------
  205.  
  206. //-----------------------------------------------------------------------------
  207. // en zarek
  208. double Single(CDetector *detector, DetectorParameters& parameters,
  209.     TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
  210. {
  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));
  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");
  230.   ray0->setPolarization(polarization);
  231.   CRay ray1;
  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.  
  242.   TVector3 outPolarization = ray1.GetP();
  243.   drawPosition = ray1.GetR();
  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;
  250.   //delete ray1;
  251.  
  252.   return (detector->GetHActive())->GetEntries() / (double)(1);
  253. }
  254. //-----------------------------------------------------------------------------
  255. // zarki, razporejeni v mrezi
  256. double Grid(CDetector *detector, DetectorParameters& parameters,
  257.     int NN = 10, double theta = 0.0)
  258. {
  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;
  262.  
  263.   //detector->Init();
  264.   if(show_3d) detector->Draw(draw_width);
  265.  
  266.   const double b = parameters.getB();
  267.   double simulated = 0;
  268.  
  269.   double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  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.           }
  293.       }
  294.  
  295.   }
  296.   double acceptance = 0.0;
  297.   /*
  298.         if( !(parameters.getPlateOn()) )
  299.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
  300.   else
  301.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
  302.    */
  303.   acceptance = (detector->GetHActive())->GetEntries() / simulated;
  304.   return acceptance;
  305. }
  306. //-----------------------------------------------------------------------------
  307. // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
  308. // vsi pod kotom (theta, phi)
  309. double RandYZ(CDetector *detector, DetectorParameters& parameters,
  310.     int NN, double theta, double phi, int show_rays)
  311. {
  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;
  319.   TRandom2 rand;
  320.   rand.SetSeed(now.Get());
  321.  
  322.   //detector->Init();
  323.   if(show_3d) detector->Draw(draw_width);
  324.  
  325.   double b = parameters.getB();
  326.   double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  327.  
  328.   for(int i = 0; i < NN; i++) {
  329.  
  330.       double rx = CENTER.x() - offset;
  331.       double ry = rand.Uniform(-b/2.0, b/2.0);
  332.       double rz = rand.Uniform(-b/2.0, b/2.0);
  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);
  339.       TVector3 k(ray0->GetK());
  340.       TVector3 polarization(k.Orthogonal());
  341.       polarization.Unit();
  342.       polarization.Rotate(phi, k);
  343.       ray0->setPolarization(polarization);
  344.       CRay ray1;
  345.  
  346.       if(i < show_rays)
  347.         detector->Propagate(*ray0, ray1, show_3d);
  348.       else
  349.         detector->Propagate(*ray0, ray1, 0);
  350.       //delete ray1;
  351.       delete ray0;
  352.   }
  353.  
  354.   double acceptance = 0.0;
  355.   /*
  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));*/
  361.   acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
  362.   return acceptance;
  363. }
  364. //-----------------------------------------------------------------------------
  365. // zarki, izotropno porazdeljeni znotraj kota theta
  366. // = nakljucni vstopni polozaj in kot
  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]
  371. double RandIso(CDetector *detector, DetectorParameters& parameters,
  372.     int NN = 1e3, double thetaMin=0.0, double thetaMax = 30.0, int show_rays = 30, int show_rand = 0)
  373. {
  374.   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  375.   double pi = 3.14159265358979312;
  376.   thetaMax = thetaMax*pi/180.0;
  377.   thetaMin = thetaMin*pi/180.0;
  378.   if(thetaMin < MARGIN) thetaMin = MARGIN;
  379.   if(thetaMax < MARGIN) thetaMax = MARGIN;
  380.  
  381.   TDatime now;
  382.   TRandom2 rand;
  383.   rand.SetSeed(now.Get());
  384.   //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  385.   double theta_min_rad = TMath::Cos(thetaMax);
  386.   double theta_max_rad = TMath::Cos(thetaMin);
  387.  
  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);
  401.  
  402.  
  403.   //detector->Init();
  404.   if (show_3d) detector->Draw(draw_width);
  405.  
  406.   double b = parameters.getB();
  407.   double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  408.  
  409.   for(int i = 0; i < NN; i++) {
  410.  
  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);
  414.  
  415.       double phi = rand.Uniform(0.0, 2.0*pi);
  416.       //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
  417.       //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  418.       double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
  419.  
  420.       double rl = TMath::Cos(theta);
  421.       double rm = TMath::Sin(theta)*TMath::Sin(phi);
  422.       double rn = TMath::Sin(theta)*TMath::Cos(phi);
  423.  
  424.       if(show_rand) {
  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);
  431.       }
  432.  
  433.       CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  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;
  440.  
  441.       if(i < show_rays) {
  442.           detector->Propagate(*ray0, ray1, show_3d);
  443.       }
  444.       else {
  445.           detector->Propagate(*ray0, ray1, 0);
  446.       }
  447.  
  448.       //delete ray1;
  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.   /*
  468.         if( !(parameters.getPlateOn()) )
  469.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  470.   else
  471.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  472.    */
  473.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  474.   return acceptance;
  475. }
  476.  
  477. // Beamtest distribution
  478. // with fixed theta and phi in interval phiMin, phiMax
  479. double beamtest(CDetector *detector, DetectorParameters& parameters,
  480.     int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
  481. {
  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;
  487.  
  488.   TDatime now;
  489.   TRandom2 rand;
  490.   rand.SetSeed(now.Get());
  491.   double rx, ry, rz, rl, rm, rn;
  492.   double phi;
  493.  
  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.  
  520.       phi = rand.Uniform(phiMin, phiMax);
  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);
  526.       rm = TMath::Sin(theta)*TMath::Sin(phi);
  527.       rn = TMath::Sin(theta)*TMath::Cos(phi);
  528.  
  529.       if(show_rand) {
  530.           //htheta->Fill(rtheta);
  531.           //hcostheta->Fill( TMath::Cos(rtheta) );
  532.           hphi->Fill(phi);
  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
  540.       TVector3 k(ray0->GetK());
  541.       TVector3 polarization(k.Orthogonal());
  542.       polarization.Unit();
  543.       polarization.Rotate(phi, k);
  544.       ray0->setPolarization(polarization);
  545.       CRay ray1;
  546.  
  547.       if(i < show_rays) {
  548.           detector->Propagate(*ray0, ray1, show_3d);
  549.       }
  550.       else {
  551.           detector->Propagate(*ray0, ray1, 0);
  552.       }
  553.  
  554.       //delete ray1;
  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.   /*
  574.         if( !(parameters.getPlateOn()) )
  575.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  576.   else
  577.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  578.    */
  579.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  580.   return acceptance;
  581. }
  582.