Go to most recent revision |
Blame |
Last modification |
View Log
| RSS feed
#include "RTUtil.h"
#include "CGuide.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TView3D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TBenchmark.h"
#include "TPolyMarker.h"
#include "TGraph.h"
#include "TF1.h"
int show_3d = 1;
int show_data = 1;
int draw_width = 2;
TVector3 center(-2.0, 0.0, 0.0);
double RandIso(int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
{
CDetector *detector = new CDetector(center);
double pi = 3.14159265358979312;
theta = theta*3.14159265358979312/180.0;
if(theta < MARGIN) theta = MARGIN;
TDatime now;
TRandom rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double rphi, rtheta;
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
double theta_min_rad = TMath::Cos(theta);
//if (show_rand) {
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
//}
detector->Init();
if (show_3d) detector->Draw(draw_width);
int fate;
double SiPM = detector->GetSiPM();
double M = detector->GetM();
CRay *ray0 = new CRay;
CRay *ray1 = new CRay;
for(int i = 0; i < NN; i++) {
rx = center.x();
ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
rphi = rand.Uniform(0.0, 2.0*pi);
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
rl = TMath::Cos(rtheta);
rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
if(show_rand) {
htheta->Fill(rtheta);
hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(rphi);
hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
}
ray0->Set(rx, ry, rz, rl, rm, rn);
if(i < show_rays)
fate = detector->Propagate(*ray0, ray1, show_3d);
else
fate = detector->Propagate(*ray0, ray1, 0);
}
if(show_rand) {
TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
else c2rand->Clear();
c2rand->Divide(2,3);
c2rand->cd(1); hphi->Draw();
c2rand->cd(3); htheta->Draw();
c2rand->cd(5); hcostheta->Draw();
c2rand->cd(2); hl->Draw();
c2rand->cd(4); hm->Draw();
c2rand->cd(6); hn->Draw();
delete c2rand;
}
delete hphi;
delete htheta;
delete hcostheta;
delete hl;
delete hm;
delete hn;
delete ray0;
delete ray1;
return (detector->GetHActive())->GetEntries() / (double)NN;
}
// a vs. d
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)
{
show_3d = 0; // don't show every simulation
show_data = 1;
//double d = detector->GetD();
double b = 5.0; // upper side of LG
double R = 0.96;
TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
printf(" d | a | Acceptance\n");
for(int i=0; i<steps; i++) {
double x = hAcceptance->GetXaxis()->GetBinCenter(i);
for(int j=0; j<steps; j++) {
double y = hAcceptance->GetYaxis()->GetBinCenter(j);
double M = b/y;
// hardcoded n0 and n1
SetLG(y, M, x);
//Init(); exclude simulation
double acceptance = RandIso(NN, theta, 5, 0);
//double acceptance = Grid(NN, theta);
//double acceptance = -1.0;
hAcceptance->Fill(x, y, acceptance);
//printf("%.2lf | %.2lf | ", x,y);
//PrintGuideStat(NN);
}
}
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
TVirtualPad *pacc = cp->cd(0);
pacc->SetRightMargin(0.10);
pacc->SetLeftMargin(0.10);
pacc->SetTopMargin(0.10);
//TFile *file = new TFile("acceptance.root","RECREATE");
//hAcceptance->Write();
//file->Close();
//hAcceptance->SetContour(100);
//gStyle->SetPalette(1,0);
hAcceptance->SetTitle(";d [mm];a [mm]");
hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
hAcceptance->Draw("COLZ");
return;
}