Rev 70 | Rev 73 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 25 | f9daq | 1 | #include "include/guide.h" |
| 2 | #include "src/raySimulator.cpp" |
||
| 3 | |||
| 4 | #include "TFile.h" |
||
| 5 | #include "TPolyMarker3D.h" |
||
| 6 | #include "TCanvas.h" |
||
| 7 | |||
| 8 | //extern int show_3d; |
||
| 9 | //extern int show_data; |
||
| 10 | |||
| 11 | // Set the simulation parameters |
||
| 70 | f9daq | 12 | void showVisual(int b) {show_3d = b;} |
| 13 | void showData(int b) {show_data = b;} |
||
| 25 | f9daq | 14 | |
| 72 | f9daq | 15 | // Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap, badCoupling) |
| 16 | DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false); |
||
| 70 | f9daq | 17 | // Print the detector parameters |
| 18 | void getParameters(); |
||
| 25 | f9daq | 19 | |
| 20 | //void SetLGType(int in = 1, int side = 1, int out = 0) |
||
| 21 | //{detector->SetLGType(in, side, out);} |
||
| 22 | |||
| 70 | f9daq | 23 | void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46) |
| 24 | { parameters.setGuide(SiPM0, b0, d0); |
||
| 25 | parameters.setIndices(n1, n2, n3); } |
||
| 25 | f9daq | 26 | |
| 70 | f9daq | 27 | void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
| 25 | f9daq | 28 | { parameters.setGap(x_gap0, y_gap0, z_gap0); } |
| 29 | |||
| 70 | f9daq | 30 | void setGlass(double glassOn, double glassD) |
| 25 | f9daq | 31 | { parameters.setGlass(glassOn, glassD); }; |
| 32 | |||
| 70 | f9daq | 33 | void setPlate(int plateOn = 1, double plateWidth = 1) |
| 25 | f9daq | 34 | { parameters.setPlate(plateOn, plateWidth); } |
| 35 | |||
| 70 | f9daq | 36 | void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); } |
| 25 | f9daq | 37 | |
| 70 | f9daq | 38 | // Refractive Indices |
| 39 | // n1 - around light guide - air |
||
| 40 | // n2 - light guide (plate) material - k9 glass |
||
| 41 | // n3 - the material at the exit - optical grease, epoxy, air, etc. |
||
| 42 | void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);} |
||
| 43 | |||
| 25 | f9daq | 44 | int save_ary = 0; |
| 45 | //------------------------------------------------------------------------------------------ |
||
| 46 | void Help() |
||
| 47 | { |
||
| 48 | printf("void SetCenter(x, y, z)\n"); |
||
| 49 | printf("void SetLGType(in, side, out)\n"); |
||
| 50 | printf("SURF_DUMMY 0, "); |
||
| 51 | printf("SURF_REFRA 1, "); |
||
| 52 | printf("SURF_REFLE 2, "); |
||
| 53 | printf("SURF_TOTAL 3, "); |
||
| 54 | printf("SURF_IMPER 4\n"); |
||
| 55 | printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n"); |
||
| 56 | printf("void SetR(R0)\n"); |
||
| 57 | printf("void SetGlass(glass_on0, glass_d0)\n"); |
||
| 58 | printf("void SetGap(x_gap0, y_gap0, z_gap0)\n"); |
||
| 59 | |||
| 60 | printf("LGG(NN, theta)\n"); |
||
| 61 | printf("LGR(NN, theta, phi)\n"); |
||
| 62 | printf("LGI(NN, theta)\n"); |
||
| 63 | |||
| 64 | printf("LGI_gap(NN, min, max, steps, theta)\n"); |
||
| 65 | printf("LGR_gap(NN, min, max, steps, theta, phi)\n"); |
||
| 66 | printf("LGR_th(NN, min, max, steps, phi)\n"); |
||
| 67 | printf("LGG_th(NN, min, max, steps)\n"); |
||
| 68 | |||
| 69 | printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)"); |
||
| 70 | } |
||
| 71 | //------------------------------------------------------------------------------------------ |
||
| 72 | int show_in_steps = 0; |
||
| 73 | void SetShowInSteps(int in = 1) |
||
| 74 | { |
||
| 75 | show_in_steps = in; |
||
| 76 | } |
||
| 77 | //------------------------------------------------------------------------------------------ |
||
| 78 | int gs_set = 0; |
||
| 79 | void SetGS() |
||
| 80 | { |
||
| 81 | const UInt_t Number = 2; |
||
| 82 | Double_t Red[Number] = { 1.00, 0.00}; |
||
| 83 | Double_t Green[Number] = { 1.00, 0.00}; |
||
| 84 | Double_t Blue[Number] = { 1.00, 0.00}; |
||
| 85 | Double_t Stops[Number] = { 0.00, 1.00}; |
||
| 86 | |||
| 87 | Int_t nb=50; |
||
| 88 | if(!gs_set) { |
||
| 89 | TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb); |
||
| 90 | gs_set = 1; |
||
| 91 | } |
||
| 92 | } |
||
| 93 | //------------------------------------------------------------------------------------------ |
||
| 94 | void Draw1(CDetector *detector, int gs = 0, double range = 0.0) |
||
| 95 | { |
||
| 96 | if(gs) SetGS(); |
||
| 97 | TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500); |
||
| 98 | gStyle->SetOptTitle(0); |
||
| 99 | gStyle->SetOptStat(0); |
||
| 100 | cdata->cd(0); |
||
| 101 | if(range > 0.1) { |
||
| 102 | ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range); |
||
| 103 | ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range); |
||
| 104 | } |
||
| 105 | (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]"); |
||
| 106 | (detector->GetHLaser())->Draw("COLZ"); |
||
| 107 | cdata->SaveAs("2dscan.gif"); |
||
| 108 | //gStyle->SetOptTitle(1); |
||
| 109 | //gStyle->SetOptStat(1); |
||
| 110 | } |
||
| 111 | |||
| 112 | //------------------------------------------------------------------------------------------ |
||
| 113 | |||
| 114 | |||
| 115 | TVector3 p_pol0(0.0, 0.0, 1.0); |
||
| 116 | void SetPol(double x, double y, double z) |
||
| 117 | {p_pol0.SetXYZ(x,y,z);} |
||
| 118 | |||
| 70 | f9daq | 119 | // Test function |
| 120 | // creates an optical boundary surface |
||
| 121 | // shows the propagation of light ray |
||
| 122 | // and polarization state of the incident ray (green) |
||
| 123 | // and surface normal (blue) |
||
| 25 | f9daq | 124 | void PolTest(double theta = 0.0) |
| 125 | { |
||
| 126 | int p_type = 1; |
||
| 70 | f9daq | 127 | double p_n1 = parameters.getN1(); |
| 128 | double p_n2 = parameters.getN2(); |
||
| 25 | f9daq | 129 | theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6; |
| 130 | TVector3 p_pol; |
||
| 131 | |||
| 132 | Init(); |
||
| 133 | |||
| 134 | double cx = 0.0; |
||
| 135 | TVector3 vodnik_edge[4]; |
||
| 136 | double t = 3.0; |
||
| 137 | vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t); |
||
| 138 | vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t); |
||
| 139 | /* |
||
| 140 | #define SURF_DUMMY 0 |
||
| 141 | #define SURF_REFRA 1 |
||
| 142 | #define SURF_REFLE 2 |
||
| 143 | #define SURF_TOTAL 3 |
||
| 144 | #define SURF_IMPER 4 |
||
| 145 | */ |
||
| 146 | CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN(); |
||
| 70 | f9daq | 147 | surf->SetFresnel(1); |
| 25 | f9daq | 148 | surf->Draw(); |
| 149 | |||
| 150 | CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta)); |
||
| 70 | f9daq | 151 | ray->SetColor(kBlack); |
| 25 | f9daq | 152 | //p_pol = rotatey(p_pol0, -theta); |
| 153 | p_pol = p_pol0; p_pol.RotateY(-theta); |
||
| 154 | //printf("p_pol = "); printv(p_pol); |
||
| 155 | ray->SetPolarization(p_pol); |
||
| 156 | |||
| 157 | ray->DrawS(cx, -5.0); |
||
| 158 | |||
| 70 | f9daq | 159 | CRay *out = new CRay; out->SetColor(kRed); |
| 25 | f9daq | 160 | TVector3 *inters = new TVector3; |
| 70 | f9daq | 161 | surf->PropagateRay(*ray, out, inters); |
| 162 | //if(fate == 1) out->DrawS(cx, 5.0); |
||
| 163 | out->DrawS(cx, 5.0); |
||
| 164 | |||
| 165 | CRay *incidentPolarization = new CRay; |
||
| 166 | incidentPolarization->SetColor(kGreen); |
||
| 167 | incidentPolarization->Set(ray->GetR(), p_pol); |
||
| 168 | incidentPolarization->DrawS(cx, 1.0); |
||
| 25 | f9daq | 169 | |
| 70 | f9daq | 170 | CRay *surfaceNormal = new CRay; |
| 171 | surfaceNormal->SetColor(kBlue); |
||
| 172 | surfaceNormal->Set(ray->GetR(), surf->GetN()); |
||
| 173 | surfaceNormal->DrawS(cx, 1.0); |
||
| 25 | f9daq | 174 | } |
| 175 | |||
| 176 | void ptt() |
||
| 177 | { |
||
| 178 | for(double th = 0.0; th < 91.0; th += 5.0) { |
||
| 179 | printf("%lf ", th); |
||
| 180 | PolTest(th); |
||
| 181 | } |
||
| 182 | } |
||
| 183 | //------------------------------------------------------------------------------------------ |
||
| 184 | // Propagate single ray and show the statistics |
||
| 185 | void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0) |
||
| 186 | { |
||
| 187 | CDetector *detector = new CDetector(CENTER, parameters); |
||
| 188 | Init(); |
||
| 189 | double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi); |
||
| 190 | PrintGuideHead(); |
||
| 54 | f9daq | 191 | PrintGuideStat(izkoristek); |
| 25 | f9daq | 192 | DrawData(detector, parameters, theta, izkoristek); |
| 193 | } |
||
| 194 | //------------------------------------------------------------------------------------------ |
||
| 195 | // Propagate NN rays generated as grid and show the statistics |
||
| 72 | f9daq | 196 | void LGG(int NN=10, double theta=0.0, bool coupling=false) |
| 25 | f9daq | 197 | { |
| 72 | f9daq | 198 | parameters.setCoupling(coupling); |
| 199 | CDetector *detector = new CDetector(CENTER, parameters); |
||
| 25 | f9daq | 200 | Init(); |
| 201 | double izkoristek = Grid(detector, parameters, NN, theta); |
||
| 202 | //printf("izkoristek = %.3lf\n", izkoristek); |
||
| 203 | PrintGuideHead(); |
||
| 54 | f9daq | 204 | PrintGuideStat(izkoristek); |
| 25 | f9daq | 205 | DrawData(detector, parameters, theta, izkoristek); |
| 206 | } |
||
| 207 | //------------------------------------------------------------------------------------------ |
||
| 208 | // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics |
||
| 209 | void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0) |
||
| 210 | { |
||
| 72 | f9daq | 211 | CDetector *detector = new CDetector(CENTER, parameters); |
| 25 | f9daq | 212 | Init(); |
| 70 | f9daq | 213 | double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30); |
| 25 | f9daq | 214 | //printf("izkoristek = %.3lf\n", izkoristek); |
| 54 | f9daq | 215 | PrintGuideHead(); PrintGuideStat(izkoristek); |
| 25 | f9daq | 216 | DrawData(detector, parameters, theta, izkoristek); |
| 217 | } |
||
| 218 | //------------------------------------------------------------------------------------------ |
||
| 219 | // Propagate NN rays isotropically generated in solid angle theta and show the statistics |
||
| 220 | void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0) |
||
| 221 | { |
||
| 222 | Init(); |
||
| 223 | CDetector *detector = new CDetector(CENTER, parameters); |
||
| 224 | //CDetector detector = new CDetector(); |
||
| 70 | f9daq | 225 | double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr); |
| 25 | f9daq | 226 | //printf("izkoristek = %.3lf\n", izkoristek); |
| 70 | f9daq | 227 | PrintGuideHead(); |
| 228 | PrintGuideStat(izkoristek); |
||
| 25 | f9daq | 229 | DrawData(detector, parameters, theta, izkoristek); |
| 230 | //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
||
| 231 | //canvasDetector->cd(); |
||
| 232 | //detector->Draw(); |
||
| 233 | } |
||
| 234 | |||
| 72 | f9daq | 235 | void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0) |
| 70 | f9daq | 236 | { |
| 237 | Init(); |
||
| 72 | f9daq | 238 | parameters.setCoupling(coupling); |
| 70 | f9daq | 239 | CDetector *detector = new CDetector(CENTER, parameters); |
| 72 | f9daq | 240 | |
| 241 | const double theta = 18.5; |
||
| 242 | double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr); |
||
| 243 | |||
| 244 | PrintGuideHead(); |
||
| 245 | PrintGuideStat(izkoristek); |
||
| 246 | DrawData(detector, parameters, 18.5, izkoristek); |
||
| 247 | |||
| 248 | //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
||
| 249 | //canvasDetector->cd(); |
||
| 250 | //detector->Draw(); |
||
| 70 | f9daq | 251 | } |
| 25 | f9daq | 252 | |
| 253 | |||
| 254 | //------------------------------------------------------------------------------------------ |
||
| 255 | void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0) |
||
| 256 | { |
||
| 257 | int tmp3d = show_3d, tmpdata = show_data; |
||
| 258 | double step[steps], acc[steps]; |
||
| 259 | |||
| 260 | show_3d = 0; show_data = 0; |
||
| 261 | |||
| 262 | //printf(" x_gap | acceptance \n"); |
||
| 263 | PrintGuideHead(); |
||
| 264 | for(int i=0; i<=steps; i++) { |
||
| 265 | step[i] = min + i*(max - min)/(double)steps; |
||
| 266 | parameters.setGap(step[i], 0.0, 0.0); |
||
| 267 | CDetector *detector = new CDetector(CENTER, parameters); |
||
| 268 | Init(); |
||
| 70 | f9daq | 269 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
| 25 | f9daq | 270 | //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0); |
| 54 | f9daq | 271 | PrintGuideStat(acc[i]); |
| 25 | f9daq | 272 | } |
| 273 | |||
| 274 | //char sbuff[256]; |
||
| 275 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
| 276 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 277 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
| 278 | |||
| 279 | DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max); |
||
| 280 | |||
| 281 | show_3d = tmp3d; show_data = tmpdata; |
||
| 282 | } |
||
| 283 | //------------------------------------------------------------------------------------------ |
||
| 284 | |||
| 70 | f9daq | 285 | void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
| 286 | { |
||
| 287 | //int tmp3d = show_3d, tmpdata = show_data; |
||
| 288 | double show_rays = 10; |
||
| 289 | double step[steps], acc[steps]; |
||
| 290 | |||
| 291 | //show_3d = 0; show_data = 0; |
||
| 292 | |||
| 293 | //printf(" theta | acceptance \n"); |
||
| 294 | PrintGuideHead(); |
||
| 295 | for(int i=0; i<=steps; i++) { |
||
| 296 | Init(); |
||
| 297 | step[i] = min + i*(max - min)/(double)steps; |
||
| 298 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
| 299 | acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays); |
||
| 300 | //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
||
| 301 | PrintGuideStat(acc[i]); |
||
| 302 | delete detector; |
||
| 303 | } |
||
| 304 | |||
| 305 | //char sbuff[256]; |
||
| 306 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)", |
||
| 307 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 308 | // (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z()); |
||
| 309 | |||
| 310 | DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
| 311 | |||
| 312 | //show_3d = tmp3d; show_data = tmpdata; |
||
| 313 | } |
||
| 314 | |||
| 315 | |||
| 25 | f9daq | 316 | //------------------------------------------------------------------------------------------ |
| 70 | f9daq | 317 | //void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
| 318 | void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15) |
||
| 25 | f9daq | 319 | { |
| 70 | f9daq | 320 | //int tmp3d = show_3d, tmpdata = show_data; |
| 321 | double step[steps], acc[steps]; |
||
| 25 | f9daq | 322 | |
| 70 | f9daq | 323 | //show_3d = 0; show_data = 0; |
| 25 | f9daq | 324 | |
| 325 | //printf(" theta | acceptance \n"); |
||
| 326 | PrintGuideHead(); |
||
| 327 | for(int i=0; i<=steps; i++) { |
||
| 328 | Init(); |
||
| 329 | step[i] = min + i*(max - min)/(double)steps; |
||
| 330 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
| 70 | f9daq | 331 | //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30); |
| 332 | acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d); |
||
| 25 | f9daq | 333 | //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
| 54 | f9daq | 334 | PrintGuideStat(acc[i]); |
| 70 | f9daq | 335 | delete detector; |
| 25 | f9daq | 336 | } |
| 337 | |||
| 338 | //char sbuff[256]; |
||
| 339 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)", |
||
| 340 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 341 | // (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z()); |
||
| 342 | |||
| 343 | DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
| 344 | |||
| 70 | f9daq | 345 | //show_3d = tmp3d; show_data = tmpdata; |
| 25 | f9daq | 346 | } |
| 347 | //------------------------------------------------------------------------------------------ |
||
| 348 | |||
| 70 | f9daq | 349 | void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10) |
| 350 | { |
||
| 351 | int tmp3d = show_3d, tmpdata = show_data; |
||
| 352 | show_3d = 1; show_data = 0; |
||
| 353 | //PrintGuideHead(); |
||
| 354 | TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi); |
||
| 355 | double min = 0.0; |
||
| 356 | printf("\nWait, this takes a while "); |
||
| 357 | for(int i=0; i<=steps; i++) { |
||
| 358 | double theta = min + i*(maxTheta - min) / (double)steps; |
||
| 359 | for (int j=0; j<=steps; j++) { |
||
| 360 | Init(); |
||
| 361 | double phi = min + j*(maxPhi - min) / (double)steps; |
||
| 362 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
| 363 | double acc = RandYZ(detector, parameters, NN, theta, phi, 30); |
||
| 364 | //PrintGuideStat(acc); |
||
| 365 | hAcceptance->Fill(i, j, acc); |
||
| 366 | //printf("Acc: %f ", acc); |
||
| 367 | delete detector; |
||
| 368 | } |
||
| 369 | printf("."); |
||
| 370 | } |
||
| 371 | printf("\n"); |
||
| 372 | //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
||
| 373 | show_3d = tmp3d; show_data = tmpdata; |
||
| 374 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000); |
||
| 375 | cp->cd(); |
||
| 376 | hAcceptance->Draw("COLZ"); |
||
| 377 | } |
||
| 378 | |||
| 25 | f9daq | 379 | // a vs. d |
| 380 | void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0) |
||
| 381 | { |
||
| 382 | //char sbuff[256]; |
||
| 383 | |||
| 70 | f9daq | 384 | //show_3d = 0; // don't show simulations |
| 385 | //show_data = 1; |
||
| 25 | f9daq | 386 | |
| 387 | //double d = detector->GetD(); |
||
| 388 | const double b = parameters.getB(); // upper side of LG |
||
| 389 | //const double SiPM = 3.0; // the length of the detector itself |
||
| 390 | //double reflectivity = detector->guide->getR(); |
||
| 391 | |||
| 392 | TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max); |
||
| 393 | |||
| 394 | // Use the Fresnel eq. instead of fixed reflectivity 96% |
||
| 395 | //detector->SetFresnel(1); |
||
| 396 | |||
| 70 | f9daq | 397 | //printf(" d | a | Acceptance\n"); |
| 398 | getParameters(); |
||
| 399 | printf("Wait, this takes a while "); |
||
| 25 | f9daq | 400 | for(int i=0; i<steps; i++) { |
| 401 | const double d = hAcceptance->GetXaxis()->GetBinCenter(i); |
||
| 402 | for(int j=0; j<steps; j++) { |
||
| 70 | f9daq | 403 | Init(); |
| 25 | f9daq | 404 | const double a = hAcceptance->GetYaxis()->GetBinCenter(j); |
| 70 | f9daq | 405 | //const double M = b/a; |
| 406 | parameters.setGuide(a, b, d); |
||
| 25 | f9daq | 407 | CDetector *detector = new CDetector(CENTER, parameters); |
| 408 | //detector->guide->setLG(y, M, x); |
||
| 409 | //Init(); exclude simulation |
||
| 70 | f9daq | 410 | double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d); |
| 25 | f9daq | 411 | //double acceptance = Grid(NN, theta); |
| 412 | //double acceptance = -1.0; |
||
| 413 | hAcceptance->Fill(d, a, acceptance); |
||
| 70 | f9daq | 414 | //printf("%.2lf | %.2lf | ", d, a); |
| 415 | //PrintGuideStat(acceptance); |
||
| 54 | f9daq | 416 | delete detector; //works fine, 50x50 grid takes ~4MB of RAM |
| 25 | f9daq | 417 | } |
| 70 | f9daq | 418 | printf("."); |
| 25 | f9daq | 419 | } |
| 70 | f9daq | 420 | printf("\n"); |
| 25 | f9daq | 421 | |
| 422 | |||
| 423 | |||
| 70 | f9daq | 424 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200); |
| 425 | cp->cd(); |
||
| 426 | //TVirtualPad *pacc = cp->cd(0); |
||
| 25 | f9daq | 427 | |
| 70 | f9daq | 428 | /* |
| 25 | f9daq | 429 | pacc->SetRightMargin(0.10); |
| 430 | pacc->SetLeftMargin(0.10); |
||
| 431 | pacc->SetTopMargin(0.10); |
||
| 70 | f9daq | 432 | pacc->SetBottomMargin(0.10); |
| 433 | */ |
||
| 25 | f9daq | 434 | |
| 435 | TFile *file = new TFile("acceptance.root","RECREATE"); |
||
| 436 | hAcceptance->Write(); |
||
| 437 | file->Close(); |
||
| 54 | f9daq | 438 | //delete file; |
| 25 | f9daq | 439 | |
| 440 | //hAcceptance->SetContour(100); |
||
| 441 | //gStyle->SetPalette(1,0); |
||
| 442 | hAcceptance->SetTitle(";d [mm];a [mm]"); |
||
| 443 | hAcceptance->GetXaxis()->SetRangeUser(minD,maxD); |
||
| 444 | hAcceptance->Draw("COLZ"); |
||
| 70 | f9daq | 445 | char filename[128]; |
| 446 | sprintf(filename,"LGI_ad%d.C", steps); |
||
| 447 | cp->SaveAs(filename); |
||
| 25 | f9daq | 448 | |
| 449 | } |
||
| 450 | |||
| 451 | //------------------------------------------------------------------------------------------ |
||
| 452 | // collection efficiency vs length |
||
| 453 | // a changes, b rests the same, d changes |
||
| 454 | // still can't use this function, it has hardcoded tan 10 deg |
||
| 455 | // !very stupid! |
||
| 456 | // new 3D graph NEEDED |
||
| 457 | // d vs a vs acceptance |
||
| 458 | void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0) |
||
| 459 | { |
||
| 460 | TVector3 gapGuideSiPM(0.3, 0, 0); |
||
| 461 | CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
||
| 462 | int tmp3d = show_3d, tmpdata = show_data; |
||
| 463 | double step[steps], acc[steps]; |
||
| 464 | double sipm_d, M0; |
||
| 465 | char sbuff[256]; |
||
| 466 | |||
| 467 | show_3d = 1; show_data = 1; |
||
| 468 | |||
| 469 | M0 = parameters.getM(); |
||
| 470 | |||
| 471 | PrintGuideHead(); |
||
| 472 | for(int i=0; i<=steps; i++) { |
||
| 473 | step[i] = min + i*(max - min)/(double)steps; |
||
| 474 | sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg |
||
| 70 | f9daq | 475 | parameters.setGuide(sipm_d, M0*sipm_d, step[i]); |
| 25 | f9daq | 476 | //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
| 477 | printf("%.1lf | %.2lf | ", step[i], sipm_d); |
||
| 478 | Init(); |
||
| 70 | f9daq | 479 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
| 54 | f9daq | 480 | PrintGuideStat(acc[i]); |
| 25 | f9daq | 481 | sprintf(sbuff, "d%2d.gif", i); |
| 482 | //c3dview->SaveAs(sbuff); |
||
| 483 | } |
||
| 484 | |||
| 485 | |||
| 486 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
| 487 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 488 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
| 489 | |||
| 490 | //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max); |
||
| 491 | |||
| 492 | TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480); |
||
| 493 | TVirtualPad *pacc = cp->cd(0); |
||
| 494 | |||
| 495 | pacc->SetRightMargin(0.05); |
||
| 496 | pacc->SetLeftMargin(0.13); |
||
| 497 | pacc->SetTopMargin(0.08); |
||
| 498 | |||
| 499 | TGraph *gacceptance = new TGraph(steps+1, step, acc); |
||
| 500 | gacceptance->SetTitle("; d [mm];Collection efficiency"); |
||
| 501 | gacceptance->SetMarkerStyle(8); |
||
| 502 | gacceptance->SetLineColor(12); |
||
| 503 | gacceptance->SetLineWidth(1); |
||
| 504 | (gacceptance->GetXaxis())->SetRangeUser(min, max); |
||
| 505 | //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
||
| 506 | |||
| 507 | gacceptance->Draw("ACP"); |
||
| 508 | |||
| 509 | |||
| 510 | show_3d = tmp3d; show_data = tmpdata; |
||
| 511 | } |
||
| 512 | //------------------------------------------------------------------------------------------ |
||
| 513 | // Acceptance vs. the length |
||
| 514 | // The magnification ratio M rests the same all the time |
||
| 515 | void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0) |
||
| 516 | { |
||
| 517 | int tmp3d = show_3d, tmpdata = show_data; |
||
| 518 | double step[steps], acc[steps]; |
||
| 519 | const double a = parameters.getA(); |
||
| 520 | const double magnif = parameters.getM(); |
||
| 521 | |||
| 522 | //show_3d = show_in_steps; show_data = 1; |
||
| 523 | |||
| 524 | // Use Fresnel equations |
||
| 525 | //detector->SetFresnel(1); |
||
| 526 | |||
| 527 | // Set glass (n=1.5) at the exit |
||
| 528 | //detector->SetGlass(1,0); |
||
| 529 | |||
| 530 | PrintGuideHead(); |
||
| 531 | for(int i=0; i<=steps; i++) { |
||
| 532 | step[i] = min + i*(max - min)/(double)steps; |
||
| 70 | f9daq | 533 | //parameters.setGuide(a, magnif, step[i]); |
| 25 | f9daq | 534 | parameters.setGuide(a, magnif, step[i]); |
| 535 | CDetector *detector = new CDetector(CENTER, parameters); |
||
| 536 | //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
||
| 537 | printf("%.1lf | ", step[i]); |
||
| 538 | Init(); |
||
| 70 | f9daq | 539 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
| 54 | f9daq | 540 | PrintGuideStat(acc[i]); |
| 25 | f9daq | 541 | delete detector; |
| 542 | } |
||
| 543 | |||
| 544 | //char sbuff[256]; |
||
| 545 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
| 546 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 547 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
| 548 | |||
| 549 | DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max); |
||
| 550 | |||
| 551 | show_3d = tmp3d; show_data = tmpdata; |
||
| 552 | } |
||
| 553 | //------------------------------------------------------------------------------------------ |
||
| 554 | // Magnification optimization. a and d are fixed, M and b change in steps |
||
| 555 | void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0) |
||
| 556 | { |
||
| 557 | int tmp3d = show_3d, tmpdata = show_data; |
||
| 558 | double step[steps], acc[steps]; |
||
| 559 | |||
| 560 | const double a = parameters.getA(); |
||
| 561 | const double d = parameters.getD(); |
||
| 562 | |||
| 563 | show_3d = 0; show_data = 0; |
||
| 564 | |||
| 565 | PrintGuideHead(); |
||
| 566 | for(int i=0; i<=steps; i++) { |
||
| 567 | step[i] = min + i*(max - min)/(double)steps; |
||
| 70 | f9daq | 568 | //parameters.setGuide(a, step[i], d); |
| 569 | parameters.setGuide(a, a*step[i], d); |
||
| 25 | f9daq | 570 | CDetector *detector = new CDetector(CENTER, parameters); |
| 571 | Init(); |
||
| 70 | f9daq | 572 | acc[i] = RandIso(detector, parameters, NN, 0, theta); |
| 54 | f9daq | 573 | PrintGuideStat(acc[i]); |
| 25 | f9daq | 574 | } |
| 575 | |||
| 576 | //char sbuff[256]; |
||
| 577 | //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
||
| 578 | // detector->GetSiPM(), detector->GetM(), detector->GetD(), |
||
| 579 | // (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
||
| 580 | |||
| 581 | DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max); |
||
| 582 | |||
| 583 | show_3d = tmp3d; show_data = tmpdata; |
||
| 584 | } |
||
| 70 | f9daq | 585 | |
| 25 | f9daq | 586 | //------------------------------------------------------------------------------------------ |
| 587 | |||
| 70 | f9daq | 588 | void getParameters() |
| 589 | { |
||
| 590 | printf("LIGHT GUIDE\n" |
||
| 591 | " b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA()); |
||
| 592 | printf("MATERIAL REFRACITVE INDICES\n" |
||
| 593 | " n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3()); |
||
| 594 | printf("PLATE\n" |
||
| 595 | " ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel()); |
||
| 596 | return; |
||
| 597 | } |
||
| 598 |