Subversion Repositories f9daq

Rev

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

  1. // own include
  2. //#include "include/raySimulator.h"
  3. #include "include/RTUtil.h"
  4. #include "include/guide.h"
  5.  
  6. // general include
  7. #include "TROOT.h"
  8. #include "TSystem.h"
  9. #include "TStyle.h"
  10. #include "TCanvas.h"
  11. #include "TMath.h"
  12. #include "TView3D.h"
  13. #include "TH1F.h"
  14. #include "TH2F.h"
  15. #include "TBenchmark.h"
  16. #include "TPolyMarker.h"
  17. #include "TGraph.h"
  18. #include "TF1.h"
  19.  
  20. //=================================================================================
  21. TCanvas *c3dview;
  22. TCanvas *cacc;
  23. int show_3d = 0, show_data = 0;
  24. int draw_width = 2;
  25.  
  26.  // calls default constructor CVodnik.cpp
  27.  
  28. //void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
  29.         //{center.SetXYZ(x,y,z); detector = new CDetector(center);}
  30.  
  31.  
  32. //void SetLGType(int in = 1, int side = 1, int out = 0)
  33.         //{detector->SetLGType(in, side, out);}
  34.  
  35. //void SetLG(double SiPM0 = 3.0, double M0 = 1.666, double d0 = 3.0, double n10 = 1.0, double n20 = 1.48, double n30 = 1.48, double R0 = 0.96)
  36.         //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
  37.        
  38. //void SetR(double R0) {detector->SetR(R0);}
  39.         /*
  40. void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)
  41.         {detector->SetGlass(glass_on0, glass_d0);}
  42.  
  43. void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
  44.         {detector->SetGap(x_gap0 , y_gap0, z_gap0);}
  45.        
  46. void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)
  47.         {detector->SetRCol(in, lg, out, gla);};
  48. void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)
  49.         {detector->SetDCol(LG0, glass0, active0);};
  50.        
  51.  
  52. void SetDWidth(int w = 1) {draw_width = w;}
  53.  
  54. void SetFresnel(int b = 1) {detector->SetFresnel(b);}
  55. void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}
  56.  
  57. void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}
  58.  
  59. void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
  60. void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
  61. */
  62. //-----------------------------------------------------------------------------
  63. void Init(double rsys_scale = 7.0)
  64. {      
  65.         RTSetStyle(gStyle);
  66.         gStyle->SetOptStat(10);
  67.        
  68.         if(show_3d) {
  69.                 double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
  70.                 double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
  71.                
  72.                 c3dview = (TCanvas*)gROOT->FindObject("c3dview");
  73.                 if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);      c3dview->SetFillColor(0);}
  74.                 else c3dview->Clear();
  75.                        
  76.                 TView3D *view = new TView3D(1, rsys0, rsys1);
  77.                 view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
  78.                 //view->SetPerspective();
  79.                 view->SetParallel();
  80.                 view->FrontView();
  81.                 view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
  82.                 //view->ToggleRulers();
  83.         }
  84. }
  85. //-----------------------------------------------------------------------------
  86. void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
  87. {
  88.         if(show_data) {
  89.                 char sbuff[256];
  90.                 sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
  91.                                 parameters.getA(), parameters.getM(), parameters.getD(),
  92.                                 parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
  93.                
  94.                 if(!only2d) {
  95.                         RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
  96.                         cdata->Divide(3,3);
  97.                         int cpc = 1;
  98.                         cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
  99.                         cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
  100.                         cdata->cd(cpc++); ((detector->GetLG())->GetHOut())->Draw("COLZ");
  101.                        
  102.                        
  103.                         cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
  104.                         cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
  105.                         cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
  106.                        
  107.                         cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
  108.                         cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
  109.                         cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
  110.                 }else{         
  111.                         RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
  112.                         cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
  113.                         cdata->SaveAs("cdata.gif");
  114.                 }
  115.         }
  116. }
  117. //-----------------------------------------------------------------------------
  118. void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
  119.                          double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
  120. {
  121.                
  122.         cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
  123.        
  124.         (cacc->cd(0))->SetRightMargin(0.05);
  125.         (cacc->cd(0))->SetLeftMargin(0.13);
  126.         (cacc->cd(0))->SetTopMargin(0.08);
  127.        
  128.         TGraph *gacceptance = new TGraph(n, x, y);
  129.         gacceptance->SetTitle(htitle);
  130.         gacceptance->SetMarkerStyle(8);
  131.         gacceptance->SetLineColor(12);
  132.         gacceptance->SetLineWidth(1);
  133.         if(xmin < xmax)
  134.                 (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
  135.         (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
  136.        
  137.         gacceptance->Draw("ACP");
  138.         //cacc->Print("acceptance.pdf");
  139.         //delete gacceptance;
  140.         //delete cacc;
  141. }
  142. //-----------------------------------------------------------------------------
  143. void PrintGlassStat(CDetector *detector)
  144. {
  145.         int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
  146.         int glass_out = (detector->GetHGlass())->GetEntries();
  147.         printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
  148. }
  149. //-----------------------------------------------------------------------------
  150. void PrintGuideHead()
  151. {
  152.         //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  153.         printf("Acceptance\n");
  154. }
  155. //-----------------------------------------------------------------------------
  156. void PrintGuideStat(double acceptance)
  157. {
  158.         /*int n_active = (detector->GetHActive())->GetEntries();
  159.         int n_enter = (detector->GetLG())->GetEnteranceHits();
  160.         double izkoristek = n_active/(double)n_enter;
  161.         int fates[7]; (detector->GetLG())->GetVFate(fates);
  162.        
  163. //      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  164.        
  165.         printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
  166.                    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
  167.         for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       
  168.        
  169.         printf("%5.1lf\n", 100.0*izkoristek);*/
  170.        
  171.         //double n_active = (detector->GetHActive())->GetEntries();
  172.         //double r_acc = 100.0 * n_active / n_in;
  173.         double r_acc = 100.0 * acceptance;
  174.         printf("%7.3lf\n", r_acc);
  175. }
  176. //-----------------------------------------------------------------------------
  177.        
  178. //-----------------------------------------------------------------------------
  179. // en zarek
  180. double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
  181. {
  182.   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  183.         theta = Pi()*theta/180.0;
  184.         if(theta < 1e-6) theta = 1e-6;
  185.         phi = phi*Pi()/180.0;
  186.        
  187.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  188.        
  189.         //detector->Init();
  190.         if(show_3d) detector->Draw(draw_width);
  191.        
  192.         CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
  193.                           TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
  194.         CRay *ray1 = new CRay;
  195.        
  196.         detector->Propagate(*ray0, ray1, show_3d);
  197.        
  198.         return (detector->GetHActive())->GetEntries() / (double)(1);
  199. }
  200. //-----------------------------------------------------------------------------
  201. // zarki, razporejeni v mrezi
  202. double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
  203. {
  204. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  205.         theta = Pi()*theta/180.0;
  206.         if(theta < 1e-6) theta = 1e-6;
  207.        
  208.         //detector->Init();
  209.         if(show_3d) detector->Draw(draw_width);
  210.  
  211.         const double b = parameters.getB();
  212.        
  213.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  214.        
  215.         for(int i = 0; i < NN; i++) {
  216.                 for(int j = 0; j < NN; j++) {
  217.                 CRay *ray0 = new CRay(CENTER.x() - offset,
  218.                               0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
  219.                               0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
  220.                                           TMath::Cos(theta),
  221.                                           0.0,
  222.                                           TMath::Sin(theta));
  223.                 CRay *ray1 = new CRay;
  224.                 if(i == (NN/2))
  225.                         detector->Propagate(*ray0, ray1, show_3d);
  226.                 else
  227.                         detector->Propagate(*ray0, ray1, 0);   
  228.                 }
  229.         }
  230.                
  231.         return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
  232. }
  233. //-----------------------------------------------------------------------------
  234. // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
  235. // vsi pod kotom (theta, phi)
  236. double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
  237. {
  238. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  239.         theta = theta*3.14159265358979312/180.0;
  240.         if(theta < MARGIN) theta = MARGIN;
  241.         phi = phi*3.14159265358979312/180.0;
  242.         if(phi < MARGIN) phi = MARGIN;
  243.        
  244.         TDatime now;
  245.         TRandom rand;
  246.         rand.SetSeed(now.Get());
  247.        
  248.         //detector->Init();
  249.         if(show_3d) detector->Draw(draw_width);
  250.        
  251.         double SiPM = parameters.getA();
  252.         double    M = parameters.getM();
  253.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  254.        
  255.         for(int i = 0; i < NN; i++) {
  256.                
  257.                 double rx = CENTER.x() - offset;
  258.                 double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  259.                 double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  260.                                
  261.                 double rl = TMath::Cos(theta);
  262.                 double rm = TMath::Sin(theta)*TMath::Sin(phi);
  263.                 double rn = TMath::Sin(theta)*TMath::Cos(phi); 
  264.                
  265.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  266.                 CRay *ray1 = new CRay;
  267.                
  268.                 if(i < show_rays)
  269.                         detector->Propagate(*ray0, ray1, show_3d);
  270.                 else
  271.                         detector->Propagate(*ray0, ray1, 0);
  272.                 //delete ray0;
  273.                 //delete ray1; 
  274.         }
  275.                
  276.         return (detector->GetHActive())->GetEntries() / (double)NN;
  277. }
  278. //-----------------------------------------------------------------------------
  279. // zarki, izotropno porazdeljeni znotraj kota theta
  280. // = nakljucni vstopni polozaj in kot
  281. double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
  282. {
  283. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  284.         double pi = 3.14159265358979312;
  285.         theta = theta*3.14159265358979312/180.0;
  286.         if(theta < MARGIN) theta = MARGIN;
  287.        
  288.         TDatime now;
  289.         TRandom rand;
  290.         rand.SetSeed(now.Get());
  291.         double rx, ry, rz, rl, rm, rn;
  292.         double rphi, rtheta;
  293.         //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  294.         double theta_min_rad = TMath::Cos(theta);
  295.  
  296.           TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
  297.           hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
  298.           hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
  299.           htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
  300.           htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
  301.           hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
  302.           hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
  303.           hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
  304.           hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
  305.           hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
  306.           hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
  307.           hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
  308.           hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
  309.        
  310.  
  311.         //detector->Init();
  312.         if (show_3d) detector->Draw(draw_width);
  313.  
  314.         double SiPM = parameters.getA();
  315.         double M = parameters.getM();
  316.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  317.        
  318.         for(int i = 0; i < NN; i++) {
  319.                
  320.                 rx = CENTER.x() - offset;
  321.                 ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  322.                 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  323.                                
  324.                 rphi = rand.Uniform(0.0, 2.0*pi);
  325.                 rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  326.                
  327.                 rl = TMath::Cos(rtheta);
  328.                 rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
  329.                 rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
  330.                
  331.                 if(show_rand) {
  332.                         htheta->Fill(rtheta);
  333.                         hcostheta->Fill( TMath::Cos(rtheta) );
  334.                         hphi->Fill(rphi);
  335.                         hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
  336.                 }
  337.                
  338.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  339.                 CRay *ray1 = new CRay;
  340.                
  341.                 if(i < show_rays) {
  342.                         detector->Propagate(*ray0, ray1, show_3d);
  343.                         }
  344.                 else {
  345.                         detector->Propagate(*ray0, ray1, 0);
  346.                         }
  347.                
  348.           //delete ray0;
  349.           //delete ray1;       
  350.         }
  351.        
  352.         if(show_rand) {
  353.                 TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
  354.                 if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
  355.                 else c2rand->Clear();
  356.                 c2rand->Divide(2,3);
  357.                
  358.                 c2rand->cd(1); hphi->Draw();
  359.                 c2rand->cd(3); htheta->Draw();
  360.                 c2rand->cd(5); hcostheta->Draw();
  361.                 c2rand->cd(2); hl->Draw();
  362.                 c2rand->cd(4); hm->Draw();
  363.                 c2rand->cd(6); hn->Draw();
  364.         }
  365.        
  366.        
  367.         double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;       
  368.         return acceptance;
  369. }
  370.  
  371.