Subversion Repositories f9daq

Rev

Rev 70 | 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 Show3D(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)
  16. DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.52, TVector3(0.3, 0, 0));
  17.  
  18. //void SetLGType(int in = 1, int side = 1, int out = 0)
  19.         //{detector->SetLGType(in, side, out);}
  20.  
  21. void SetLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n10 = 1.0, double n20 = 1.53, double n30 = 1.52)
  22.         { parameters.setGuide(SiPM0, b0, d0, n10, n20, n30); }
  23.        
  24. void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
  25.   { parameters.setGap(x_gap0, y_gap0, z_gap0); }
  26.  
  27. void SetGlass(double glassOn, double glassD)
  28.   { parameters.setGlass(glassOn, glassD); };
  29.  
  30. void SetPlate(int plateOn = 1, double plateWidth = 1)
  31.   { parameters.setPlate(plateOn, plateWidth); }
  32.  
  33. void SetFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
  34.  
  35. int save_ary = 0;
  36. //------------------------------------------------------------------------------------------
  37. void Help()
  38. {
  39.         printf("void SetCenter(x, y, z)\n");
  40.         printf("void SetLGType(in, side, out)\n");
  41.         printf("SURF_DUMMY 0, ");
  42.         printf("SURF_REFRA 1, ");
  43.         printf("SURF_REFLE 2, ");
  44.         printf("SURF_TOTAL 3, ");
  45.         printf("SURF_IMPER 4\n");
  46.         printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n");
  47.         printf("void SetR(R0)\n");
  48.         printf("void SetGlass(glass_on0, glass_d0)\n");
  49.         printf("void SetGap(x_gap0, y_gap0, z_gap0)\n");
  50.        
  51.         printf("LGG(NN, theta)\n");
  52.         printf("LGR(NN, theta, phi)\n");
  53.         printf("LGI(NN, theta)\n");
  54.        
  55.         printf("LGI_gap(NN, min, max, steps, theta)\n");
  56.         printf("LGR_gap(NN, min, max, steps, theta, phi)\n");
  57.         printf("LGR_th(NN, min, max, steps, phi)\n");
  58.         printf("LGG_th(NN, min, max, steps)\n");
  59.        
  60.         printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)");  
  61. }
  62. //------------------------------------------------------------------------------------------
  63. int show_in_steps = 0;
  64. void SetShowInSteps(int in = 1)
  65. {                                                                                
  66.         show_in_steps = in;
  67. }
  68. //------------------------------------------------------------------------------------------
  69. int gs_set = 0;
  70. void SetGS()
  71. {                                                                                
  72.         const UInt_t Number = 2;
  73.         Double_t Red[Number]   = { 1.00, 0.00};
  74.         Double_t Green[Number] = { 1.00, 0.00};
  75.         Double_t Blue[Number]  = { 1.00, 0.00};
  76.         Double_t Stops[Number] = { 0.00, 1.00};
  77.                                                                                
  78.         Int_t nb=50;
  79.         if(!gs_set) {
  80.                 TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
  81.                 gs_set = 1;
  82.         }
  83. }
  84. //------------------------------------------------------------------------------------------
  85. void Draw1(CDetector *detector, int gs = 0, double range = 0.0)
  86. {
  87.         if(gs) SetGS();
  88.         TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500);
  89.         gStyle->SetOptTitle(0);
  90.         gStyle->SetOptStat(0);
  91.         cdata->cd(0);
  92.                 if(range > 0.1) {
  93.                         ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range);
  94.                         ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range);
  95.                 }
  96.                 (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]");
  97.                 (detector->GetHLaser())->Draw("COLZ");
  98.         cdata->SaveAs("2dscan.gif");
  99.         //gStyle->SetOptTitle(1);
  100.         //gStyle->SetOptStat(1);
  101. }
  102.  
  103. //------------------------------------------------------------------------------------------
  104.  
  105.  
  106. TVector3 p_pol0(0.0, 0.0, 1.0);
  107. void SetPol(double x, double y, double z)
  108. {p_pol0.SetXYZ(x,y,z);}
  109.  
  110. void PolTest(double theta = 0.0)
  111. {
  112.   int p_type = 1;
  113.   double p_n1 = 1.5;
  114.   double p_n2 = 1.0;
  115.         theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
  116.         TVector3 p_pol;
  117.        
  118.         Init();
  119.        
  120.         double cx = 0.0;
  121.         TVector3 vodnik_edge[4];
  122.         double t = 3.0;
  123.         vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t);
  124.         vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t);
  125.         /*
  126. #define SURF_DUMMY 0
  127. #define SURF_REFRA 1
  128. #define SURF_REFLE 2
  129. #define SURF_TOTAL 3
  130. #define SURF_IMPER 4
  131.         */
  132.         CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
  133.         surf->SetFresnel(0);
  134.         surf->Draw();
  135.        
  136.         CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
  137.         ray->SetColor(1);
  138.         //p_pol = rotatey(p_pol0, -theta);
  139.         p_pol = p_pol0; p_pol.RotateY(-theta);
  140.         //printf("p_pol = "); printv(p_pol);
  141.         ray->SetPolarization(p_pol);
  142.        
  143.         ray->DrawS(cx, -5.0);
  144.        
  145.         CRay *out = new CRay; out->SetColor(2);
  146.         TVector3 *inters = new TVector3;
  147.         int fate = surf->PropagateRay(*ray, out, inters);
  148.         if(fate == 1) out->DrawS(cx, 5.0);
  149.        
  150.         CRay *pp = new CRay; pp->SetColor(3);
  151.         pp->Set(ray->GetR(), p_pol);
  152.         pp->DrawS(cx, 2.0);
  153.        
  154.         CRay *pn = new CRay; pn->SetColor(4);
  155.         pn->Set(ray->GetR(), surf->GetN());
  156.         pn->DrawS(cx, 3.0);
  157. }
  158.  
  159. void ptt()
  160. {
  161.         for(double th = 0.0; th < 91.0; th += 5.0) {
  162.                 printf("%lf ", th);
  163.                 PolTest(th);
  164.         }
  165. }
  166. //------------------------------------------------------------------------------------------
  167. // Propagate single ray and show the statistics
  168. void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
  169. {
  170.   CDetector *detector = new CDetector(CENTER, parameters);
  171.         Init();
  172.         double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
  173.         PrintGuideHead();
  174.         PrintGuideStat(detector, 1);
  175.         DrawData(detector, parameters, theta, izkoristek);
  176. }
  177. //------------------------------------------------------------------------------------------
  178. // Propagate NN rays generated as grid and show the statistics
  179. void LGG(int NN = 10, double theta = 0.0)
  180. {
  181.   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  182.         Init();
  183.         double izkoristek = Grid(detector, parameters, NN, theta);
  184.         //printf("izkoristek = %.3lf\n", izkoristek);
  185.         PrintGuideHead();
  186.         PrintGuideStat(detector, (NN+1)*(NN+1));
  187.         DrawData(detector, parameters, theta, izkoristek);
  188. }
  189. //------------------------------------------------------------------------------------------
  190. // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
  191. void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
  192. {
  193.   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  194.         Init();
  195.         double izkoristek = RandYZ(detector, parameters, NN, theta, phi);
  196.         //printf("izkoristek = %.3lf\n", izkoristek);
  197.         PrintGuideHead(); PrintGuideStat(detector, NN);
  198.         DrawData(detector, parameters, theta, izkoristek);
  199. }
  200. //------------------------------------------------------------------------------------------
  201. // Propagate NN rays isotropically generated in solid angle theta and show the statistics
  202. void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
  203. {
  204.   Init();
  205.   CDetector *detector = new CDetector(CENTER, parameters);
  206.   //CDetector detector = new CDetector();
  207.         double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr);
  208.         //printf("izkoristek = %.3lf\n", izkoristek);
  209.         PrintGuideHead(); PrintGuideStat(detector, NN);
  210.         DrawData(detector, parameters, theta, izkoristek);
  211.         //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
  212.         //canvasDetector->cd();
  213.         //detector->Draw();
  214.         //delete detector;
  215. }
  216.  
  217.        
  218.        
  219. //------------------------------------------------------------------------------------------
  220. void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
  221. {
  222. int tmp3d = show_3d, tmpdata = show_data;
  223. double step[steps], acc[steps];
  224.  
  225.         show_3d = 0; show_data = 0;
  226.        
  227.         //printf(" x_gap | acceptance \n");
  228.         PrintGuideHead();
  229.         for(int i=0; i<=steps; i++) {
  230.                 step[i] = min + i*(max - min)/(double)steps;
  231.                 parameters.setGap(step[i], 0.0, 0.0);
  232.                 CDetector *detector = new CDetector(CENTER, parameters);
  233.                 Init();
  234.                 acc[i] = RandIso(detector, parameters, NN, theta);
  235.                 //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
  236.                 PrintGuideStat(detector, NN);
  237.         }
  238.        
  239.         //char sbuff[256];
  240.         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  241.                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
  242.                 //      (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  243.        
  244.         DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
  245.        
  246.         show_3d = tmp3d; show_data = tmpdata;
  247. }
  248. //------------------------------------------------------------------------------------------
  249.  
  250. //------------------------------------------------------------------------------------------
  251. void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 5.0)
  252. {
  253. int tmp3d = show_3d, tmpdata = show_data;
  254. double step[steps], acc[steps];
  255.  
  256.         show_3d = 0; show_data = 0;
  257.        
  258.         //printf(" theta | acceptance \n");
  259.         PrintGuideHead();
  260.         for(int i=0; i<=steps; i++) {
  261.                 Init();
  262.                 step[i] = min + i*(max - min)/(double)steps;
  263.                 CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  264.                 acc[i] = RandYZ(detector, parameters, NN, step[i], phi);
  265.                 //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
  266.                 PrintGuideStat(detector, NN);
  267.                 delete detector;
  268.         }
  269.        
  270.         //char sbuff[256];
  271.         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
  272.                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
  273.                 //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
  274.        
  275.         DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
  276.        
  277.         show_3d = tmp3d; show_data = tmpdata;
  278. }
  279. //------------------------------------------------------------------------------------------
  280.  
  281. // a vs. d
  282. 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)
  283. {
  284.   int tmp3d = show_3d;
  285.         int tmpdata = show_data;
  286.         //char sbuff[256];
  287.  
  288.         show_3d = 0; // don't show simulations
  289.         show_data = 1;
  290.        
  291.         //double d  = detector->GetD();
  292.         const double b = parameters.getB(); // upper side of LG
  293.         //const double SiPM = 3.0; // the length of the detector itself
  294.         //double reflectivity = detector->guide->getR();
  295.        
  296.         TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
  297.        
  298.         // Use the Fresnel eq. instead of fixed reflectivity 96%
  299.         //detector->SetFresnel(1);
  300.        
  301.         printf("   d |   a  | Acceptance\n");
  302.         for(int i=0; i<steps; i++) {
  303.                 const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
  304.                 for(int j=0; j<steps; j++) {
  305.                   const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
  306.                   const double M = b/a;
  307.       parameters.setGuide(a, M, d);
  308.                   CDetector *detector = new CDetector(CENTER, parameters);
  309.                   //detector->guide->setLG(y, M, x);
  310.                   //Init(); exclude simulation
  311.                   double acceptance = RandIso(detector, parameters, NN, theta, 0, 0);
  312.                   //double acceptance = Grid(NN, theta);
  313.                   //double acceptance = -1.0;
  314.                   hAcceptance->Fill(d, a, acceptance);
  315.                   //printf("%.2lf | %.2lf | ", x,y);
  316.                   //PrintGuideStat(detector, NN);
  317.                   //delete detector; //works fine, 50x50 grid takes ~4MB of RAM
  318.                   }
  319.         }
  320.        
  321.  
  322.        
  323.         TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
  324.         TVirtualPad *pacc = cp->cd(0);
  325.        
  326.         pacc->SetRightMargin(0.10);
  327.         pacc->SetLeftMargin(0.10);
  328.         pacc->SetTopMargin(0.10);
  329.        
  330.         TFile *file = new TFile("acceptance.root","RECREATE");
  331.         hAcceptance->Write();
  332.         file->Close();
  333.         delete file;
  334.        
  335.         //hAcceptance->SetContour(100);
  336.         //gStyle->SetPalette(1,0);
  337.         hAcceptance->SetTitle(";d [mm];a [mm]");
  338.         hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
  339.         hAcceptance->Draw("COLZ");
  340.        
  341.         show_3d = tmp3d; show_data = tmpdata;
  342. }
  343.  
  344. //------------------------------------------------------------------------------------------
  345. // collection efficiency vs length
  346. // a changes, b rests the same, d changes
  347. // still can't use this function, it has hardcoded tan 10 deg
  348. // !very stupid!
  349. // new 3D graph NEEDED
  350. // d vs a vs acceptance
  351. void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
  352. {
  353.   TVector3 gapGuideSiPM(0.3, 0, 0);
  354.   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
  355.         int tmp3d = show_3d, tmpdata = show_data;
  356.         double step[steps], acc[steps];
  357.         double sipm_d, M0;
  358.         char sbuff[256];
  359.  
  360.         show_3d = 1; show_data = 1;
  361.        
  362.         M0 = parameters.getM();
  363.        
  364.         PrintGuideHead();
  365.         for(int i=0; i<=steps; i++) {
  366.                 step[i] = min + i*(max - min)/(double)steps;           
  367.                 sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
  368.                 parameters.setGuide(sipm_d, M0/sipm_d, step[i]);
  369.                 //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  370.                 printf("%.1lf | %.2lf | ", step[i], sipm_d);
  371.                 Init();
  372.                 acc[i] = RandIso(detector, parameters, NN, theta);
  373.                 PrintGuideStat(detector, NN);
  374.                 sprintf(sbuff, "d%2d.gif", i);
  375.                 //c3dview->SaveAs(sbuff);
  376.         }
  377.        
  378.        
  379.         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  380.                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
  381.                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  382.        
  383.         //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
  384.        
  385.         TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
  386.         TVirtualPad *pacc = cp->cd(0);
  387.        
  388.         pacc->SetRightMargin(0.05);
  389.         pacc->SetLeftMargin(0.13);
  390.         pacc->SetTopMargin(0.08);
  391.        
  392.         TGraph *gacceptance = new TGraph(steps+1, step, acc);
  393.         gacceptance->SetTitle("; d [mm];Collection efficiency");
  394.         gacceptance->SetMarkerStyle(8);
  395.         gacceptance->SetLineColor(12);
  396.         gacceptance->SetLineWidth(1);
  397.         (gacceptance->GetXaxis())->SetRangeUser(min, max);
  398.         //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
  399.        
  400.         gacceptance->Draw("ACP");
  401.        
  402.        
  403.         show_3d = tmp3d; show_data = tmpdata;
  404. }
  405. //------------------------------------------------------------------------------------------
  406. // Acceptance vs. the length
  407. // The magnification ratio M rests the same all the time
  408. void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
  409. {
  410.         int tmp3d = show_3d, tmpdata = show_data;
  411.         double step[steps], acc[steps];
  412.         const double a = parameters.getA();
  413.         const double magnif = parameters.getM();
  414.  
  415.         //show_3d = show_in_steps; show_data = 1;
  416.        
  417.         // Use Fresnel equations
  418.         //detector->SetFresnel(1);
  419.        
  420.         // Set glass (n=1.5) at the exit
  421.         //detector->SetGlass(1,0);
  422.        
  423.         PrintGuideHead();
  424.         for(int i=0; i<=steps; i++) {
  425.                 step[i] = min + i*(max - min)/(double)steps;
  426.                 parameters.setGuide(a, magnif, step[i]);
  427.                 CDetector *detector = new CDetector(CENTER, parameters);
  428.                 //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
  429.                 printf("%.1lf | ", step[i]);
  430.                 Init();
  431.                 acc[i] = RandIso(detector, parameters, NN, theta);
  432.                 PrintGuideStat(detector, NN);
  433.                 delete detector;
  434.         }
  435.        
  436.         //char sbuff[256];
  437.         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  438.                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
  439.                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  440.        
  441.         DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
  442.        
  443.         show_3d = tmp3d; show_data = tmpdata;
  444. }
  445. //------------------------------------------------------------------------------------------
  446. // Magnification optimization. a and d are fixed, M and b change in steps
  447. void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
  448. {
  449.         int tmp3d = show_3d, tmpdata = show_data;
  450.         double step[steps], acc[steps];
  451.        
  452.         const double a = parameters.getA();
  453.         const double d = parameters.getD();
  454.  
  455.         show_3d = 0; show_data = 0;
  456.        
  457.         PrintGuideHead();
  458.         for(int i=0; i<=steps; i++) {
  459.                 step[i] = min + i*(max - min)/(double)steps;
  460.                 parameters.setGuide(a, step[i], d);
  461.                 CDetector *detector = new CDetector(CENTER, parameters);
  462.                 Init();
  463.                 acc[i] = RandIso(detector, parameters, NN, theta);
  464.                 PrintGuideStat(detector, NN);
  465.         }
  466.        
  467.         //char sbuff[256];
  468.         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
  469.                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
  470.                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
  471.        
  472.         DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
  473.        
  474.         show_3d = tmp3d; show_data = tmpdata;
  475. }
  476. //------------------------------------------------------------------------------------------
  477.  
  478.