Subversion Repositories f9daq

Rev

Rev 54 | 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(CDetector *detector, double n_in = 1)
  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.         printf("%7.3lf\n", r_acc);
  174. }
  175. //-----------------------------------------------------------------------------
  176.        
  177. //-----------------------------------------------------------------------------
  178. // en zarek
  179. double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
  180. {
  181.   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  182.         theta = Pi()*theta/180.0;
  183.         if(theta < 1e-6) theta = 1e-6;
  184.         phi = phi*Pi()/180.0;
  185.        
  186.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  187.        
  188.         //detector->Init();
  189.         if(show_3d) detector->Draw(draw_width);
  190.        
  191.         CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
  192.                           TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
  193.         CRay *ray1 = new CRay;
  194.        
  195.         detector->Propagate(*ray0, ray1, show_3d);
  196.        
  197.         return (detector->GetHActive())->GetEntries() / (double)(1);
  198. }
  199. //-----------------------------------------------------------------------------
  200. // zarki, razporejeni v mrezi
  201. double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
  202. {
  203. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  204.         theta = Pi()*theta/180.0;
  205.         if(theta < 1e-6) theta = 1e-6;
  206.        
  207.         //detector->Init();
  208.         if(show_3d) detector->Draw(draw_width);
  209.  
  210.         const double b = parameters.getB();
  211.        
  212.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  213.        
  214.         for(int i = 0; i < NN; i++) {
  215.                 for(int j = 0; j < NN; j++) {
  216.                 CRay *ray0 = new CRay(CENTER.x() - offset,
  217.                               0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
  218.                               0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
  219.                                           TMath::Cos(theta),
  220.                                           0.0,
  221.                                           TMath::Sin(theta));
  222.                 CRay *ray1 = new CRay;
  223.                 if(i == (NN/2))
  224.                         detector->Propagate(*ray0, ray1, show_3d);
  225.                 else
  226.                         detector->Propagate(*ray0, ray1, 0);   
  227.                 }
  228.         }
  229.                
  230.         return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
  231. }
  232. //-----------------------------------------------------------------------------
  233. // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
  234. // vsi pod kotom (theta, phi)
  235. double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
  236. {
  237. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  238.         theta = theta*3.14159265358979312/180.0;
  239.         if(theta < MARGIN) theta = MARGIN;
  240.         phi = phi*3.14159265358979312/180.0;
  241.         if(phi < MARGIN) phi = MARGIN;
  242.        
  243.         TDatime now;
  244.         TRandom rand;
  245.         rand.SetSeed(now.Get());
  246.        
  247.         //detector->Init();
  248.         if(show_3d) detector->Draw(draw_width);
  249.        
  250.         double SiPM = parameters.getA();
  251.         double    M = parameters.getM();
  252.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  253.        
  254.         for(int i = 0; i < NN; i++) {
  255.                
  256.                 double rx = CENTER.x() - offset;
  257.                 double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  258.                 double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  259.                                
  260.                 double rl = TMath::Cos(theta);
  261.                 double rm = TMath::Sin(theta)*TMath::Sin(phi);
  262.                 double rn = TMath::Sin(theta)*TMath::Cos(phi); 
  263.                
  264.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  265.                 CRay *ray1 = new CRay;
  266.                
  267.                 if(i < show_rays)
  268.                         detector->Propagate(*ray0, ray1, show_3d);
  269.                 else
  270.                         detector->Propagate(*ray0, ray1, 0);
  271.                 delete ray0;
  272.                 delete ray1;   
  273.         }
  274.                
  275.         return (detector->GetHActive())->GetEntries() / (double)NN;
  276. }
  277. //-----------------------------------------------------------------------------
  278. // zarki, izotropno porazdeljeni znotraj kota theta
  279. // = nakljucni vstopni polozaj in kot
  280. double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
  281. {
  282. //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
  283.         double pi = 3.14159265358979312;
  284.         theta = theta*3.14159265358979312/180.0;
  285.         if(theta < MARGIN) theta = MARGIN;
  286.        
  287.         TDatime now;
  288.         TRandom rand;
  289.         rand.SetSeed(now.Get());
  290.         double rx, ry, rz, rl, rm, rn;
  291.         double rphi, rtheta;
  292.         //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  293.         double theta_min_rad = TMath::Cos(theta);
  294.  
  295.           TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
  296.           hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
  297.           hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
  298.           htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
  299.           htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
  300.           hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
  301.           hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
  302.           hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
  303.           hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
  304.           hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
  305.           hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
  306.           hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
  307.           hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
  308.        
  309.  
  310.         //detector->Init();
  311.         if (show_3d) detector->Draw(draw_width);
  312.  
  313.         double SiPM = parameters.getA();
  314.         double M = parameters.getM();
  315.         double offset = (parameters._plateOn ? parameters._plateWidth : 0);
  316.        
  317.         for(int i = 0; i < NN; i++) {
  318.                
  319.                 rx = CENTER.x() - offset;
  320.                 ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  321.                 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  322.                                
  323.                 rphi = rand.Uniform(0.0, 2.0*pi);
  324.                 rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  325.                
  326.                 rl = TMath::Cos(rtheta);
  327.                 rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
  328.                 rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
  329.                
  330.                 if(show_rand) {
  331.                         htheta->Fill(rtheta);
  332.                         hcostheta->Fill( TMath::Cos(rtheta) );
  333.                         hphi->Fill(rphi);
  334.                         hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
  335.                 }
  336.                
  337.                 CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
  338.                 CRay *ray1 = new CRay;
  339.                
  340.                 if(i < show_rays) {
  341.                         detector->Propagate(*ray0, ray1, show_3d);
  342.                         }
  343.                 else {
  344.                         detector->Propagate(*ray0, ray1, 0);
  345.                         }
  346.                
  347.           delete ray0;
  348.           delete ray1; 
  349.         }
  350.        
  351.         if(show_rand) {
  352.                 TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
  353.                 if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
  354.                 else c2rand->Clear();
  355.                 c2rand->Divide(2,3);
  356.                
  357.                 c2rand->cd(1); hphi->Draw();
  358.                 c2rand->cd(3); htheta->Draw();
  359.                 c2rand->cd(5); hcostheta->Draw();
  360.                 c2rand->cd(2); hl->Draw();
  361.                 c2rand->cd(4); hm->Draw();
  362.                 c2rand->cd(6); hn->Draw();
  363.         }
  364.        
  365.        
  366.         double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;       
  367.         return acceptance;
  368. }
  369.  
  370.