Subversion Repositories f9daq

Rev

Rev 73 | 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 = 4.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;
  138.   //surf = CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
  139.   surf.Set(p_type, vodnik_edge, p_n1, p_n2, 0.96);
  140.   surf.FlipN();
  141.   surf.SetFresnel(1);
  142.   surf.Draw();
  143.  
  144.   CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
  145.   ray->SetColor(kRed);
  146.   //p_pol = rotatey(p_pol0, -theta);
  147.   TVector3 sE(0,1,0);
  148.   TVector3 pE(0,0,1);
  149.   TVector3 p_pol = pE;
  150.   p_pol.RotateY(-theta);
  151.   //printf("p_pol = "); printv(p_pol);
  152.   ray->setPolarization(p_pol);
  153.  
  154.   ray->DrawS(cx, -5.0);
  155.  
  156.   CRay out;
  157.   out.SetColor(kBlack);
  158.   TVector3 *inters = new TVector3;
  159.   surf.PropagateRay(*ray, out, inters);
  160.   printf(" n1 = %f, n2 = %f\n", p_n1, p_n2);
  161.   //if(fate == 1) out->DrawS(cx, 5.0);
  162.   out.DrawS(cx, 5.0);
  163.  
  164.   CRay *incidentPolarization = new CRay;
  165.   incidentPolarization->SetColor(kGreen);
  166.   incidentPolarization->Set(ray->GetR(), p_pol);
  167.   incidentPolarization->DrawS(cx, 1.0);
  168.   printf(" GREEN: polarization vector\n");
  169.  
  170.   CRay *surfaceNormal = new CRay;
  171.   surfaceNormal->SetColor(kBlue);
  172.   surfaceNormal->Set(ray->GetR(), surf.GetN());
  173.   surfaceNormal->DrawS(cx, 1.0);
  174.   printf(" BLUE: surface normal vector\n");
  175. }
  176.  
  177. void ptt()
  178. {
  179.   for(double th = 0.0; th < 91.0; th += 5.0) {
  180.       printf("%lf ", th);
  181.       PolTest(th);
  182.   }
  183. }
  184. //------------------------------------------------------------------------------------------
  185. // Propagate single ray and show the statistics
  186. void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
  187. {
  188.   CDetector *detector = new CDetector(CENTER, parameters);
  189.   Init();
  190.   double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
  191.   PrintGuideHead();
  192.   PrintGuideStat(izkoristek);
  193.   DrawData(detector, parameters, theta, izkoristek);
  194. }
  195. //------------------------------------------------------------------------------------------
  196. // Propagate NN rays generated as grid and show the statistics
  197. void LGG(int NN=10, double theta=0.0, bool coupling=false)
  198. {
  199.   parameters.setCoupling(coupling);
  200.   CDetector *detector = new CDetector(CENTER, parameters);
  201.   Init();
  202.   double izkoristek = Grid(detector, parameters, NN, theta);
  203.   //printf("izkoristek = %.3lf\n", izkoristek);
  204.   PrintGuideHead();
  205.   PrintGuideStat(izkoristek);
  206.   DrawData(detector, parameters, theta, izkoristek);
  207. }
  208. //------------------------------------------------------------------------------------------
  209. // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
  210. void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
  211. {
  212.   CDetector *detector = new CDetector(CENTER, parameters);
  213.   Init();
  214.   double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
  215.   //printf("izkoristek = %.3lf\n", izkoristek);
  216.   PrintGuideHead(); PrintGuideStat(izkoristek);
  217.   DrawData(detector, parameters, theta, izkoristek);
  218. }
  219. //------------------------------------------------------------------------------------------
  220. // Propagate NN rays isotropically generated in solid angle theta and show the statistics
  221. void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
  222. {
  223.   Init();
  224.   CDetector *detector = new CDetector(CENTER, parameters);
  225.   //CDetector detector = new CDetector();
  226.   double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
  227.   //printf("izkoristek = %.3lf\n", izkoristek);
  228.   PrintGuideHead();
  229.   PrintGuideStat(izkoristek);
  230.   DrawData(detector, parameters, theta, izkoristek);
  231.   //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
  232.   //canvasDetector->cd();
  233.   //detector->Draw();
  234. }
  235.  
  236. void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
  237. {
  238.   Init();
  239.   parameters.setCoupling(coupling);
  240.   CDetector *detector = new CDetector(CENTER, parameters);
  241.  
  242.   const double theta = 18.5;
  243.   double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr);
  244.  
  245.   PrintGuideHead();
  246.   PrintGuideStat(izkoristek);
  247.   DrawData(detector, parameters, 18.5, izkoristek);
  248.  
  249.   //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
  250.   //canvasDetector->cd();
  251.   //detector->Draw();
  252. }
  253.  
  254.  
  255. //------------------------------------------------------------------------------------------
  256. void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
  257. {
  258.   int tmp3d = show_3d, tmpdata = show_data;
  259.   double step[steps], acc[steps];
  260.  
  261.   show_3d = 0; show_data = 0;
  262.  
  263.   //printf(" x_gap | acceptance \n");
  264.   PrintGuideHead();
  265.   for(int i=0; i<=steps; i++) {
  266.       step[i] = min + i*(max - min)/(double)steps;
  267.       parameters.setGap(step[i], 0.0, 0.0);
  268.       CDetector *detector = new CDetector(CENTER, parameters);
  269.       Init();
  270.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  271.       //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
  272.       PrintGuideStat(acc[i]);
  273.   }
  274.  
  275.   //char sbuff[256];
  276.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  277.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  278.   //    (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  279.  
  280.   DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
  281.  
  282.   show_3d = tmp3d; show_data = tmpdata;
  283. }
  284. //------------------------------------------------------------------------------------------
  285.  
  286. void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
  287. {
  288.   //int tmp3d = show_3d, tmpdata = show_data;
  289.   double show_rays = 10;
  290.   double step[steps], acc[steps];
  291.  
  292.   //show_3d = 0; show_data = 0;
  293.  
  294.   //printf(" theta | acceptance \n");
  295.   PrintGuideHead();
  296.   for(int i=0; i<=steps; i++) {
  297.       Init();
  298.       step[i] = min + i*(max - min)/(double)steps;
  299.       CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  300.       acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
  301.       //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
  302.       PrintGuideStat(acc[i]);
  303.       delete detector;
  304.   }
  305.  
  306.   //char sbuff[256];
  307.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
  308.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  309.   //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
  310.  
  311.   DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  312.  
  313.   //show_3d = tmp3d; show_data = tmpdata;
  314. }
  315.  
  316.  
  317. //------------------------------------------------------------------------------------------
  318. //void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
  319. void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
  320. {
  321.   //int tmp3d = show_3d, tmpdata = show_data;
  322.   double step[steps], acc[steps];
  323.  
  324.   //show_3d = 0; show_data = 0;
  325.  
  326.   //printf(" theta | acceptance \n");
  327.   PrintGuideHead();
  328.   for(int i=0; i<=steps; i++) {
  329.       Init();
  330.       step[i] = min + i*(max - min)/(double)steps;
  331.       CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  332.       //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
  333.       acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
  334.       //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
  335.       PrintGuideStat(acc[i]);
  336.       delete detector;
  337.   }
  338.  
  339.   //char sbuff[256];
  340.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
  341.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  342.   //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
  343.  
  344.   DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  345.  
  346.   //show_3d = tmp3d; show_data = tmpdata;
  347. }
  348. //------------------------------------------------------------------------------------------
  349.  
  350. void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
  351. {
  352.   int tmp3d = show_3d, tmpdata = show_data;
  353.   show_3d = 1; show_data = 0;
  354.   //PrintGuideHead();
  355.   TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
  356.   double min = 0.0;
  357.   printf("\nWait, this takes a while ");
  358.   for(int i=0; i<=steps; i++) {
  359.       double theta = min + i*(maxTheta - min) / (double)steps;
  360.       for (int j=0; j<=steps; j++) {
  361.           Init();
  362.           double phi = min + j*(maxPhi - min) / (double)steps;
  363.           CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  364.           double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
  365.           //PrintGuideStat(acc);
  366.           hAcceptance->Fill(i, j, acc);
  367.           //printf("Acc: %f ", acc);
  368.           delete detector;
  369.       }
  370.       printf(".");
  371.   }
  372.   printf("\n");
  373.   //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  374.   show_3d = tmp3d; show_data = tmpdata;
  375.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
  376.   cp->cd();
  377.   hAcceptance->Draw("COLZ");
  378. }
  379.  
  380. // a vs. d
  381. 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)
  382. {
  383.   //char sbuff[256];
  384.  
  385.   //show_3d = 0; // don't show simulations
  386.   //show_data = 1;
  387.  
  388.   //double d  = detector->GetD();
  389.   const double b = parameters.getB(); // upper side of LG
  390.   //const double SiPM = 3.0; // the length of the detector itself
  391.   //double reflectivity = detector->guide->getR();
  392.  
  393.   TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
  394.  
  395.   // Use the Fresnel eq. instead of fixed reflectivity 96%
  396.   //detector->SetFresnel(1);
  397.  
  398.   //printf("   d |   a  | Acceptance\n");
  399.   getParameters();
  400.   printf("Wait, this takes a while ");
  401.   for(int i=0; i<steps; i++) {
  402.       const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
  403.       for(int j=0; j<steps; j++) {
  404.           Init();
  405.           const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
  406.           //const double M = b/a;
  407.           parameters.setGuide(a, b, d);
  408.           CDetector *detector = new CDetector(CENTER, parameters);
  409.           //detector->guide->setLG(y, M, x);
  410.           //Init(); exclude simulation
  411.           double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
  412.           //double acceptance = Grid(NN, theta);
  413.           //double acceptance = -1.0;
  414.           hAcceptance->Fill(d, a, acceptance);
  415.           //printf("%.2lf | %.2lf | ", d, a);
  416.           //PrintGuideStat(acceptance);
  417.           delete detector; //works fine, 50x50 grid takes ~4MB of RAM
  418.       }
  419.       printf(".");
  420.   }
  421.   printf("\n");
  422.  
  423.  
  424.  
  425.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
  426.   cp->cd();
  427.   //TVirtualPad *pacc = cp->cd(0);
  428.  
  429.   /*
  430.       pacc->SetRightMargin(0.10);
  431.       pacc->SetLeftMargin(0.10);
  432.       pacc->SetTopMargin(0.10);
  433.       pacc->SetBottomMargin(0.10);
  434.    */
  435.  
  436.   TFile *file = new TFile("acceptance.root","RECREATE");
  437.   hAcceptance->Write();
  438.   file->Close();
  439.   //delete file;
  440.  
  441.   //hAcceptance->SetContour(100);
  442.   //gStyle->SetPalette(1,0);
  443.   hAcceptance->SetTitle(";d [mm];a [mm]");
  444.   hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
  445.   hAcceptance->Draw("COLZ");
  446.   char filename[128];
  447.   sprintf(filename,"LGI_ad%d.C", steps);
  448.   cp->SaveAs(filename);
  449.  
  450. }
  451.  
  452. //------------------------------------------------------------------------------------------
  453. // collection efficiency vs length
  454. // a changes, b rests the same, d changes
  455. // still can't use this function, it has hardcoded tan 10 deg
  456. // !very stupid!
  457. // new 3D graph NEEDED
  458. // d vs a vs acceptance
  459. void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
  460. {
  461.   TVector3 gapGuideSiPM(0.3, 0, 0);
  462.   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  463.   int tmp3d = show_3d, tmpdata = show_data;
  464.   double step[steps], acc[steps];
  465.   double sipm_d, M0;
  466.   char sbuff[256];
  467.  
  468.   show_3d = 1; show_data = 1;
  469.  
  470.   M0 = parameters.getM();
  471.  
  472.   PrintGuideHead();
  473.   for(int i=0; i<=steps; i++) {
  474.       step[i] = min + i*(max - min)/(double)steps;
  475.       sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
  476.       parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
  477.       //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  478.       printf("%.1lf | %.2lf | ", step[i], sipm_d);
  479.       Init();
  480.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  481.       PrintGuideStat(acc[i]);
  482.       sprintf(sbuff, "d%2d.gif", i);
  483.       //c3dview->SaveAs(sbuff);
  484.   }
  485.  
  486.  
  487.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  488.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  489.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  490.  
  491.   //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
  492.  
  493.   TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
  494.   TVirtualPad *pacc = cp->cd(0);
  495.  
  496.   pacc->SetRightMargin(0.05);
  497.   pacc->SetLeftMargin(0.13);
  498.   pacc->SetTopMargin(0.08);
  499.  
  500.   TGraph *gacceptance = new TGraph(steps+1, step, acc);
  501.   gacceptance->SetTitle("; d [mm];Collection efficiency");
  502.   gacceptance->SetMarkerStyle(8);
  503.   gacceptance->SetLineColor(12);
  504.   gacceptance->SetLineWidth(1);
  505.   (gacceptance->GetXaxis())->SetRangeUser(min, max);
  506.   //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
  507.  
  508.   gacceptance->Draw("ACP");
  509.  
  510.  
  511.   show_3d = tmp3d; show_data = tmpdata;
  512. }
  513. //------------------------------------------------------------------------------------------
  514. // Acceptance vs. the length
  515. // The magnification ratio M rests the same all the time
  516. void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
  517. {
  518.   int tmp3d = show_3d, tmpdata = show_data;
  519.   double step[steps], acc[steps];
  520.   const double a = parameters.getA();
  521.   const double magnif = parameters.getM();
  522.  
  523.   //show_3d = show_in_steps; show_data = 1;
  524.  
  525.   // Use Fresnel equations
  526.   //detector->SetFresnel(1);
  527.  
  528.   // Set glass (n=1.5) at the exit
  529.   //detector->SetGlass(1,0);
  530.  
  531.   PrintGuideHead();
  532.   for(int i=0; i<=steps; i++) {
  533.       step[i] = min + i*(max - min)/(double)steps;
  534.       //parameters.setGuide(a, magnif, step[i]);
  535.       parameters.setGuide(a, magnif, step[i]);
  536.       CDetector *detector = new CDetector(CENTER, parameters);
  537.       //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  538.       printf("%.1lf | ", step[i]);
  539.       Init();
  540.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  541.       PrintGuideStat(acc[i]);
  542.       delete detector;
  543.   }
  544.  
  545.   //char sbuff[256];
  546.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  547.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  548.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  549.  
  550.   DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
  551.  
  552.   show_3d = tmp3d; show_data = tmpdata;
  553. }
  554. //------------------------------------------------------------------------------------------
  555. // Magnification optimization. a and d are fixed, M and b change in steps
  556. void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
  557. {
  558.   int tmp3d = show_3d, tmpdata = show_data;
  559.   double step[steps], acc[steps];
  560.  
  561.   const double a = parameters.getA();
  562.   const double d = parameters.getD();
  563.  
  564.   show_3d = 0; show_data = 0;
  565.  
  566.   PrintGuideHead();
  567.   for(int i=0; i<=steps; i++) {
  568.       step[i] = min + i*(max - min)/(double)steps;
  569.       //parameters.setGuide(a, step[i], d);
  570.       parameters.setGuide(a, a*step[i], d);
  571.       CDetector *detector = new CDetector(CENTER, parameters);
  572.       Init();
  573.       acc[i] = RandIso(detector, parameters, NN, 0, theta);
  574.       PrintGuideStat(acc[i]);
  575.   }
  576.  
  577.   //char sbuff[256];
  578.   //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  579.   //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
  580.   //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  581.  
  582.   DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
  583.  
  584.   show_3d = tmp3d; show_data = tmpdata;
  585. }
  586.  
  587. //------------------------------------------------------------------------------------------
  588.  
  589. void getParameters()
  590. {
  591.   printf("LIGHT GUIDE\n"
  592.       "   b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
  593.   printf("MATERIAL REFRACITVE INDICES\n"
  594.       "   n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
  595.   printf("PLATE\n"
  596.       "   ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
  597.   return;
  598. }
  599.  
  600.