Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1. #include "RTUtil.h"
  2. #include "CGuide.h"
  3.  
  4. #include "TROOT.h"
  5. #include "TSystem.h"
  6. #include "TStyle.h"
  7. #include "TCanvas.h"
  8. #include "TMath.h"
  9. #include "TView3D.h"
  10. #include "TH1F.h"
  11. #include "TH2F.h"
  12. #include "TBenchmark.h"
  13. #include "TPolyMarker.h"
  14. #include "TGraph.h"
  15. #include "TF1.h"
  16.  
  17. int show_3d = 1;
  18. int show_data = 1;
  19. int draw_width = 2;
  20. TVector3 center(-2.0, 0.0, 0.0);
  21.  
  22. double RandIso(int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
  23. {
  24.   CDetector *detector = new CDetector(center);
  25.         double pi = 3.14159265358979312;
  26.         theta = theta*3.14159265358979312/180.0;
  27.         if(theta < MARGIN) theta = MARGIN;
  28.        
  29.         TDatime now;
  30.         TRandom rand;
  31.         rand.SetSeed(now.Get());
  32.         double rx, ry, rz, rl, rm, rn;
  33.         double rphi, rtheta;
  34.         //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
  35.         double theta_min_rad = TMath::Cos(theta);
  36.         //if (show_rand) {
  37.           TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
  38.           hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
  39.           hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
  40.           htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
  41.           htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
  42.           hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
  43.           hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
  44.           hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
  45.           hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
  46.           hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
  47.           hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
  48.           hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
  49.           hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
  50.         //}
  51.        
  52.  
  53.         detector->Init();
  54.         if (show_3d) detector->Draw(draw_width);
  55.        
  56.         int fate;
  57.         double SiPM = detector->GetSiPM();
  58.         double M = detector->GetM();
  59.         CRay *ray0 = new CRay;
  60.         CRay *ray1 = new CRay;
  61.        
  62.         for(int i = 0; i < NN; i++) {
  63.                
  64.                 rx = center.x();
  65.                 ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  66.                 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
  67.                                
  68.                 rphi = rand.Uniform(0.0, 2.0*pi);
  69.                 rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
  70.                
  71.                 rl = TMath::Cos(rtheta);
  72.                 rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
  73.                 rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
  74.                
  75.                 if(show_rand) {
  76.                         htheta->Fill(rtheta);
  77.                         hcostheta->Fill( TMath::Cos(rtheta) );
  78.                         hphi->Fill(rphi);
  79.                         hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
  80.                 }
  81.                
  82.                 ray0->Set(rx, ry, rz, rl, rm, rn);
  83.                
  84.                 if(i < show_rays)
  85.                         fate = detector->Propagate(*ray0, ray1, show_3d);
  86.                 else
  87.                         fate = detector->Propagate(*ray0, ray1, 0);    
  88.         }
  89.        
  90.         if(show_rand) {
  91.                 TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
  92.                 if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
  93.                 else c2rand->Clear();
  94.                 c2rand->Divide(2,3);
  95.                
  96.                 c2rand->cd(1); hphi->Draw();
  97.                 c2rand->cd(3); htheta->Draw();
  98.                 c2rand->cd(5); hcostheta->Draw();
  99.                 c2rand->cd(2); hl->Draw();
  100.                 c2rand->cd(4); hm->Draw();
  101.                 c2rand->cd(6); hn->Draw();
  102.                 delete c2rand;
  103.         }
  104.        
  105.         delete hphi;
  106.         delete htheta;
  107.         delete hcostheta;
  108.         delete hl;
  109.         delete hm;
  110.         delete hn;
  111.        
  112.         delete ray0;
  113.         delete ray1;
  114.        
  115.                
  116.         return (detector->GetHActive())->GetEntries() / (double)NN;
  117. }
  118.  
  119. // a vs. d
  120. 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)
  121. {
  122.         show_3d = 0; // don't show every simulation
  123.         show_data = 1;
  124.        
  125.         //double d  = detector->GetD();
  126.         double b = 5.0; // upper side of LG
  127.         double R = 0.96;
  128.        
  129.         TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
  130.        
  131.         printf("   d |   a  | Acceptance\n");
  132.         for(int i=0; i<steps; i++) {
  133.                 double x = hAcceptance->GetXaxis()->GetBinCenter(i);
  134.                 for(int j=0; j<steps; j++) {
  135.                   double y = hAcceptance->GetYaxis()->GetBinCenter(j);
  136.                   double M = b/y;
  137.                   // hardcoded n0 and n1
  138.                   SetLG(y, M, x);
  139.                   //Init(); exclude simulation
  140.                   double acceptance = RandIso(NN, theta, 5, 0);
  141.                   //double acceptance = Grid(NN, theta);
  142.                   //double acceptance = -1.0;
  143.                   hAcceptance->Fill(x, y, acceptance);
  144.                   //printf("%.2lf | %.2lf | ", x,y);
  145.                   //PrintGuideStat(NN);
  146.                   }
  147.         }
  148.        
  149.  
  150.        
  151.         TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
  152.         TVirtualPad *pacc = cp->cd(0);
  153.        
  154.         pacc->SetRightMargin(0.10);
  155.         pacc->SetLeftMargin(0.10);
  156.         pacc->SetTopMargin(0.10);
  157.        
  158.         //TFile *file = new TFile("acceptance.root","RECREATE");
  159.         //hAcceptance->Write();
  160.         //file->Close();
  161.        
  162.         //hAcceptance->SetContour(100);
  163.         //gStyle->SetPalette(1,0);
  164.         hAcceptance->SetTitle(";d [mm];a [mm]");
  165.         hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
  166.         hAcceptance->Draw("COLZ");
  167.        
  168.  return;
  169. }
  170.