Rev 54 | Rev 71 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 25 | f9daq | 1 | // own include |
| 2 | //#include "include/raySimulator.h" |
||
| 3 | #include "include/RTUtil.h" |
||
| 4 | #include "include/guide.h" |
||
| 5 | |||
| 6 | // general include |
||
| 7 | #include "TROOT.h" |
||
| 8 | #include "TSystem.h" |
||
| 9 | #include "TStyle.h" |
||
| 10 | #include "TCanvas.h" |
||
| 11 | #include "TMath.h" |
||
| 12 | #include "TView3D.h" |
||
| 13 | #include "TH1F.h" |
||
| 14 | #include "TH2F.h" |
||
| 15 | #include "TBenchmark.h" |
||
| 16 | #include "TPolyMarker.h" |
||
| 17 | #include "TGraph.h" |
||
| 18 | #include "TF1.h" |
||
| 19 | |||
| 20 | //================================================================================= |
||
| 21 | TCanvas *c3dview; |
||
| 22 | TCanvas *cacc; |
||
| 23 | int show_3d = 0, show_data = 0; |
||
| 24 | int draw_width = 2; |
||
| 25 | |||
| 26 | // calls default constructor CVodnik.cpp |
||
| 27 | |||
| 28 | //void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0) |
||
| 29 | //{center.SetXYZ(x,y,z); detector = new CDetector(center);} |
||
| 30 | |||
| 31 | |||
| 32 | //void SetLGType(int in = 1, int side = 1, int out = 0) |
||
| 33 | //{detector->SetLGType(in, side, out);} |
||
| 34 | |||
| 35 | //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) |
||
| 36 | //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);} |
||
| 37 | |||
| 38 | //void SetR(double R0) {detector->SetR(R0);} |
||
| 39 | /* |
||
| 40 | void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3) |
||
| 41 | {detector->SetGlass(glass_on0, glass_d0);} |
||
| 42 | |||
| 43 | void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
||
| 44 | {detector->SetGap(x_gap0 , y_gap0, z_gap0);} |
||
| 45 | |||
| 46 | void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6) |
||
| 47 | {detector->SetRCol(in, lg, out, gla);}; |
||
| 48 | void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3) |
||
| 49 | {detector->SetDCol(LG0, glass0, active0);}; |
||
| 50 | |||
| 51 | |||
| 52 | void SetDWidth(int w = 1) {draw_width = w;} |
||
| 53 | |||
| 54 | void SetFresnel(int b = 1) {detector->SetFresnel(b);} |
||
| 55 | void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);} |
||
| 56 | |||
| 57 | void SetGuideOn(int b = 0) {detector->SetGuideOn(b);} |
||
| 58 | |||
| 59 | void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);} |
||
| 60 | void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);} |
||
| 61 | */ |
||
| 62 | //----------------------------------------------------------------------------- |
||
| 63 | void Init(double rsys_scale = 7.0) |
||
| 64 | { |
||
| 65 | RTSetStyle(gStyle); |
||
| 66 | gStyle->SetOptStat(10); |
||
| 67 | |||
| 68 | if(show_3d) { |
||
| 69 | double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale}; |
||
| 70 | double rsys1[]={ rsys_scale, rsys_scale, rsys_scale}; |
||
| 71 | |||
| 72 | c3dview = (TCanvas*)gROOT->FindObject("c3dview"); |
||
| 73 | if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723); c3dview->SetFillColor(0);} |
||
| 74 | else c3dview->Clear(); |
||
| 75 | |||
| 76 | TView3D *view = new TView3D(1, rsys0, rsys1); |
||
| 77 | view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]); |
||
| 78 | //view->SetPerspective(); |
||
| 79 | view->SetParallel(); |
||
| 80 | view->FrontView(); |
||
| 81 | view->Zoom();view->Zoom();view->Zoom();//view->Zoom(); |
||
| 82 | //view->ToggleRulers(); |
||
| 83 | } |
||
| 84 | } |
||
| 85 | //----------------------------------------------------------------------------- |
||
| 86 | void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0) |
||
| 87 | { |
||
| 88 | if(show_data) { |
||
| 89 | char sbuff[256]; |
||
| 70 | f9daq | 90 | sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf", |
| 91 | parameters.getA(), parameters.getLightYield()*acc, parameters.getD(), |
||
| 25 | f9daq | 92 | parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc); |
| 93 | |||
| 94 | if(!only2d) { |
||
| 95 | RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
||
| 96 | cdata->Divide(3,3); |
||
| 97 | int cpc = 1; |
||
| 98 | cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
||
| 99 | cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
||
| 100 | cdata->cd(cpc++); ((detector->GetLG())->GetHOut())->Draw("COLZ"); |
||
| 101 | |||
| 102 | |||
| 103 | cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
||
| 104 | cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
||
| 105 | cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
||
| 106 | |||
| 107 | cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
||
| 108 | cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
||
| 109 | cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
||
| 110 | }else{ |
||
| 111 | RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500); |
||
| 112 | cdata->cd(0); (detector->GetHLaser())->Draw("COLZ"); |
||
| 113 | cdata->SaveAs("cdata.gif"); |
||
| 114 | } |
||
| 115 | } |
||
| 116 | } |
||
| 117 | //----------------------------------------------------------------------------- |
||
| 118 | void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance", |
||
| 119 | double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0) |
||
| 120 | { |
||
| 121 | |||
| 122 | cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480); |
||
| 123 | |||
| 124 | (cacc->cd(0))->SetRightMargin(0.05); |
||
| 125 | (cacc->cd(0))->SetLeftMargin(0.13); |
||
| 126 | (cacc->cd(0))->SetTopMargin(0.08); |
||
| 127 | |||
| 128 | TGraph *gacceptance = new TGraph(n, x, y); |
||
| 129 | gacceptance->SetTitle(htitle); |
||
| 130 | gacceptance->SetMarkerStyle(8); |
||
| 131 | gacceptance->SetLineColor(12); |
||
| 132 | gacceptance->SetLineWidth(1); |
||
| 133 | if(xmin < xmax) |
||
| 134 | (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax); |
||
| 135 | (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
||
| 136 | |||
| 137 | gacceptance->Draw("ACP"); |
||
| 138 | //cacc->Print("acceptance.pdf"); |
||
| 139 | //delete gacceptance; |
||
| 140 | //delete cacc; |
||
| 141 | } |
||
| 142 | //----------------------------------------------------------------------------- |
||
| 143 | void PrintGlassStat(CDetector *detector) |
||
| 144 | { |
||
| 145 | int glass_in = ((detector->GetLG())->GetHOut())->GetEntries(); |
||
| 146 | int glass_out = (detector->GetHGlass())->GetEntries(); |
||
| 147 | printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in); |
||
| 148 | } |
||
| 149 | //----------------------------------------------------------------------------- |
||
| 150 | void PrintGuideHead() |
||
| 151 | { |
||
| 152 | //printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
||
| 153 | printf("Acceptance\n"); |
||
| 154 | } |
||
| 155 | //----------------------------------------------------------------------------- |
||
| 54 | f9daq | 156 | void PrintGuideStat(double acceptance) |
| 25 | f9daq | 157 | { |
| 158 | /*int n_active = (detector->GetHActive())->GetEntries(); |
||
| 159 | int n_enter = (detector->GetLG())->GetEnteranceHits(); |
||
| 160 | double izkoristek = n_active/(double)n_enter; |
||
| 161 | int fates[7]; (detector->GetLG())->GetVFate(fates); |
||
| 162 | |||
| 163 | // printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
||
| 164 | |||
| 165 | printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |", |
||
| 166 | (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi); |
||
| 167 | for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter); |
||
| 168 | |||
| 169 | printf("%5.1lf\n", 100.0*izkoristek);*/ |
||
| 170 | |||
| 54 | f9daq | 171 | //double n_active = (detector->GetHActive())->GetEntries(); |
| 172 | //double r_acc = 100.0 * n_active / n_in; |
||
| 173 | double r_acc = 100.0 * acceptance; |
||
| 25 | f9daq | 174 | printf("%7.3lf\n", r_acc); |
| 175 | } |
||
| 176 | //----------------------------------------------------------------------------- |
||
| 177 | |||
| 178 | //----------------------------------------------------------------------------- |
||
| 179 | // en zarek |
||
| 180 | double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0) |
||
| 181 | { |
||
| 182 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
| 183 | theta = Pi()*theta/180.0; |
||
| 184 | if(theta < 1e-6) theta = 1e-6; |
||
| 185 | phi = phi*Pi()/180.0; |
||
| 70 | f9daq | 186 | if(phi < 1e-6) phi = 1e-6; |
| 25 | f9daq | 187 | |
| 188 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
| 189 | |||
| 190 | //detector->Init(); |
||
| 191 | if(show_3d) detector->Draw(draw_width); |
||
| 192 | |||
| 193 | CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(), |
||
| 194 | TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi)); |
||
| 195 | CRay *ray1 = new CRay; |
||
| 196 | |||
| 197 | detector->Propagate(*ray0, ray1, show_3d); |
||
| 198 | |||
| 70 | f9daq | 199 | delete ray0; |
| 200 | delete ray1; |
||
| 201 | |||
| 25 | f9daq | 202 | return (detector->GetHActive())->GetEntries() / (double)(1); |
| 203 | } |
||
| 204 | //----------------------------------------------------------------------------- |
||
| 205 | // zarki, razporejeni v mrezi |
||
| 206 | double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0) |
||
| 207 | { |
||
| 208 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
| 209 | theta = Pi()*theta/180.0; |
||
| 210 | if(theta < 1e-6) theta = 1e-6; |
||
| 211 | |||
| 212 | //detector->Init(); |
||
| 213 | if(show_3d) detector->Draw(draw_width); |
||
| 214 | |||
| 215 | const double b = parameters.getB(); |
||
| 216 | |||
| 217 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
| 218 | |||
| 219 | for(int i = 0; i < NN; i++) { |
||
| 220 | for(int j = 0; j < NN; j++) { |
||
| 221 | CRay *ray0 = new CRay(CENTER.x() - offset, |
||
| 222 | 0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
||
| 223 | 0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
||
| 224 | TMath::Cos(theta), |
||
| 225 | 0.0, |
||
| 226 | TMath::Sin(theta)); |
||
| 227 | CRay *ray1 = new CRay; |
||
| 228 | if(i == (NN/2)) |
||
| 229 | detector->Propagate(*ray0, ray1, show_3d); |
||
| 230 | else |
||
| 70 | f9daq | 231 | detector->Propagate(*ray0, ray1, 0); |
| 232 | delete ray0; |
||
| 233 | delete ray1; |
||
| 25 | f9daq | 234 | } |
| 70 | f9daq | 235 | |
| 25 | f9daq | 236 | } |
| 70 | f9daq | 237 | double acceptance = 0.0; |
| 238 | /* |
||
| 239 | if( !(parameters.getPlateOn()) ) |
||
| 240 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN; |
||
| 241 | else |
||
| 242 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
| 243 | */ |
||
| 244 | acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
||
| 245 | return acceptance; |
||
| 25 | f9daq | 246 | } |
| 247 | //----------------------------------------------------------------------------- |
||
| 248 | // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika) |
||
| 249 | // vsi pod kotom (theta, phi) |
||
| 70 | f9daq | 250 | double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays) |
| 25 | f9daq | 251 | { |
| 252 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
| 253 | theta = theta*3.14159265358979312/180.0; |
||
| 254 | if(theta < MARGIN) theta = MARGIN; |
||
| 255 | phi = phi*3.14159265358979312/180.0; |
||
| 256 | if(phi < MARGIN) phi = MARGIN; |
||
| 257 | |||
| 258 | TDatime now; |
||
| 259 | TRandom rand; |
||
| 260 | rand.SetSeed(now.Get()); |
||
| 261 | |||
| 262 | //detector->Init(); |
||
| 263 | if(show_3d) detector->Draw(draw_width); |
||
| 264 | |||
| 265 | double SiPM = parameters.getA(); |
||
| 266 | double M = parameters.getM(); |
||
| 267 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
| 268 | |||
| 269 | for(int i = 0; i < NN; i++) { |
||
| 270 | |||
| 271 | double rx = CENTER.x() - offset; |
||
| 272 | double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
| 273 | double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
| 274 | |||
| 275 | double rl = TMath::Cos(theta); |
||
| 276 | double rm = TMath::Sin(theta)*TMath::Sin(phi); |
||
| 277 | double rn = TMath::Sin(theta)*TMath::Cos(phi); |
||
| 278 | |||
| 279 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
| 280 | CRay *ray1 = new CRay; |
||
| 281 | |||
| 282 | if(i < show_rays) |
||
| 283 | detector->Propagate(*ray0, ray1, show_3d); |
||
| 284 | else |
||
| 285 | detector->Propagate(*ray0, ray1, 0); |
||
| 54 | f9daq | 286 | //delete ray0; |
| 287 | //delete ray1; |
||
| 25 | f9daq | 288 | } |
| 70 | f9daq | 289 | |
| 290 | double acceptance = 0.0; |
||
| 291 | /* |
||
| 292 | if( !(parameters.getPlateOn()) ) |
||
| 293 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
| 294 | else |
||
| 295 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
| 296 | //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/ |
||
| 297 | acceptance = (detector->GetHActive())->GetEntries() / (double) NN; |
||
| 298 | return acceptance; |
||
| 25 | f9daq | 299 | } |
| 300 | //----------------------------------------------------------------------------- |
||
| 301 | // zarki, izotropno porazdeljeni znotraj kota theta |
||
| 302 | // = nakljucni vstopni polozaj in kot |
||
| 70 | f9daq | 303 | // NN - number of rays to be simulated with angles theta distributed uniformly |
| 304 | // in cos theta and phi uniformly: |
||
| 305 | // theta [0, theta] |
||
| 306 | // phi [0,360] |
||
| 307 | 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) |
||
| 25 | f9daq | 308 | { |
| 309 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
| 310 | double pi = 3.14159265358979312; |
||
| 311 | theta = theta*3.14159265358979312/180.0; |
||
| 312 | if(theta < MARGIN) theta = MARGIN; |
||
| 313 | |||
| 314 | TDatime now; |
||
| 315 | TRandom rand; |
||
| 316 | rand.SetSeed(now.Get()); |
||
| 317 | double rx, ry, rz, rl, rm, rn; |
||
| 318 | double rphi, rtheta; |
||
| 319 | //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
||
| 320 | double theta_min_rad = TMath::Cos(theta); |
||
| 70 | f9daq | 321 | double theta_max_rad = TMath::Cos(thetaMin); |
| 25 | f9daq | 322 | |
| 323 | TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
||
| 324 | hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
||
| 325 | hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi); |
||
| 326 | htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
||
| 327 | htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
||
| 328 | hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
||
| 329 | hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
||
| 330 | hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
||
| 331 | hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
||
| 332 | hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
||
| 333 | hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
||
| 334 | hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
||
| 335 | hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
||
| 336 | |||
| 337 | |||
| 338 | //detector->Init(); |
||
| 339 | if (show_3d) detector->Draw(draw_width); |
||
| 340 | |||
| 341 | double SiPM = parameters.getA(); |
||
| 342 | double M = parameters.getM(); |
||
| 343 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
| 344 | |||
| 345 | for(int i = 0; i < NN; i++) { |
||
| 346 | |||
| 347 | rx = CENTER.x() - offset; |
||
| 348 | ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
| 349 | rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
| 350 | |||
| 351 | rphi = rand.Uniform(0.0, 2.0*pi); |
||
| 70 | f9daq | 352 | //double theta_max_rad = TMath::Cos(17.0*pi/180.0); |
| 353 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
||
| 354 | rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad))); |
||
| 25 | f9daq | 355 | |
| 356 | rl = TMath::Cos(rtheta); |
||
| 357 | rm = TMath::Sin(rtheta)*TMath::Sin(rphi); |
||
| 358 | rn = TMath::Sin(rtheta)*TMath::Cos(rphi); |
||
| 359 | |||
| 360 | if(show_rand) { |
||
| 361 | htheta->Fill(rtheta); |
||
| 362 | hcostheta->Fill( TMath::Cos(rtheta) ); |
||
| 363 | hphi->Fill(rphi); |
||
| 364 | hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
||
| 365 | } |
||
| 366 | |||
| 367 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
| 368 | CRay *ray1 = new CRay; |
||
| 369 | |||
| 370 | if(i < show_rays) { |
||
| 371 | detector->Propagate(*ray0, ray1, show_3d); |
||
| 372 | } |
||
| 373 | else { |
||
| 374 | detector->Propagate(*ray0, ray1, 0); |
||
| 375 | } |
||
| 376 | |||
| 70 | f9daq | 377 | delete ray0; |
| 378 | delete ray1; |
||
| 25 | f9daq | 379 | } |
| 380 | |||
| 381 | if(show_rand) { |
||
| 382 | TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
||
| 383 | if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
||
| 384 | else c2rand->Clear(); |
||
| 385 | c2rand->Divide(2,3); |
||
| 386 | |||
| 387 | c2rand->cd(1); hphi->Draw(); |
||
| 388 | c2rand->cd(3); htheta->Draw(); |
||
| 389 | c2rand->cd(5); hcostheta->Draw(); |
||
| 390 | c2rand->cd(2); hl->Draw(); |
||
| 391 | c2rand->cd(4); hm->Draw(); |
||
| 392 | c2rand->cd(6); hn->Draw(); |
||
| 393 | } |
||
| 394 | |||
| 70 | f9daq | 395 | double acceptance = 0.0; |
| 396 | /* |
||
| 397 | if( !(parameters.getPlateOn()) ) |
||
| 398 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
| 399 | else |
||
| 400 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
| 401 | */ |
||
| 402 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
| 403 | return acceptance; |
||
| 404 | } |
||
| 405 | |||
| 406 | // Beamtest distribution |
||
| 407 | // with fixed theta and phi in interval phiMin, phiMax |
||
| 408 | double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand) |
||
| 409 | { |
||
| 410 | double pi = 3.14159265358979312; |
||
| 411 | theta *= pi/180.0; |
||
| 412 | if(theta < MARGIN) theta = MARGIN; |
||
| 413 | phiMin *= pi/180.0; |
||
| 414 | phiMax *= pi/180.0; |
||
| 25 | f9daq | 415 | |
| 70 | f9daq | 416 | TDatime now; |
| 417 | TRandom rand; |
||
| 418 | rand.SetSeed(now.Get()); |
||
| 419 | double rx, ry, rz, rl, rm, rn; |
||
| 420 | double rphi; |
||
| 421 | //double rtheta; |
||
| 422 | //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
||
| 423 | //double theta_min_rad = TMath::Cos(theta); |
||
| 424 | |||
| 425 | TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
||
| 426 | hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
||
| 427 | hphi = new TH1F("hphi", "hphi", 100, -pi, pi); |
||
| 428 | htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
||
| 429 | htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
||
| 430 | hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
||
| 431 | hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
||
| 432 | hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
||
| 433 | hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
||
| 434 | hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
||
| 435 | hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
||
| 436 | hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
||
| 437 | hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
||
| 438 | |||
| 439 | //detector->Init(); |
||
| 440 | if (show_3d) detector->Draw(draw_width); |
||
| 441 | |||
| 442 | double b = parameters.getB(); |
||
| 443 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
| 444 | |||
| 445 | for(int i = 0; i < NN; i++) { |
||
| 446 | |||
| 447 | rx = CENTER.x() - offset; |
||
| 448 | ry = rand.Uniform(-b/2.0, b/2.0); |
||
| 449 | rz = rand.Uniform(-b/2.0, b/2.0); |
||
| 450 | |||
| 451 | rphi = rand.Uniform(phiMin, phiMax); |
||
| 452 | //double theta_max_rad = TMath::Cos(17.0*pi/180.0); |
||
| 453 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
||
| 454 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad))); |
||
| 455 | |||
| 456 | rl = TMath::Cos(theta); |
||
| 457 | rm = TMath::Sin(theta)*TMath::Sin(rphi); |
||
| 458 | rn = TMath::Sin(theta)*TMath::Cos(rphi); |
||
| 459 | |||
| 460 | if(show_rand) { |
||
| 461 | //htheta->Fill(rtheta); |
||
| 462 | //hcostheta->Fill( TMath::Cos(rtheta) ); |
||
| 463 | hphi->Fill(rphi); |
||
| 464 | hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
||
| 465 | } |
||
| 466 | |||
| 467 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
| 468 | CRay *ray1 = new CRay; |
||
| 469 | |||
| 470 | if(i < show_rays) { |
||
| 471 | detector->Propagate(*ray0, ray1, show_3d); |
||
| 472 | } |
||
| 473 | else { |
||
| 474 | detector->Propagate(*ray0, ray1, 0); |
||
| 475 | } |
||
| 476 | |||
| 477 | delete ray0; |
||
| 478 | delete ray1; |
||
| 479 | } |
||
| 480 | |||
| 481 | if(show_rand) { |
||
| 482 | TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
||
| 483 | if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
||
| 484 | else c2rand->Clear(); |
||
| 485 | c2rand->Divide(2,3); |
||
| 486 | |||
| 487 | c2rand->cd(1); hphi->Draw(); |
||
| 488 | c2rand->cd(3); htheta->Draw(); |
||
| 489 | c2rand->cd(5); hcostheta->Draw(); |
||
| 490 | c2rand->cd(2); hl->Draw(); |
||
| 491 | c2rand->cd(4); hm->Draw(); |
||
| 492 | c2rand->cd(6); hn->Draw(); |
||
| 493 | } |
||
| 494 | |||
| 495 | double acceptance = 0.0; |
||
| 496 | /* |
||
| 497 | if( !(parameters.getPlateOn()) ) |
||
| 498 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
| 499 | else |
||
| 500 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
| 501 | */ |
||
| 502 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
| 25 | f9daq | 503 | return acceptance; |
| 504 | } |