Subversion Repositories f9daq

Rev

Rev 70 | 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 = 7.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);      c3dview->SetFillColor(0);}
  74.                 else c3dview->Clear();
  75.                        
  76.                 TView3D *view = new TView3D(1, rsys0, rsys1);
  77.                 view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
  78.                 //view->SetPerspective();
  79.                 view->SetParallel();
  80.                 view->FrontView();
  81.                 view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
  82.                 //view->ToggleRulers();
  83.         }
  84. }
  85. //-----------------------------------------------------------------------------
  86. void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
  87. {
  88.         if(show_data) {
  89.                 char sbuff[256];
  90.                 sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
  91.                                 parameters.getA(), parameters.getLightYield()*acc, parameters.getD(),
  92.                                 parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
  93.                
  94.                 if(!only2d) {
  95.                         RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
  96.                         cdata->Divide(3,3);
  97.                         int cpc = 1;
  98.                         cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
  99.                         cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
  100.                         cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
  101.                         (detector->GetHGlass())->Draw("COLZ");
  102.                        
  103.                         cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
  104.                         cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
  105.                         cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
  106.                        
  107.                         cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
  108.                         cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
  109.                         cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
  110.                 }else{         
  111.                         RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
  112.                         cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
  113.                         cdata->SaveAs("cdata.gif");
  114.                 }
  115.         }
  116. }
  117. //-----------------------------------------------------------------------------
  118. void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
  119.                          double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
  120. {
  121.                
  122.         cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
  123.        
  124.         (cacc->cd(0))->SetRightMargin(0.05);
  125.         (cacc->cd(0))->SetLeftMargin(0.13);
  126.         (cacc->cd(0))->SetTopMargin(0.08);
  127.        
  128.         TGraph *gacceptance = new TGraph(n, x, y);
  129.         gacceptance->SetTitle(htitle);
  130.         gacceptance->SetMarkerStyle(8);
  131.         gacceptance->SetLineColor(12);
  132.         gacceptance->SetLineWidth(1);
  133.         if(xmin < xmax)
  134.                 (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
  135.         (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
  136.        
  137.         gacceptance->Draw("ACP");
  138.         //cacc->Print("acceptance.pdf");
  139.         //delete gacceptance;
  140.         //delete cacc;
  141. }
  142. //-----------------------------------------------------------------------------
  143. void PrintGlassStat(CDetector *detector)
  144. {
  145.         int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
  146.         int glass_out = (detector->GetHGlass())->GetEntries();
  147.         printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
  148. }
  149. //-----------------------------------------------------------------------------
  150. void PrintGuideHead()
  151. {
  152.         //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  153.         printf("Acceptance\n");
  154. }
  155. //-----------------------------------------------------------------------------
  156. void PrintGuideStat(double acceptance)
  157. {
  158.         /*int n_active = (detector->GetHActive())->GetEntries();
  159.         int n_enter = (detector->GetLG())->GetEnteranceHits();
  160.         double izkoristek = n_active/(double)n_enter;
  161.         int fates[7]; (detector->GetLG())->GetVFate(fates);
  162.        
  163. //      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  164.        
  165.         printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
  166.                    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
  167.         for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       
  168.        
  169.         printf("%5.1lf\n", 100.0*izkoristek);*/
  170.        
  171.         //double n_active = (detector->GetHActive())->GetEntries();
  172.         //double r_acc = 100.0 * n_active / n_in;
  173.         double r_acc = 100.0 * acceptance;
  174.         printf("%7.3lf\n", r_acc);
  175. }
  176. //-----------------------------------------------------------------------------
  177.        
  178. //-----------------------------------------------------------------------------
  179. // en zarek
  180. double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
  181. {
  182.   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  183.         theta = Pi()*theta/180.0;
  184.         if(theta < 1e-6) theta = 1e-6;
  185.         phi = phi*Pi()/180.0;
  186.         if(phi < 1e-6) phi = 1e-6;
  187.        
  188.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  189.        
  190.         //detector->Init();
  191.         if(show_3d) detector->Draw(draw_width);
  192.        
  193.         CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
  194.                           TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
  195.         CRay *ray1 = new CRay;
  196.        
  197.         detector->Propagate(*ray0, ray1, show_3d);
  198.        
  199.         delete ray0;
  200.         delete ray1;
  201.        
  202.         return (detector->GetHActive())->GetEntries() / (double)(1);
  203. }
  204. //-----------------------------------------------------------------------------
  205. // zarki, razporejeni v mrezi
  206. double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
  207. {
  208. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  209.         theta = Pi()*theta/180.0;
  210.         if(theta < 1e-6) theta = 1e-6;
  211.        
  212.         //detector->Init();
  213.         if(show_3d) detector->Draw(draw_width);
  214.  
  215.         const double b = parameters.getB();
  216.        
  217.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  218.        
  219.         for(int i = 0; i < NN; i++) {
  220.                 for(int j = 0; j < NN; j++) {
  221.                 CRay *ray0 = new CRay(CENTER.x() - offset,
  222.                               0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
  223.                               0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
  224.                                           TMath::Cos(theta),
  225.                                           0.0,
  226.                                           TMath::Sin(theta));
  227.                 CRay *ray1 = new CRay;
  228.                 if(i == (NN/2))
  229.                         detector->Propagate(*ray0, ray1, show_3d);
  230.                 else
  231.                         detector->Propagate(*ray0, ray1, 0);
  232.                 delete ray0;
  233.                 delete ray1;   
  234.                 }
  235.                
  236.         }
  237.         double acceptance = 0.0;
  238.         /*
  239.         if( !(parameters.getPlateOn()) )
  240.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
  241.   else
  242.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
  243.     */
  244.         acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
  245.         return acceptance;
  246. }
  247. //-----------------------------------------------------------------------------
  248. // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
  249. // vsi pod kotom (theta, phi)
  250. double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays)
  251. {
  252. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  253.         theta = theta*3.14159265358979312/180.0;
  254.         if(theta < MARGIN) theta = MARGIN;
  255.         phi = phi*3.14159265358979312/180.0;
  256.         if(phi < MARGIN) phi = MARGIN;
  257.        
  258.         TDatime now;
  259.         TRandom rand;
  260.         rand.SetSeed(now.Get());
  261.        
  262.         //detector->Init();
  263.         if(show_3d) detector->Draw(draw_width);
  264.        
  265.         double SiPM = parameters.getA();
  266.         double    M = parameters.getM();
  267.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  268.        
  269.         for(int i = 0; i < NN; i++) {
  270.                
  271.                 double rx = CENTER.x() - offset;
  272.                 double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  273.                 double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  274.                                
  275.                 double rl = TMath::Cos(theta);
  276.                 double rm = TMath::Sin(theta)*TMath::Sin(phi);
  277.                 double rn = TMath::Sin(theta)*TMath::Cos(phi); 
  278.                
  279.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  280.                 CRay *ray1 = new CRay;
  281.                
  282.                 if(i < show_rays)
  283.                         detector->Propagate(*ray0, ray1, show_3d);
  284.                 else
  285.                         detector->Propagate(*ray0, ray1, 0);
  286.                 //delete ray0;
  287.                 //delete ray1; 
  288.         }
  289.        
  290.         double acceptance = 0.0;
  291.         /*
  292.         if( !(parameters.getPlateOn()) )
  293.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  294.   else
  295.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
  296.         //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/
  297.         acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
  298.         return acceptance;     
  299. }
  300. //-----------------------------------------------------------------------------
  301. // zarki, izotropno porazdeljeni znotraj kota theta
  302. // = nakljucni vstopni polozaj in kot
  303. // NN - number of rays to be simulated with angles theta distributed uniformly
  304. //      in cos theta and phi uniformly:
  305. //      theta [0, theta]
  306. //      phi [0,360]
  307. double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
  308. {
  309. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  310.         double pi = 3.14159265358979312;
  311.         theta = theta*3.14159265358979312/180.0;
  312.         if(theta < MARGIN) theta = MARGIN;
  313.        
  314.         TDatime now;
  315.         TRandom rand;
  316.         rand.SetSeed(now.Get());
  317.         double rx, ry, rz, rl, rm, rn;
  318.         double rphi, rtheta;
  319.         //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  320.         double theta_min_rad = TMath::Cos(theta);
  321.         double theta_max_rad = TMath::Cos(thetaMin);
  322.  
  323.           TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
  324.           hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
  325.           hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
  326.           htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
  327.           htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
  328.           hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
  329.           hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
  330.           hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
  331.           hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
  332.           hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
  333.           hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
  334.           hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
  335.           hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
  336.        
  337.  
  338.         //detector->Init();
  339.         if (show_3d) detector->Draw(draw_width);
  340.  
  341.         double SiPM = parameters.getA();
  342.         double M = parameters.getM();
  343.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  344.        
  345.         for(int i = 0; i < NN; i++) {
  346.                
  347.                 rx = CENTER.x() - offset;
  348.                 ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  349.                 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  350.                                
  351.                 rphi = rand.Uniform(0.0, 2.0*pi);
  352.                 //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
  353.                 //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  354.                 rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
  355.                
  356.                 rl = TMath::Cos(rtheta);
  357.                 rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
  358.                 rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
  359.                
  360.                 if(show_rand) {
  361.                         htheta->Fill(rtheta);
  362.                         hcostheta->Fill( TMath::Cos(rtheta) );
  363.                         hphi->Fill(rphi);
  364.                         hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
  365.                 }
  366.                
  367.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  368.                 CRay *ray1 = new CRay;
  369.                
  370.                 if(i < show_rays) {
  371.                         detector->Propagate(*ray0, ray1, show_3d);
  372.                         }
  373.                 else {
  374.                         detector->Propagate(*ray0, ray1, 0);
  375.                         }
  376.                
  377.           delete ray0;
  378.           delete ray1; 
  379.         }
  380.        
  381.         if(show_rand) {
  382.                 TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
  383.                 if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
  384.                 else c2rand->Clear();
  385.                 c2rand->Divide(2,3);
  386.                
  387.                 c2rand->cd(1); hphi->Draw();
  388.                 c2rand->cd(3); htheta->Draw();
  389.                 c2rand->cd(5); hcostheta->Draw();
  390.                 c2rand->cd(2); hl->Draw();
  391.                 c2rand->cd(4); hm->Draw();
  392.                 c2rand->cd(6); hn->Draw();
  393.         }
  394.        
  395.         double acceptance = 0.0;
  396.         /*
  397.         if( !(parameters.getPlateOn()) )
  398.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  399.   else
  400.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  401.     */
  402.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  403.         return acceptance;
  404. }
  405.  
  406. // Beamtest distribution
  407. // with fixed theta and phi in interval phiMin, phiMax
  408. double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
  409. {
  410.         double pi = 3.14159265358979312;
  411.         theta *= pi/180.0;
  412.         if(theta < MARGIN) theta = MARGIN;
  413.         phiMin *= pi/180.0;
  414.         phiMax *= pi/180.0;
  415.        
  416.         TDatime now;
  417.         TRandom rand;
  418.         rand.SetSeed(now.Get());
  419.         double rx, ry, rz, rl, rm, rn;
  420.         double rphi;
  421.         //double rtheta;
  422.         //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  423.         //double theta_min_rad = TMath::Cos(theta);
  424.  
  425.         TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
  426.         hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
  427.         hphi = new TH1F("hphi", "hphi", 100, -pi, pi);  
  428.         htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
  429.         htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);       
  430.         hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
  431.         hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); 
  432.         hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
  433.         hl = new TH1F("hl", "hl", 100, 0.0, 1.0);      
  434.         hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
  435.         hm = new TH1F("hm", "hm", 100, 0.0, 1.0);      
  436.         hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
  437.         hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
  438.        
  439.         //detector->Init();
  440.         if (show_3d) detector->Draw(draw_width);
  441.  
  442.         double b = parameters.getB();
  443.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  444.        
  445.         for(int i = 0; i < NN; i++) {
  446.                
  447.                 rx = CENTER.x() - offset;
  448.                 ry = rand.Uniform(-b/2.0, b/2.0);
  449.                 rz = rand.Uniform(-b/2.0, b/2.0);
  450.                                
  451.                 rphi = rand.Uniform(phiMin, phiMax);
  452.                 //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
  453.                 //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  454.                 //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
  455.                
  456.                 rl = TMath::Cos(theta);
  457.                 rm = TMath::Sin(theta)*TMath::Sin(rphi);
  458.                 rn = TMath::Sin(theta)*TMath::Cos(rphi);       
  459.                
  460.                 if(show_rand) {
  461.                         //htheta->Fill(rtheta);
  462.                         //hcostheta->Fill( TMath::Cos(rtheta) );
  463.                         hphi->Fill(rphi);
  464.                         hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
  465.                 }
  466.                
  467.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  468.                 CRay *ray1 = new CRay;
  469.                
  470.                 if(i < show_rays) {
  471.                         detector->Propagate(*ray0, ray1, show_3d);
  472.                         }
  473.                 else {
  474.                         detector->Propagate(*ray0, ray1, 0);
  475.                         }
  476.                
  477.           delete ray0;
  478.           delete ray1; 
  479.         }
  480.        
  481.         if(show_rand) {
  482.                 TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
  483.                 if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
  484.                 else c2rand->Clear();
  485.                 c2rand->Divide(2,3);
  486.                
  487.                 c2rand->cd(1); hphi->Draw();
  488.                 c2rand->cd(3); htheta->Draw();
  489.                 c2rand->cd(5); hcostheta->Draw();
  490.                 c2rand->cd(2); hl->Draw();
  491.                 c2rand->cd(4); hm->Draw();
  492.                 c2rand->cd(6); hn->Draw();
  493.         }
  494.        
  495.         double acceptance = 0.0;
  496.         /*
  497.         if( !(parameters.getPlateOn()) )
  498.           acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  499.   else
  500.     acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
  501.     */
  502.   acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  503.         return acceptance;
  504. }
  505.