Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 84 → Rev 83

/lightguide/trunk/src/raySimulator.cpp
96,42 → 96,18
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();
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->histoPlate)->Draw("COLZ");
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
TH2F* histoGlass = (detector->GetHGlass());
histoGlass->Draw("COLZ");
(detector->GetHGlass())->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->GetHActive())->Draw("COLZ");
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
cdata->cd(cpc++); (detector->GetHDetector())->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();
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");
203,7 → 179,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;
218,14 → 194,13
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");
// 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;
CRay *ray1 = new CRay;
 
detector->Propagate(*ray0, ray1, show_3d);
 
236,8 → 211,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);
244,7 → 219,7
rayPol->DrawS(drawPosition.X(), 1);
 
delete ray0;
//delete ray1;
delete ray1;
 
return (detector->GetHActive())->GetEntries() / (double)(1);
}
251,7 → 226,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;
261,32 → 236,27
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;
}
 
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;
}
 
}
297,7 → 267,7
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / simulated;
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
return acceptance;
}
//-----------------------------------------------------------------------------
304,7 → 274,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;
313,20 → 283,21
if(phi < MARGIN) phi = MARGIN;
 
TDatime now;
TRandom2 rand;
TRandom rand;
rand.SetSeed(now.Get());
 
//detector->Init();
if(show_3d) detector->Draw(draw_width);
 
double b = parameters.getB();
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(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
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);
333,18 → 304,17
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);
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateZ(phi);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays)
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
//delete ray1;
delete ray1;
delete ray0;
}
 
366,20 → 336,20
// 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)
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;
thetaMax = thetaMax*pi/180.0;
thetaMin = thetaMin*pi/180.0;
if(thetaMin < MARGIN) thetaMin = MARGIN;
if(thetaMax < MARGIN) thetaMax = MARGIN;
theta = theta*3.14159265358979312/180.0;
if(theta < MARGIN) theta = MARGIN;
 
TDatime now;
TRandom2 rand;
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(thetaMax);
double theta_min_rad = TMath::Cos(theta);
double theta_max_rad = TMath::Cos(thetaMin);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
400,40 → 370,34
//detector->Init();
if (show_3d) detector->Draw(draw_width);
 
double b = parameters.getB();
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(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
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 phi = rand.Uniform(0.0, 2.0*pi);
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)));
double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
rtheta = 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);
rl = TMath::Cos(rtheta);
rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
 
if(show_rand) {
htheta->Fill(theta);
hcostheta->Fill( TMath::Cos(theta) );
hphi->Fill(phi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
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);
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
442,8 → 406,8
detector->Propagate(*ray0, ray1, 0);
}
 
//delete ray1;
delete ray0;
delete ray1;
}
 
if(show_rand) {
474,7 → 438,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;
483,10 → 447,10
phiMax *= pi/180.0;
 
TDatime now;
TRandom2 rand;
TRandom rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double phi;
double rphi;
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
514,19 → 478,19
ry = rand.Uniform(-b/2.0, b/2.0);
rz = rand.Uniform(-b/2.0, b/2.0);
 
phi = rand.Uniform(phiMin, phiMax);
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(phi);
rn = TMath::Sin(theta)*TMath::Cos(phi);
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(phi);
hphi->Fill(rphi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
534,12 → 498,11
 
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);
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateX(rphi);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
548,8 → 511,8
detector->Propagate(*ray0, ray1, 0);
}
 
//delete ray1;
delete ray0;
delete ray1;
}
 
if(show_rand) {