Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 69 → Rev 70

/lightguide/trunk/src/raySimulator.cpp
87,8 → 87,8
{
if(show_data) {
char sbuff[256];
sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
parameters.getA(), parameters.getM(), parameters.getD(),
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) {
183,6 → 183,7
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);
195,6 → 196,9
detector->Propagate(*ray0, ray1, show_3d);
delete ray0;
delete ray1;
return (detector->GetHActive())->GetEntries() / (double)(1);
}
//-----------------------------------------------------------------------------
224,16 → 228,26
if(i == (NN/2))
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
detector->Propagate(*ray0, ray1, 0);
delete ray0;
delete ray1;
}
}
return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
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;
}
//-----------------------------------------------------------------------------
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
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;
272,13 → 286,25
//delete ray0;
//delete ray1;
}
return (detector->GetHActive())->GetEntries() / (double)NN;
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
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
// 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 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;
292,6 → 318,7
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);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
322,7 → 349,9
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)));
//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);
345,8 → 374,8
detector->Propagate(*ray0, ray1, 0);
}
//delete ray0;
//delete ray1;
delete ray0;
delete ray1;
}
if(show_rand) {
363,8 → 392,113
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;
double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
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);
 
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);
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;
}