WebSVN – f9daq – Diff – /lightguide/trunk/src/userFunctions.cpp

 //extern int show_3d;
 //extern int show_data;
 // Set the simulation parameters
-void Show3D(int b) {show_3d = b;}
+void showVisual(int b) {show_3d = b;}
-void ShowData(int b) {show_data = b;}
+void showData(int b) {show_data = b;}
 // Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap)
-DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.52, TVector3(0.3, 0, 0));
+DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0));
+// Print the detector parameters
+void getParameters();
 //void SetLGType(int in = 1, int side = 1, int out = 0)
         //{detector->SetLGType(in, side, out);}
-void SetLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n10 = 1.0, double n20 = 1.53, double n30 = 1.52)
+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)
-        { parameters.setGuide(SiPM0, b0, d0, n10, n20, n30); }
+        { parameters.setGuide(SiPM0, b0, d0);
+          parameters.setIndices(n1, n2, n3); }
-void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
+void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
   { parameters.setGap(x_gap0, y_gap0, z_gap0); }
-void SetGlass(double glassOn, double glassD)
+void setGlass(double glassOn, double glassD)
   { parameters.setGlass(glassOn, glassD); };
-void SetPlate(int plateOn = 1, double plateWidth = 1)
+void setPlate(int plateOn = 1, double plateWidth = 1)
   { parameters.setPlate(plateOn, plateWidth); }
-void SetFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
+void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
+// Refractive Indices
+// n1 - around light guide - air
+// n2 - light guide (plate) material - k9 glass
+// n3 - the material at the exit - optical grease, epoxy, air, etc.
+void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);}
 int save_ary = 0;
 //------------------------------------------------------------------------------------------
 void Help()
+{
 TVector3 p_pol0(0.0, 0.0, 1.0);
 void SetPol(double x, double y, double z)
 {p_pol0.SetXYZ(x,y,z);}
+// Test function
+// creates an optical boundary surface
+// shows the propagation of light ray
+// and polarization state of the incident ray (green)
+// and surface normal (blue)
 void PolTest(double theta = 0.0)
+{
   int p_type = 1;
-  double p_n1 = 1.5;
+  double p_n1 = parameters.getN1();
-  double p_n2 = 1.0;
+  double p_n2 = parameters.getN2();
         theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
         TVector3 p_pol;
         Init();
 #define SURF_REFLE 2
 #define SURF_TOTAL 3
 #define SURF_IMPER 4
         */
         CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
-        surf->SetFresnel(0);
+        surf->SetFresnel(1);
         surf->Draw();
         CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
-        ray->SetColor(1);
+        ray->SetColor(kBlack);
         //p_pol = rotatey(p_pol0, -theta);
         p_pol = p_pol0; p_pol.RotateY(-theta);
         //printf("p_pol = "); printv(p_pol);
         ray->SetPolarization(p_pol);
         ray->DrawS(cx, -5.0);
-        CRay *out = new CRay; out->SetColor(2);
+        CRay *out = new CRay; out->SetColor(kRed);
         TVector3 *inters = new TVector3;
-        int fate = surf->PropagateRay(*ray, out, inters);
+        surf->PropagateRay(*ray, out, inters);
-        if(fate == 1) out->DrawS(cx, 5.0);
+        //if(fate == 1) out->DrawS(cx, 5.0);
+  out->DrawS(cx, 5.0);
+        CRay *incidentPolarization = new CRay;
-        CRay *pp = new CRay; pp->SetColor(3);
+        incidentPolarization->SetColor(kGreen);
-        pp->Set(ray->GetR(), p_pol);
+        incidentPolarization->Set(ray->GetR(), p_pol);
-        pp->DrawS(cx, 2.0);
+        incidentPolarization->DrawS(cx, 1.0);
+        CRay *surfaceNormal = new CRay;
-        CRay *pn = new CRay; pn->SetColor(4);
+        surfaceNormal->SetColor(kBlue);
-        pn->Set(ray->GetR(), surf->GetN());
+        surfaceNormal->Set(ray->GetR(), surf->GetN());
-        pn->DrawS(cx, 3.0);
+        surfaceNormal->DrawS(cx, 1.0);
+}
 void ptt()
+{
         for(double th = 0.0; th < 91.0; th += 5.0) {
 // Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
 void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
+{
   CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
         Init();
-        double izkoristek = RandYZ(detector, parameters, NN, theta, phi);
+        double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
         //printf("izkoristek = %.3lf\n", izkoristek);
         PrintGuideHead(); PrintGuideStat(izkoristek);
         DrawData(detector, parameters, theta, izkoristek);
+}
 //------------------------------------------------------------------------------------------
 void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
+{
   Init();
   CDetector *detector = new CDetector(CENTER, parameters);
   //CDetector detector = new CDetector();
-        double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr);
+        double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
         //printf("izkoristek = %.3lf\n", izkoristek);
+        PrintGuideHead();
-        PrintGuideHead(); PrintGuideStat(izkoristek);
+        PrintGuideStat(izkoristek);
         DrawData(detector, parameters, theta, izkoristek);
         //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
         //canvasDetector->cd();
         //detector->Draw();
+}
+void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, int nnrays = 30, int showr = 0)
+{
+  Init();
+  CDetector *detector = new CDetector(CENTER, parameters);
+  //CDetector detector = new CDetector();
+        double izkoristek = beamtest(detector, parameters, NN, 18.5, phiMin, phiMax, nnrays, showr);
+        //printf("izkoristek = %.3lf\n", izkoristek);
+        PrintGuideHead();
+        PrintGuideStat(izkoristek);
+        DrawData(detector, parameters, 18.5, izkoristek);
+        //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
+        //canvasDetector->cd();
+        //detector->Draw();
+}
 //------------------------------------------------------------------------------------------
 void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
 {
 int tmp3d = show_3d, tmpdata = show_data;
 double step[steps], acc[steps];
         show_3d = 0; show_data = 0;
         //printf(" x_gap | acceptance \n");
         PrintGuideHead();
         for(int i=0; i<=steps; i++) {
                 step[i] = min + i*(max - min)/(double)steps;
                 parameters.setGap(step[i], 0.0, 0.0);
                 CDetector *detector = new CDetector(CENTER, parameters);
                 Init();
-                acc[i] = RandIso(detector, parameters, NN, theta);
+                acc[i] = RandIso(detector, parameters, NN, 0, theta);
                 //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
                 PrintGuideStat(acc[i]);
+        }
         //char sbuff[256];
         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
                 //      (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
         DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
         show_3d = tmp3d; show_data = tmpdata;
 }
 //------------------------------------------------------------------------------------------
+void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
+{
+  //int tmp3d = show_3d, tmpdata = show_data;
+  double show_rays = 10;
+  double step[steps], acc[steps];
+        //show_3d = 0; show_data = 0;
+        //printf(" theta | acceptance \n");
+        PrintGuideHead();
+        for(int i=0; i<=steps; i++) {
+                Init();
+                step[i] = min + i*(max - min)/(double)steps;
+                CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
+                acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
+                //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
+                PrintGuideStat(acc[i]);
+                delete detector;
+        }
+        //char sbuff[256];
+        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
+                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
+                //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
+        DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
+        //show_3d = tmp3d; show_data = tmpdata;
+}
 //------------------------------------------------------------------------------------------
-void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 5.0)
+//void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
+void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
 {
-int tmp3d = show_3d, tmpdata = show_data;
+  //int tmp3d = show_3d, tmpdata = show_data;
-double step[steps], acc[steps];
+  double step[steps], acc[steps];
-        show_3d = 0; show_data = 0;
+        //show_3d = 0; show_data = 0;
         //printf(" theta | acceptance \n");
         PrintGuideHead();
         for(int i=0; i<=steps; i++) {
                 Init();
                 step[i] = min + i*(max - min)/(double)steps;
                 CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
-                acc[i] = RandYZ(detector, parameters, NN, step[i], phi);
+                //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
+                acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
                 //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
                 PrintGuideStat(acc[i]);
-                //delete detector;
+                delete detector;
+        }
         //char sbuff[256];
         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
                 //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
         DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
-        show_3d = tmp3d; show_data = tmpdata;
+        //show_3d = tmp3d; show_data = tmpdata;
 }
 //------------------------------------------------------------------------------------------
