Subversion Repositories f9daq

Rev

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