Rev 84 |
Go to most recent revision |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
// own include
//#include "include/raySimulator.h"
#include "include/RTUtil.h"
#include "include/guide.h"
// general include
#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"
//=================================================================================
TCanvas *c3dview;
TCanvas *cacc;
int show_3d = 0, show_data = 0;
int draw_width = 2;
// calls default constructor CVodnik.cpp
//void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
//{center.SetXYZ(x,y,z); detector = new CDetector(center);}
//void SetLGType(int in = 1, int side = 1, int out = 0)
//{detector->SetLGType(in, side, out);}
//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)
//{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
//void SetR(double R0) {detector->SetR(R0);}
/*
void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)
{detector->SetGlass(glass_on0, glass_d0);}
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
{detector->SetGap(x_gap0 , y_gap0, z_gap0);}
void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)
{detector->SetRCol(in, lg, out, gla);};
void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)
{detector->SetDCol(LG0, glass0, active0);};
void SetDWidth(int w = 1) {draw_width = w;}
void SetFresnel(int b = 1) {detector->SetFresnel(b);}
void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}
void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
*/
//-----------------------------------------------------------------------------
void Init(double rsys_scale = 10.0)
{
RTSetStyle(gStyle);
gStyle->SetOptStat(10);
if(show_3d) {
double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
c3dview = (TCanvas*)gROOT->FindObject("c3dview");
if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
c3dview->SetFillColor(0);}
else c3dview->Clear();
TView3D *view = new TView3D(1, rsys0, rsys1);
view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
//view->SetPerspective();
view->SetParallel();
view->FrontView();
view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
//view->ToggleRulers();
}
}
//-----------------------------------------------------------------------------
void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
{
if(show_data) {
char sbuff[256];
sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
parameters.getA(), parameters.getLightYield()*acc, parameters.getD(),
parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc*100);
if(!only2d) {
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
cdata->Divide(3,3);
int cpc = 1;
cdata->cd(cpc++);
TH2F* generated = detector->GetGenerated();
double nGenerated = generated->GetEntries();
//int minimum = generated->GetBinContent(generated->GetMinimumBin());
int maximum = generated->GetBinContent(generated->GetMaximumBin());
generated->GetZaxis()->SetRangeUser(0, 1.05*maximum);
//double variation = (maximum-minimum)/(double)nGenerated;
generated->SetTitle("Generated");
generated->Draw("colz");
cdata->cd(cpc++); //((detector->GetLG())->GetHIn())->Draw("COLZ");
TH2F* histoPlate = detector->GetHPlate();
histoPlate->Draw("COLZ");
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
TH2F* histoGlass = (detector->GetHGlass());
histoGlass->Draw("COLZ");
cdata->cd(cpc++);
TH2F* histoActive = (detector->GetHActive());
//double nAccepted = histoActive->GetEntries();
//double variation = 1 / sqrt(nAccepted);
double variation = sqrt(acc*(1-acc)/nGenerated);
printf("Statistical error = %f perc.\n", variation*100);
histoActive->Draw("COLZ");
cdata->cd(cpc++);
TH2F* histoLaser = (detector->GetHLaser());
histoLaser->Draw("COLZ");
cdata->cd(cpc++);
TH2F* histoDetector = (detector->GetHDetector());
histoDetector->Draw("COLZ");
cdata->cd(cpc++);
TH1F* fate=((detector->GetLG())->GetHFate());
fate->Draw();
cdata->cd(cpc++);
TH1F* reflections = ((detector->GetLG())->GetHNOdb_all());
reflections->Draw();
cdata->cd(cpc++);
TH1F *reflectedExit = ((detector->GetLG())->GetHNOdb_exit());
reflectedExit->Draw();
}else{
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
cdata->SaveAs("cdata.gif");
}
}
}
//-----------------------------------------------------------------------------
void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
{
cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
(cacc->cd(0))->SetRightMargin(0.05);
(cacc->cd(0))->SetLeftMargin(0.13);
(cacc->cd(0))->SetTopMargin(0.08);
TGraph *gacceptance = new TGraph(n, x, y);
gacceptance->SetTitle(htitle);
gacceptance->SetMarkerStyle(8);
gacceptance->SetLineColor(12);
gacceptance->SetLineWidth(1);
if(xmin < xmax)
(gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
gacceptance->Draw("ACP");
//cacc->Print("acceptance.pdf");
//delete gacceptance;
//delete cacc;
}
//-----------------------------------------------------------------------------
void PrintGlassStat(CDetector *detector)
{
int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
int glass_out = (detector->GetHGlass())->GetEntries();
printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
}
//-----------------------------------------------------------------------------
void PrintGuideHead()
{
//printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
printf("Acceptance\n");
}
//-----------------------------------------------------------------------------
void PrintGuideStat(double acceptance)
{
/*int n_active = (detector->GetHActive())->GetEntries();
int n_enter = (detector->GetLG())->GetEnteranceHits();
double izkoristek = n_active/(double)n_enter;
int fates[7]; (detector->GetLG())->GetVFate(fates);
// printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
(detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);
printf("%5.1lf\n", 100.0*izkoristek);*/
//double n_active = (detector->GetHActive())->GetEntries();
//double r_acc = 100.0 * n_active / n_in;
double r_acc = 100.0 * acceptance;
printf("%7.3lf\n", r_acc);
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
// en zarek
double Single(CDetector *detector, DetectorParameters& parameters,
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
{
theta = Pi()*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
phi = phi*Pi()/180.0;
if(phi < 1e-6) phi = 1e-6;
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
//detector->Init();
if(show_3d) detector->Draw(draw_width);
CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
TMath::Sin(theta)*TMath::Cos(phi));
// Set y-polarization == vertical
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
if (polarization.Dot(ray0->GetK()) > 1e-2) printf("ERROR: pol not perep\n");
ray0->setPolarization(polarization);
CRay ray1;
detector->Propagate(*ray0, ray1, show_3d);
CRay *incidentPolarization = new CRay;
incidentPolarization->SetColor(kGreen);
TVector3 drawPosition = ray0->GetR();
drawPosition.SetX(drawPosition.X() - 5);
incidentPolarization->Set(drawPosition, polarization);
incidentPolarization->DrawS(drawPosition.X(), 1);
TVector3 outPolarization = ray1.GetP();
drawPosition = ray1.GetR();
drawPosition.SetX(drawPosition.X() + 5);
CRay* rayPol = new CRay(drawPosition, outPolarization);
rayPol->SetColor(kBlack);
rayPol->DrawS(drawPosition.X(), 1);
delete ray0;
//delete ray1;
return (detector->GetHActive())->GetEntries() / (double)(1);
}
//-----------------------------------------------------------------------------
// zarki, razporejeni v mrezi
double Grid(CDetector *detector, DetectorParameters& parameters,
int NN = 10, double theta = 0.0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = Pi()*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
//detector->Init();
if(show_3d) detector->Draw(draw_width);
const double b = parameters.getB();
double simulated = 0;
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
const int GRID = 100;
for(int i = 0; i < GRID; i++) {
for(int j = 0; j < GRID; j++) {
for(int jk = 0; jk<NN; jk++){
CRay *ray0 = new CRay(CENTER.x() - offset,
(i*b/GRID + b/(2.0*GRID) - b/2.0),
(j*b/GRID + b/(2.0*GRID) - b/2.0),
TMath::Cos(theta),
0.0,
TMath::Sin(theta));
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
ray0->setPolarization(polarization);
CRay ray1;
if(i == (GRID/2) and jk == 0)
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
delete ray0;
simulated++;
//delete ray1;
}
}
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / simulated;
return acceptance;
}
//-----------------------------------------------------------------------------
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phi, int show_rays)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = theta*3.14159265358979312/180.0;
if(theta < MARGIN) theta = MARGIN;
phi = phi*3.14159265358979312/180.0;
if(phi < MARGIN) phi = MARGIN;
TDatime now;
TRandom2 rand;
rand.SetSeed(now.Get());
//detector->Init();
if(show_3d) detector->Draw(draw_width);
double b = parameters.getB();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
for(int i = 0; i < NN; i++) {
double rx = CENTER.x() - offset;
double ry = rand.Uniform(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
double rl = TMath::Cos(theta);
double rm = TMath::Sin(theta)*TMath::Sin(phi);
double rn = TMath::Sin(theta)*TMath::Cos(phi);
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
ray0->setPolarization(polarization);
CRay ray1;
if(i < show_rays)
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
//delete ray1;
delete ray0;
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
//return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/
acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
return acceptance;
}
//-----------------------------------------------------------------------------
// zarki, izotropno porazdeljeni znotraj kota theta
// = nakljucni vstopni polozaj in kot
// NN - number of rays to be simulated with angles theta distributed uniformly
// in cos theta and phi uniformly:
// theta [0, theta]
// phi [0,360]
double RandIso(CDetector *detector, DetectorParameters& parameters,
int NN = 1e3, double thetaMin=0.0, double thetaMax = 30.0, int show_rays = 30, int show_rand = 0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
double pi = 3.14159265358979312;
thetaMax = thetaMax*pi/180.0;
thetaMin = thetaMin*pi/180.0;
if(thetaMin < MARGIN) thetaMin = MARGIN;
if(thetaMax < MARGIN) thetaMax = MARGIN;
TDatime now;
TRandom2 rand;
rand.SetSeed(now.Get());
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
double theta_min_rad = TMath::Cos(thetaMax);
double theta_max_rad = TMath::Cos(thetaMin);
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);
double b = parameters.getB();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
for(int i = 0; i < NN; i++) {
double rx = CENTER.x() - offset;
double ry = rand.Uniform(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
double phi = rand.Uniform(0.0, 2.0*pi);
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
double rl = TMath::Cos(theta);
double rm = TMath::Sin(theta)*TMath::Sin(phi);
double rn = TMath::Sin(theta)*TMath::Cos(phi);
if(show_rand) {
htheta->Fill(theta);
hcostheta->Fill( TMath::Cos(theta) );
hphi->Fill(phi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
}
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
ray0->setPolarization(polarization);
CRay ray1;
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
}
else {
detector->Propagate(*ray0, ray1, 0);
}
//delete ray1;
delete ray0;
}
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();
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
return acceptance;
}
// Beamtest distribution
// with fixed theta and phi in interval phiMin, phiMax
double beamtest(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
{
double pi = 3.14159265358979312;
theta *= pi/180.0;
if(theta < MARGIN) theta = MARGIN;
phiMin *= pi/180.0;
phiMax *= pi/180.0;
TDatime now;
TRandom2 rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double phi;
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
hphi = new TH1F("hphi", "hphi", 100, -pi, 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);
double b = parameters.getB();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
for(int i = 0; i < NN; i++) {
rx = CENTER.x() - offset;
ry = rand.Uniform(-b/2.0, b/2.0);
rz = rand.Uniform(-b/2.0, b/2.0);
phi = rand.Uniform(phiMin, phiMax);
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
rl = TMath::Cos(theta);
rm = TMath::Sin(theta)*TMath::Sin(phi);
rn = TMath::Sin(theta)*TMath::Cos(phi);
if(show_rand) {
//htheta->Fill(rtheta);
//hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(phi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
}
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
// inital polarizaton
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
ray0->setPolarization(polarization);
CRay ray1;
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
}
else {
detector->Propagate(*ray0, ray1, 0);
}
//delete ray1;
delete ray0;
}
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();
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
return acceptance;
}