23,32 → 23,32 |
int show_3d = 0, show_data = 0; |
int draw_width = 2; |
|
// calls default constructor CVodnik.cpp |
// 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);} |
//{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);} |
//{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);} |
|
//{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);} |
58,244 → 58,275 |
|
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);} |
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);} |
*/ |
*/ |
//----------------------------------------------------------------------------- |
void Init(double rsys_scale = 7.0) |
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(); |
} |
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); |
|
if(!only2d) { |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
cdata->Divide(3,3); |
int cpc = 1; |
cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ"); |
(detector->GetHGlass())->Draw("COLZ"); |
|
cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
|
cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
}else{ |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500); |
cdata->cd(0); (detector->GetHLaser())->Draw("COLZ"); |
cdata->SaveAs("cdata.gif"); |
} |
} |
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); |
|
if(!only2d) { |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
cdata->Divide(3,3); |
int cpc = 1; |
cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ"); |
(detector->GetHGlass())->Draw("COLZ"); |
|
cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
|
cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->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) |
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; |
|
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); |
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"); |
//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_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); |
|
//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) |
double Single(CDetector *detector, DetectorParameters& parameters, |
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 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; |
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)); |
CRay *ray1 = new CRay; |
|
detector->Propagate(*ray0, ray1, show_3d); |
|
delete ray0; |
delete ray1; |
|
return (detector->GetHActive())->GetEntries() / (double)(1); |
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 z-polarization == vertical |
TVector3 polarization(0,0,1); |
polarization.RotateY(-theta); |
polarization.RotateZ(phi); |
if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n"); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
|
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) |
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); |
//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; |
|
const double b = parameters.getB(); |
|
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
for(int j = 0; j < NN; j++) { |
CRay *ray0 = new CRay(CENTER.x() - offset, |
0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
TMath::Cos(theta), |
0.0, |
TMath::Sin(theta)); |
CRay *ray1 = new CRay; |
if(i == (NN/2)) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray0; |
delete ray1; |
} |
|
} |
double acceptance = 0.0; |
/* |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
|
const double b = parameters.getB(); |
|
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
for(int j = 0; j < NN; j++) { |
CRay *ray0 = new CRay(CENTER.x() - offset, |
0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
TMath::Cos(theta), |
0.0, |
TMath::Sin(theta)); |
TVector3 polarization(0, 1, 0); |
polarization.RotateY(-theta); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
if(i == (NN/2)) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray0; |
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() / (double)((NN)*(NN)); |
return acceptance; |
*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
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) |
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; |
TRandom rand; |
rand.SetSeed(now.Get()); |
|
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
|
double rx = CENTER.x() - offset; |
double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/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); |
CRay *ray1 = new CRay; |
|
if(i < show_rays) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
//delete ray0; |
//delete ray1; |
} |
|
double acceptance = 0.0; |
/* |
//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; |
TRandom rand; |
rand.SetSeed(now.Get()); |
|
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
|
double rx = CENTER.x() - offset; |
double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/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 polarization(0, 0, 1); |
polarization.RotateY(-theta); |
polarization.RotateZ(phi); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
|
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; |
acceptance = (detector->GetHActive())->GetEntries() / (double) NN; |
return acceptance; |
} |
//----------------------------------------------------------------------------- |
// zarki, izotropno porazdeljeni znotraj kota theta |
304,201 → 335,207 |
// 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 theta = 30.0, int show_rays = 30, int show_rand = 0) |
double RandIso(CDetector *detector, DetectorParameters& parameters, |
int NN = 1e3, double thetaMin=0.0, double theta = 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; |
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); |
double theta_max_rad = TMath::Cos(thetaMin); |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
double pi = 3.14159265358979312; |
theta = theta*3.14159265358979312/180.0; |
if(theta < MARGIN) theta = MARGIN; |
|
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); |
|
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); |
double theta_max_rad = TMath::Cos(thetaMin); |
|
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
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); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
|
rx = CENTER.x() - offset; |
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); |
//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(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); |
} |
|
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
CRay *ray1 = new CRay; |
|
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
} |
else { |
detector->Propagate(*ray0, ray1, 0); |
} |
|
delete ray0; |
delete ray1; |
} |
|
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; |
/* |
|
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
|
for(int i = 0; i < NN; i++) { |
|
rx = CENTER.x() - offset; |
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); |
//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(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); |
} |
|
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
CRay *ray1 = new CRay; |
|
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
} |
else { |
detector->Propagate(*ray0, ray1, 0); |
} |
|
delete ray0; |
delete ray1; |
} |
|
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; |
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 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; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
double rphi; |
//double rtheta; |
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
//double theta_min_rad = TMath::Cos(theta); |
double pi = 3.14159265358979312; |
theta *= pi/180.0; |
if(theta < MARGIN) theta = MARGIN; |
phiMin *= pi/180.0; |
phiMax *= pi/180.0; |
|
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); |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
double rphi; |
|
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); |
|
rphi = 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(rphi); |
rn = TMath::Sin(theta)*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); |
} |
|
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
CRay *ray1 = new CRay; |
|
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
} |
else { |
detector->Propagate(*ray0, ray1, 0); |
} |
|
delete ray0; |
delete ray1; |
} |
|
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; |
/* |
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); |
|
rphi = 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(rphi); |
rn = TMath::Sin(theta)*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); |
} |
|
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
// inital polarizaton |
TVector3 polarization(0, 0, 1); |
polarization.RotateY(-theta); |
polarization.RotateX(rphi); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
|
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
} |
else { |
detector->Propagate(*ray0, ray1, 0); |
} |
|
delete ray0; |
delete ray1; |
} |
|
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; |
return acceptance; |
} |