+void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
+{
+  int tmp3d = show_3d, tmpdata = show_data;
+  show_3d = 1; show_data = 0;
+  //PrintGuideHead();
+  TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
+  double min = 0.0;
+  printf("\nWait, this takes a while ");
+        for(int i=0; i<=steps; i++) {
+          double theta = min + i*(maxTheta - min) / (double)steps;
+          for (int j=0; j<=steps; j++) {
+            Init();
+                  double phi = min + j*(maxPhi - min) / (double)steps;
+                  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
+      double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
+      //PrintGuideStat(acc);
+      hAcceptance->Fill(i, j, acc);
+      //printf("Acc: %f ", acc);
+      delete detector;
+      }
+    printf(".");
+  }
+  printf("\n");
+  //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
+  show_3d = tmp3d; show_data = tmpdata;
+  TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
+  cp->cd();
+  hAcceptance->Draw("COLZ");
+}
 // a vs. d
 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)
+{
-  int tmp3d = show_3d;
-        int tmpdata = show_data;
         //char sbuff[256];
-        show_3d = 0; // don't show simulations
+        //show_3d = 0; // don't show simulations
-        show_data = 1;
+        //show_data = 1;
         //double d  = detector->GetD();
         const double b = parameters.getB(); // upper side of LG
         //const double SiPM = 3.0; // the length of the detector itself
         //double reflectivity = detector->guide->getR();
         TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
         // Use the Fresnel eq. instead of fixed reflectivity 96%
         //detector->SetFresnel(1);
-        printf("   d |   a  | Acceptance\n");
+        //printf("   d |   a  | Acceptance\n");
+        getParameters();
+        printf("Wait, this takes a while ");
         for(int i=0; i<steps; i++) {
                 const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
                 for(int j=0; j<steps; j++) {
+                  Init();
                   const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
-                  const double M = b/a;
+                  //const double M = b/a;
-      parameters.setGuide(a, M, d);
+      parameters.setGuide(a, b, d);
                   CDetector *detector = new CDetector(CENTER, parameters);
                   //detector->guide->setLG(y, M, x);
                   //Init(); exclude simulation
-                  double acceptance = RandIso(detector, parameters, NN, theta, 0, 0);
+                  double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
                   //double acceptance = Grid(NN, theta);
                   //double acceptance = -1.0;
                   hAcceptance->Fill(d, a, acceptance);
-                  //printf("%.2lf | %.2lf | ", x,y);
+                  //printf("%.2lf | %.2lf | ", d, a);
-                  //PrintGuideStat(izkoristek);
+                  //PrintGuideStat(acceptance);
                   delete detector; //works fine, 50x50 grid takes ~4MB of RAM
+                  }
+                printf(".");
         }
+        printf("\n");
-        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
+        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
+        cp->cd();
-        TVirtualPad *pacc = cp->cd(0);
+        //TVirtualPad *pacc = cp->cd(0);
+        /*
         pacc->SetRightMargin(0.10);
         pacc->SetLeftMargin(0.10);
         pacc->SetTopMargin(0.10);
+        pacc->SetBottomMargin(0.10);
+        */
         TFile *file = new TFile("acceptance.root","RECREATE");
         hAcceptance->Write();
         file->Close();
         //delete file;
         //hAcceptance->SetContour(100);
         //gStyle->SetPalette(1,0);
         hAcceptance->SetTitle(";d [mm];a [mm]");
         hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
         hAcceptance->Draw("COLZ");
+        char filename[128];
+        sprintf(filename,"LGI_ad%d.C", steps);
+        cp->SaveAs(filename);
-        show_3d = tmp3d; show_data = tmpdata;
+}
 //------------------------------------------------------------------------------------------
 // collection efficiency vs length
 // a changes, b rests the same, d changes
         PrintGuideHead();
         for(int i=0; i<=steps; i++) {
                 step[i] = min + i*(max - min)/(double)steps;
                 sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
-                parameters.setGuide(sipm_d, M0/sipm_d, step[i]);
+                parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
                 //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
                 printf("%.1lf | %.2lf | ", step[i], sipm_d);
                 Init();
-                acc[i] = RandIso(detector, parameters, NN, theta);
+                acc[i] = RandIso(detector, parameters, NN, 0, theta);
                 PrintGuideStat(acc[i]);
                 sprintf(sbuff, "d%2d.gif", i);
                 //c3dview->SaveAs(sbuff);
+        }
         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
         //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
         TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
         TVirtualPad *pacc = cp->cd(0);
         pacc->SetRightMargin(0.05);
         pacc->SetLeftMargin(0.13);
         // Use Fresnel equations
         //detector->SetFresnel(1);
         // Set glass (n=1.5) at the exit
         //detector->SetGlass(1,0);
         PrintGuideHead();
         for(int i=0; i<=steps; i++) {
                 step[i] = min + i*(max - min)/(double)steps;
+                //parameters.setGuide(a, magnif, step[i]);
                 parameters.setGuide(a, magnif, step[i]);
                 CDetector *detector = new CDetector(CENTER, parameters);
                 //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
                 printf("%.1lf | ", step[i]);
                 Init();
-                acc[i] = RandIso(detector, parameters, NN, theta);
+                acc[i] = RandIso(detector, parameters, NN, 0, theta);
                 PrintGuideStat(acc[i]);
                 delete detector;
         }
         //char sbuff[256];
         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
         DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
         show_3d = tmp3d; show_data = tmpdata;
 }
 //------------------------------------------------------------------------------------------
 // Magnification optimization. a and d are fixed, M and b change in steps
 void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
