96,18 → 96,42 |
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++); |
TH2F* generated = detector->GetGenerated(); |
int 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; |
printf("Statistical variation (max-min)/all = %f perc. \n", variation*100); |
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"); |
(detector->GetHGlass())->Draw("COLZ"); |
TH2F* histoGlass = (detector->GetHGlass()); |
histoGlass->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++); |
TH2F* histoActive = (detector->GetHActive()); |
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++); ((detector->GetLG())->GetHFate())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
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"); |
179,7 → 203,7 |
//----------------------------------------------------------------------------- |
// 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) |
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; |
194,13 → 218,14 |
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"); |
// 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 = new CRay; |
CRay ray1; |
|
detector->Propagate(*ray0, ray1, show_3d); |
|
211,8 → 236,8 |
incidentPolarization->Set(drawPosition, polarization); |
incidentPolarization->DrawS(drawPosition.X(), 1); |
|
TVector3 outPolarization = ray1->GetP(); |
drawPosition = ray1->GetR(); |
TVector3 outPolarization = ray1.GetP(); |
drawPosition = ray1.GetR(); |
drawPosition.SetX(drawPosition.X() + 5); |
CRay* rayPol = new CRay(drawPosition, outPolarization); |
rayPol->SetColor(kBlack); |
219,7 → 244,7 |
rayPol->DrawS(drawPosition.X(), 1); |
|
delete ray0; |
delete ray1; |
//delete ray1; |
|
return (detector->GetHActive())->GetEntries() / (double)(1); |
} |
226,7 → 251,7 |
//----------------------------------------------------------------------------- |
// zarki, razporejeni v mrezi |
double Grid(CDetector *detector, DetectorParameters& parameters, |
int NN = 10, double theta = 0.0) |
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; |
236,27 → 261,32 |
if(show_3d) detector->Draw(draw_width); |
|
const double b = parameters.getB(); |
double simulated = 0; |
|
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; |
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; |
} |
} |
|
} |
267,7 → 297,7 |
else |
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
acceptance = (detector->GetHActive())->GetEntries() / simulated; |
return acceptance; |
} |
//----------------------------------------------------------------------------- |
274,7 → 304,7 |
// 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) |
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; |
283,21 → 313,20 |
if(phi < MARGIN) phi = MARGIN; |
|
TDatime now; |
TRandom rand; |
TRandom2 rand; |
rand.SetSeed(now.Get()); |
|
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
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(-M*SiPM/2.0, M*SiPM/2.0); |
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
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); |
304,17 → 333,18 |
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); |
TVector3 k(ray0->GetK()); |
TVector3 polarization(k.Orthogonal()); |
polarization.Unit(); |
polarization.Rotate(phi, k); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
CRay ray1; |
|
if(i < show_rays) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray1; |
//delete ray1; |
delete ray0; |
} |
|
336,20 → 366,20 |
// 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) |
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; |
theta = theta*3.14159265358979312/180.0; |
if(theta < MARGIN) theta = MARGIN; |
thetaMax = thetaMax*pi/180.0; |
thetaMin = thetaMin*pi/180.0; |
if(thetaMin < MARGIN) thetaMin = MARGIN; |
if(thetaMax < MARGIN) thetaMax = MARGIN; |
|
TDatime now; |
TRandom rand; |
TRandom2 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_min_rad = TMath::Cos(thetaMax); |
double theta_max_rad = TMath::Cos(thetaMin); |
|
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
370,34 → 400,40 |
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
|
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
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(-M*SiPM/2.0, M*SiPM/2.0); |
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
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); |
|
rphi = rand.Uniform(0.0, 2.0*pi); |
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))); |
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad))); |
double theta = 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); |
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(rtheta); |
hcostheta->Fill( TMath::Cos(rtheta) ); |
hphi->Fill(rphi); |
hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
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); |
CRay *ray1 = new CRay; |
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); |
406,8 → 442,8 |
detector->Propagate(*ray0, ray1, 0); |
} |
|
//delete ray1; |
delete ray0; |
delete ray1; |
} |
|
if(show_rand) { |
438,7 → 474,7 |
// 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) |
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand) |
{ |
double pi = 3.14159265358979312; |
theta *= pi/180.0; |
447,10 → 483,10 |
phiMax *= pi/180.0; |
|
TDatime now; |
TRandom rand; |
TRandom2 rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
double rphi; |
double phi; |
|
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
478,19 → 514,19 |
ry = rand.Uniform(-b/2.0, b/2.0); |
rz = rand.Uniform(-b/2.0, b/2.0); |
|
rphi = rand.Uniform(phiMin, phiMax); |
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(rphi); |
rn = TMath::Sin(theta)*TMath::Cos(rphi); |
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(rphi); |
hphi->Fill(phi); |
hl->Fill(rl); |
hm->Fill(rm); |
hn->Fill(rn); |
498,11 → 534,12 |
|
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
// inital polarizaton |
TVector3 polarization(0, 0, 1); |
polarization.RotateY(-theta); |
polarization.RotateX(rphi); |
TVector3 k(ray0->GetK()); |
TVector3 polarization(k.Orthogonal()); |
polarization.Unit(); |
polarization.Rotate(phi, k); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
CRay ray1; |
|
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
511,8 → 548,8 |
detector->Propagate(*ray0, ray1, 0); |
} |
|
//delete ray1; |
delete ray0; |
delete ray1; |
} |
|
if(show_rand) { |