Subversion Repositories f9daq

Rev

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

  1. #include "include/guide.h"
  2. #include "src/raySimulator.cpp"
  3.  
  4. #include "TFile.h"
  5. #include "TPolyMarker3D.h"
  6. #include "TCanvas.h"
  7.  
  8. //extern int show_3d;
  9. //extern int show_data;
  10.  
  11. // Set the simulation parameters
  12. void showVisual(int b) {show_3d = b;}
  13. void showData(int b) {show_data = b;}
  14.  
  15. // Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap, badCoupling)
  16. DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false);
  17. // Print the detector parameters
  18. void getParameters();
  19.  
  20. //void SetLGType(int in = 1, int side = 1, int out = 0)
  21. //{detector->SetLGType(in, side, out);}
  22.  
  23. void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46)
  24. { parameters.setGuide(SiPM0, b0, d0);
  25. parameters.setIndices(n1, n2, n3); }
  26.  
  27. void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
  28. { parameters.setGap(x_gap0, y_gap0, z_gap0); }
  29.  
  30. void setGlass(double glassOn, double glassD)
  31. { parameters.setGlass(glassOn, glassD); };
  32.  
  33. void setPlate(int plateOn = 1, double plateWidth = 1)
  34. { parameters.setPlate(plateOn, plateWidth); }
  35.  
  36. void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
  37.  
  38. // Refractive Indices
  39. // n1 - around light guide - air
  40. // n2 - light guide (plate) material - k9 glass
  41. // n3 - the material at the exit - optical grease, epoxy, air, etc.
  42. void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);}
  43.  
  44. int save_ary = 0;
  45. //------------------------------------------------------------------------------------------
  46. void Help()
  47. {
  48.   printf("void SetCenter(x, y, z)\n");
  49.   printf("void SetLGType(in, side, out)\n");
  50.   printf("SURF_DUMMY 0, ");
  51.   printf("SURF_REFRA 1, ");
  52.   printf("SURF_REFLE 2, ");
  53.   printf("SURF_TOTAL 3, ");
  54.   printf("SURF_IMPER 4\n");
  55.   printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n");
  56.   printf("void SetR(R0)\n");
  57.   printf("void SetGlass(glass_on0, glass_d0)\n");
  58.   printf("void SetGap(x_gap0, y_gap0, z_gap0)\n");
  59.  
  60.   printf("LGG(NN, theta)\n");
  61.   printf("LGR(NN, theta, phi)\n");
  62.   printf("LGI(NN, theta)\n");
  63.  
  64.   printf("LGI_gap(NN, min, max, steps, theta)\n");
  65.   printf("LGR_gap(NN, min, max, steps, theta, phi)\n");
  66.   printf("LGR_th(NN, min, max, steps, phi)\n");
  67.   printf("LGG_th(NN, min, max, steps)\n");
  68.  
  69.   printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)");
  70. }
  71. //------------------------------------------------------------------------------------------
  72. int show_in_steps = 0;
  73. void SetShowInSteps(int in = 1)
  74. {                                                                                
  75.   show_in_steps = in;
  76. }
  77. //------------------------------------------------------------------------------------------
  78. int gs_set = 0;
  79. void SetGS()
  80. {                                                                                
  81.   const UInt_t Number = 2;
  82.   Double_t Red[Number]   = { 1.00, 0.00};
  83.   Double_t Green[Number] = { 1.00, 0.00};
  84.   Double_t Blue[Number]  = { 1.00, 0.00};
  85.   Double_t Stops[Number] = { 0.00, 1.00};
  86.  
  87.   Int_t nb=50;
  88.   if(!gs_set) {
  89.       TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
  90.       gs_set = 1;
  91.   }
  92. }
  93. //------------------------------------------------------------------------------------------
  94. void Draw1(CDetector *detector, int gs = 0, double range = 0.0)
  95. {
  96.   if(gs) SetGS();
  97.   TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500);
  98.   gStyle->SetOptTitle(0);
  99.   gStyle->SetOptStat(0);
  100.   cdata->cd(0);
  101.   if(range > 0.1) {
  102.       ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range);
  103.       ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range);
  104.   }
  105.   (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]");
  106.   (detector->GetHLaser())->Draw("COLZ");
  107.   cdata->SaveAs("2dscan.gif");
  108.   //gStyle->SetOptTitle(1);
  109.   //gStyle->SetOptStat(1);
  110. }
  111.  
  112. //------------------------------------------------------------------------------------------
  113.  
  114. // Test function
  115. // creates an optical boundary surface
  116. // shows the propagation of light ray
  117. // and polarization state of the incident ray (green)
  118. // and surface normal (blue)
  119. void PolTest(double theta = 0.0)
  120. {
  121.   int p_type = SURF_REFRA;
  122.   double p_n1 = parameters.getN1();
  123.   double p_n2 = parameters.getN2();
  124.   theta = 3.141593*theta/180.0;
  125.   if(theta < 1e-6) theta = 1e-6;
  126.  
  127.   Init();
  128.  
  129.   double cx = 0.0;
  130.   TVector3 vodnik_edge[4];
  131.   double t = 3.0;
  132.   vodnik_edge[0].SetXYZ(cx, t,-t);
  133.   vodnik_edge[1].SetXYZ(cx, t, t);
  134.   vodnik_edge[2].SetXYZ(cx,-t, t);
  135.   vodnik_edge[3].SetXYZ(cx,-t,-t);
  136.  
  137.   CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
  138.   surf->FlipN();
  139.   surf->SetFresnel(1);
  140.   surf->Draw();
  141.  
  142.   CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
  143.   ray->SetColor(kRed);
  144.   //p_pol = rotatey(p_pol0, -theta);
  145.   TVector3 sE(0,1,0);
  146.   TVector3 pE(0,0,1);
  147.   TVector3 p_pol = pE;
  148.   p_pol.RotateY(-theta);
  149.   //printf("p_pol = "); printv(p_pol);
  150.   ray->setPolarization(p_pol);
  151.  
  152.   ray->DrawS(cx, -5.0);
  153.  
  154.   CRay *out = new CRay;
  155.   out->SetColor(kBlack);
  156.   TVector3 *inters = new TVector3;
  157.   surf->PropagateRay(*ray, out, inters);
  158.   printf(" n1 = %f, n2 = %f\n", p_n1, p_n2);
  159.   //if(fate == 1) out->DrawS(cx, 5.0);
  160.   out->DrawS(cx, 5.0);
  161.  
  162.   CRay *incidentPolarization = new CRay;
  163.   incidentPolarization->SetColor(kGreen);
  164.   incidentPolarization->Set(ray->GetR(), p_pol);
  165.   incidentPolarization->DrawS(cx, 1.0);
  166.   printf(" GREEN: polarization vector\n");
  167.  
  168.   CRay *surfaceNormal = new CRay;
  169.   surfaceNormal->SetColor(kBlue);
  170.   surfaceNormal->Set(ray->GetR(), surf->GetN());
  171.   surfaceNormal->DrawS(cx, 1.0);
  172.   printf(" BLUE: surface normal vector\n");
  173. }
  174.  
  175. void ptt()
  176. {
  177.   for(double th = 0.0; th < 91.0; th += 5.0) {
  178.       printf("%lf ", th);
  179.       PolTest(th);
  180.   }
  181. }
  182. //------------------------------------------------------------------------------------------
  183. // Propagate single ray and show the statistics
  184. void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
  185. {
  186.   CDetector *detector = new CDetector(CENTER, parameters);
  187.   Init();
  188.   double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
  189.   PrintGuideHead();
  190.   PrintGuideStat(izkoristek);
  191.   DrawData(detector, parameters, theta, izkoristek);
  192. }
  193. //------------------------------------------------------------------------------------------
  194. // Propagate NN rays generated as grid and show the statistics
  195. void LGG(int NN=10, double theta=0.0, bool coupling=false)
  196. {
  197.   parameters.setCoupling(coupling);
  198.   CDetector *detector = new CDetector(CENTER, parameters);
  199.   Init();
  200.   double izkoristek = Grid(detector, parameters, NN, theta);
  201.   //printf("izkoristek = %.3lf\n", izkoristek);
  202.   PrintGuideHead();
  203.   PrintGuideStat(izkoristek);
  204.   DrawData(detector, parameters, theta, izkoristek);
  205. }
  206. //------------------------------------------------------------------------------------------
  207. // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
  208. void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
  209. {
  210.   CDetector *detector = new CDetector(CENTER, parameters);
  211.   Init();
  212.   double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
  213.   //printf("izkoristek = %.3lf\n", izkoristek);
  214.   PrintGuideHead(); PrintGuideStat(izkoristek);
  215.   DrawData(detector, parameters, theta, izkoristek);
  216. }
  217. //------------------------------------------------------------------------------------------
  218. // Propagate NN rays isotropically generated in solid angle theta and show the statistics
  219. void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
  220. {
  221.   Init();
  222.   CDetector *detector = new CDetector(CENTER, parameters);
  223.   //CDetector detector = new CDetector();
  224.   double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
  225.   //printf("izkoristek = %.3lf\n", izkoristek);
  226.   PrintGuideHead();
  227.   PrintGuideStat(izkoristek);
  228.   DrawData(detector, parameters, theta, izkoristek);
  229.   //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
  230.   //canvasDetector->cd();
  231.   //detector->Draw();
  232. }
  233.  
  234. void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
  235. {
  236.   Init();
  237.   parameters.setCoupling(coupling);
  238.   CDetector *detector = new CDetector(CENTER, parameters);
  239.  
  240.   const double theta = 18.5;
  241.   double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr);
  242.  
  243.   PrintGuideHead();
  244.   PrintGuideStat(izkoristek);
  245.   DrawData(detector, parameters, 18.5, izkoristek);
  246.  
  247.   //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
  248.   //canvasDetector->cd();
  249.   //detector->Draw();
  250. }
  251.  
  252.  
  253. //------------------------------------------------------------------------------------------
  254. void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
  255. {
  256.   int tmp3d = show_3d, tmpdata = show_data;
  257.   double step[steps], acc[steps];
  258.  
  259.   show_3d = 0; show_data = 0;
  260.  
  261.   //printf(" x_gap | acceptance \n");
  262.   PrintGuideHead();
  263.   for(int i=0; i<=steps; i++) {
  264.       step[i] = min + i*(max - min)/(double)steps;
  265.       parameters.setGap(step[i], 0.0, 0.0);
  266.       CDetector *detector = new CDetector(CENTER, parameters);
  267.       Init();
  268.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  269.       //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
  270.       PrintGuideStat(acc[i]);
  271.   }
  272.  
  273.   //char sbuff[256];
  274.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  275.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  276.   //    (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  277.  
  278.   DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
  279.  
  280.   show_3d = tmp3d; show_data = tmpdata;
  281. }
  282. //------------------------------------------------------------------------------------------
  283.  
  284. void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
  285. {
  286.   //int tmp3d = show_3d, tmpdata = show_data;
  287.   double show_rays = 10;
  288.   double step[steps], acc[steps];
  289.  
  290.   //show_3d = 0; show_data = 0;
  291.  
  292.   //printf(" theta | acceptance \n");
  293.   PrintGuideHead();
  294.   for(int i=0; i<=steps; i++) {
  295.       Init();
  296.       step[i] = min + i*(max - min)/(double)steps;
  297.       CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  298.       acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
  299.       //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
  300.       PrintGuideStat(acc[i]);
  301.       delete detector;
  302.   }
  303.  
  304.   //char sbuff[256];
  305.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
  306.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  307.   //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
  308.  
  309.   DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  310.  
  311.   //show_3d = tmp3d; show_data = tmpdata;
  312. }
  313.  
  314.  
  315. //------------------------------------------------------------------------------------------
  316. //void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
  317. void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
  318. {
  319.   //int tmp3d = show_3d, tmpdata = show_data;
  320.   double step[steps], acc[steps];
  321.  
  322.   //show_3d = 0; show_data = 0;
  323.  
  324.   //printf(" theta | acceptance \n");
  325.   PrintGuideHead();
  326.   for(int i=0; i<=steps; i++) {
  327.       Init();
  328.       step[i] = min + i*(max - min)/(double)steps;
  329.       CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  330.       //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
  331.       acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
  332.       //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
  333.       PrintGuideStat(acc[i]);
  334.       delete detector;
  335.   }
  336.  
  337.   //char sbuff[256];
  338.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
  339.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  340.   //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
  341.  
  342.   DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  343.  
  344.   //show_3d = tmp3d; show_data = tmpdata;
  345. }
  346. //------------------------------------------------------------------------------------------
  347.  
  348. void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
  349. {
  350.   int tmp3d = show_3d, tmpdata = show_data;
  351.   show_3d = 1; show_data = 0;
  352.   //PrintGuideHead();
  353.   TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
  354.   double min = 0.0;
  355.   printf("\nWait, this takes a while ");
  356.   for(int i=0; i<=steps; i++) {
  357.       double theta = min + i*(maxTheta - min) / (double)steps;
  358.       for (int j=0; j<=steps; j++) {
  359.           Init();
  360.           double phi = min + j*(maxPhi - min) / (double)steps;
  361.           CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  362.           double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
  363.           //PrintGuideStat(acc);
  364.           hAcceptance->Fill(i, j, acc);
  365.           //printf("Acc: %f ", acc);
  366.           delete detector;
  367.       }
  368.       printf(".");
  369.   }
  370.   printf("\n");
  371.   //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  372.   show_3d = tmp3d; show_data = tmpdata;
  373.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
  374.   cp->cd();
  375.   hAcceptance->Draw("COLZ");
  376. }
  377.  
  378. // a vs. d
  379. void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0)
  380. {
  381.   //char sbuff[256];
  382.  
  383.   //show_3d = 0; // don't show simulations
  384.   //show_data = 1;
  385.  
  386.   //double d  = detector->GetD();
  387.   const double b = parameters.getB(); // upper side of LG
  388.   //const double SiPM = 3.0; // the length of the detector itself
  389.   //double reflectivity = detector->guide->getR();
  390.  
  391.   TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
  392.  
  393.   // Use the Fresnel eq. instead of fixed reflectivity 96%
  394.   //detector->SetFresnel(1);
  395.  
  396.   //printf("   d |   a  | Acceptance\n");
  397.   getParameters();
  398.   printf("Wait, this takes a while ");
  399.   for(int i=0; i<steps; i++) {
  400.       const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
  401.       for(int j=0; j<steps; j++) {
  402.           Init();
  403.           const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
  404.           //const double M = b/a;
  405.           parameters.setGuide(a, b, d);
  406.           CDetector *detector = new CDetector(CENTER, parameters);
  407.           //detector->guide->setLG(y, M, x);
  408.           //Init(); exclude simulation
  409.           double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
  410.           //double acceptance = Grid(NN, theta);
  411.           //double acceptance = -1.0;
  412.           hAcceptance->Fill(d, a, acceptance);
  413.           //printf("%.2lf | %.2lf | ", d, a);
  414.           //PrintGuideStat(acceptance);
  415.           delete detector; //works fine, 50x50 grid takes ~4MB of RAM
  416.       }
  417.       printf(".");
  418.   }
  419.   printf("\n");
  420.  
  421.  
  422.  
  423.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
  424.   cp->cd();
  425.   //TVirtualPad *pacc = cp->cd(0);
  426.  
  427.   /*
  428.       pacc->SetRightMargin(0.10);
  429.       pacc->SetLeftMargin(0.10);
  430.       pacc->SetTopMargin(0.10);
  431.       pacc->SetBottomMargin(0.10);
  432.    */
  433.  
  434.   TFile *file = new TFile("acceptance.root","RECREATE");
  435.   hAcceptance->Write();
  436.   file->Close();
  437.   //delete file;
  438.  
  439.   //hAcceptance->SetContour(100);
  440.   //gStyle->SetPalette(1,0);
  441.   hAcceptance->SetTitle(";d [mm];a [mm]");
  442.   hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
  443.   hAcceptance->Draw("COLZ");
  444.   char filename[128];
  445.   sprintf(filename,"LGI_ad%d.C", steps);
  446.   cp->SaveAs(filename);
  447.  
  448. }
  449.  
  450. //------------------------------------------------------------------------------------------
  451. // collection efficiency vs length
  452. // a changes, b rests the same, d changes
  453. // still can't use this function, it has hardcoded tan 10 deg
  454. // !very stupid!
  455. // new 3D graph NEEDED
  456. // d vs a vs acceptance
  457. void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
  458. {
  459.   TVector3 gapGuideSiPM(0.3, 0, 0);
  460.   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  461.   int tmp3d = show_3d, tmpdata = show_data;
  462.   double step[steps], acc[steps];
  463.   double sipm_d, M0;
  464.   char sbuff[256];
  465.  
  466.   show_3d = 1; show_data = 1;
  467.  
  468.   M0 = parameters.getM();
  469.  
  470.   PrintGuideHead();
  471.   for(int i=0; i<=steps; i++) {
  472.       step[i] = min + i*(max - min)/(double)steps;
  473.       sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
  474.       parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
  475.       //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  476.       printf("%.1lf | %.2lf | ", step[i], sipm_d);
  477.       Init();
  478.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  479.       PrintGuideStat(acc[i]);
  480.       sprintf(sbuff, "d%2d.gif", i);
  481.       //c3dview->SaveAs(sbuff);
  482.   }
  483.  
  484.  
  485.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  486.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  487.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  488.  
  489.   //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
  490.  
  491.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
  492.   TVirtualPad *pacc = cp->cd(0);
  493.  
  494.   pacc->SetRightMargin(0.05);
  495.   pacc->SetLeftMargin(0.13);
  496.   pacc->SetTopMargin(0.08);
  497.  
  498.   TGraph *gacceptance = new TGraph(steps+1, step, acc);
  499.   gacceptance->SetTitle("; d [mm];Collection efficiency");
  500.   gacceptance->SetMarkerStyle(8);
  501.   gacceptance->SetLineColor(12);
  502.   gacceptance->SetLineWidth(1);
  503.   (gacceptance->GetXaxis())->SetRangeUser(min, max);
  504.   //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
  505.  
  506.   gacceptance->Draw("ACP");
  507.  
  508.  
  509.   show_3d = tmp3d; show_data = tmpdata;
  510. }
  511. //------------------------------------------------------------------------------------------
  512. // Acceptance vs. the length
  513. // The magnification ratio M rests the same all the time
  514. void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
  515. {
  516.   int tmp3d = show_3d, tmpdata = show_data;
  517.   double step[steps], acc[steps];
  518.   const double a = parameters.getA();
  519.   const double magnif = parameters.getM();
  520.  
  521.   //show_3d = show_in_steps; show_data = 1;
  522.  
  523.   // Use Fresnel equations
  524.   //detector->SetFresnel(1);
  525.  
  526.   // Set glass (n=1.5) at the exit
  527.   //detector->SetGlass(1,0);
  528.  
  529.   PrintGuideHead();
  530.   for(int i=0; i<=steps; i++) {
  531.       step[i] = min + i*(max - min)/(double)steps;
  532.       //parameters.setGuide(a, magnif, step[i]);
  533.       parameters.setGuide(a, magnif, step[i]);
  534.       CDetector *detector = new CDetector(CENTER, parameters);
  535.       //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  536.       printf("%.1lf | ", step[i]);
  537.       Init();
  538.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  539.       PrintGuideStat(acc[i]);
  540.       delete detector;
  541.   }
  542.  
  543.   //char sbuff[256];
  544.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  545.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  546.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  547.  
  548.   DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
  549.  
  550.   show_3d = tmp3d; show_data = tmpdata;
  551. }
  552. //------------------------------------------------------------------------------------------
  553. // Magnification optimization. a and d are fixed, M and b change in steps
  554. void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
  555. {
  556.   int tmp3d = show_3d, tmpdata = show_data;
  557.   double step[steps], acc[steps];
  558.  
  559.   const double a = parameters.getA();
  560.   const double d = parameters.getD();
  561.  
  562.   show_3d = 0; show_data = 0;
  563.  
  564.   PrintGuideHead();
  565.   for(int i=0; i<=steps; i++) {
  566.       step[i] = min + i*(max - min)/(double)steps;
  567.       //parameters.setGuide(a, step[i], d);
  568.       parameters.setGuide(a, a*step[i], d);
  569.       CDetector *detector = new CDetector(CENTER, parameters);
  570.       Init();
  571.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  572.       PrintGuideStat(acc[i]);
  573.   }
  574.  
  575.   //char sbuff[256];
  576.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  577.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  578.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  579.  
  580.   DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
  581.  
  582.   show_3d = tmp3d; show_data = tmpdata;
  583. }
  584.  
  585. //------------------------------------------------------------------------------------------
  586.  
  587. void getParameters()
  588. {
  589.   printf("LIGHT GUIDE\n"
  590.       "   b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
  591.   printf("MATERIAL REFRACITVE INDICES\n"
  592.       "   n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
  593.   printf("PLATE\n"
  594.       "   ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
  595.   return;
  596. }
  597.  
  598.