+{
         int tmp3d = show_3d, tmpdata = show_data;
         double step[steps], acc[steps];
         const double a = parameters.getA();
         const double d = parameters.getD();
         show_3d = 0; show_data = 0;
         PrintGuideHead();
         for(int i=0; i<=steps; i++) {
                 step[i] = min + i*(max - min)/(double)steps;
+                //parameters.setGuide(a, step[i], d);
-                parameters.setGuide(a, step[i], d);
+                parameters.setGuide(a, a*step[i], d);
                 CDetector *detector = new CDetector(CENTER, parameters);
                 Init();
-                acc[i] = RandIso(detector, parameters, NN, theta);
+                acc[i] = RandIso(detector, parameters, NN, 0, theta);
                 PrintGuideStat(acc[i]);
+        }
         //char sbuff[256];
         //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
                 //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
                         //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
         DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
         show_3d = tmp3d; show_data = tmpdata;
 }
 //------------------------------------------------------------------------------------------
+void getParameters()
+{
+  printf("LIGHT GUIDE\n"
+    "   b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
+  printf("MATERIAL REFRACITVE INDICES\n"
+    "   n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
+  printf("PLATE\n"
+    "   ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
+  return;
+}

Subversion Repositories f9daq

(root)/lightguide/trunk/src/userFunctions.cpp – Rev 54 → 70