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

 TCanvas *c3dview;
 TCanvas *cacc;
 int show_3d = 0, show_data = 0;
 int draw_width = 2;
- // calls default constructor CVodnik.cpp
+// calls default constructor CVodnik.cpp
 //void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
-        //{center.SetXYZ(x,y,z); detector = new CDetector(center);}
+//{center.SetXYZ(x,y,z); detector = new CDetector(center);}
 //void SetLGType(int in = 1, int side = 1, int out = 0)
-        //{detector->SetLGType(in, side, out);}
+//{detector->SetLGType(in, side, out);}
 //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)
-        //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
+//{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
 //void SetR(double R0) {detector->SetR(R0);}
-        /*
+/*
 void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)
         {detector->SetGlass(glass_on0, glass_d0);}
 void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
         {detector->SetGap(x_gap0 , y_gap0, z_gap0);}
 void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)
         {detector->SetRCol(in, lg, out, gla);};
 void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)
         {detector->SetDCol(LG0, glass0, active0);};
 void SetDWidth(int w = 1) {draw_width = w;}
 void SetFresnel(int b = 1) {detector->SetFresnel(b);}
 void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}
 void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}
 void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
 void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
-*/
+ */
 //-----------------------------------------------------------------------------
-void Init(double rsys_scale = 7.0)
+void Init(double rsys_scale = 10.0)
+{
-        RTSetStyle(gStyle);
+  RTSetStyle(gStyle);
-        gStyle->SetOptStat(10);
+  gStyle->SetOptStat(10);
-        if(show_3d) {
+  if(show_3d) {
-                double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
+      double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
-                double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
+      double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
-                c3dview = (TCanvas*)gROOT->FindObject("c3dview");
+      c3dview = (TCanvas*)gROOT->FindObject("c3dview");
-                if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);      c3dview->SetFillColor(0);}
+      if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
+      c3dview->SetFillColor(0);}
-                else c3dview->Clear();
+      else c3dview->Clear();
-                TView3D *view = new TView3D(1, rsys0, rsys1);
+      TView3D *view = new TView3D(1, rsys0, rsys1);
-                view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
+      view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
-                //view->SetPerspective();
+      //view->SetPerspective();
-                view->SetParallel();
+      view->SetParallel();
-                view->FrontView();
+      view->FrontView();
-                view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
+      view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
-                //view->ToggleRulers();
+      //view->ToggleRulers();
-        }
+  }
+}
 //-----------------------------------------------------------------------------
 void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
+{
-        if(show_data) {
+  if(show_data) {
-                char sbuff[256];
+      char sbuff[256];
-                sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
+      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.getA(), parameters.getLightYield()*acc, parameters.getD(),
-                                parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
+          parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
-                if(!only2d) {
+      if(!only2d) {
-                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
+          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
-                        cdata->Divide(3,3);
+          cdata->Divide(3,3);
-                        int cpc = 1;
+          int cpc = 1;
-                        cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
+          cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
-                        cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
+          cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
-                        cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
+          cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
-                        (detector->GetHGlass())->Draw("COLZ");
+          (detector->GetHGlass())->Draw("COLZ");
-                        cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
+          cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
-                        cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
+          cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
-                        cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
+          cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
-                        cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
+          cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
-                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
+          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
-                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
+          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
-                }else{
+      }else{
-                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
+          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
-                        cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
+          cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
-                        cdata->SaveAs("cdata.gif");
+          cdata->SaveAs("cdata.gif");
-                }
+      }
-        }
+  }
+}
 //-----------------------------------------------------------------------------
 void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
-                         double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
+    double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
+{
-        cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
+  cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
-        (cacc->cd(0))->SetRightMargin(0.05);
+  (cacc->cd(0))->SetRightMargin(0.05);
-        (cacc->cd(0))->SetLeftMargin(0.13);
+  (cacc->cd(0))->SetLeftMargin(0.13);
-        (cacc->cd(0))->SetTopMargin(0.08);
+  (cacc->cd(0))->SetTopMargin(0.08);
-        TGraph *gacceptance = new TGraph(n, x, y);
+  TGraph *gacceptance = new TGraph(n, x, y);
-        gacceptance->SetTitle(htitle);
+  gacceptance->SetTitle(htitle);
-        gacceptance->SetMarkerStyle(8);
+  gacceptance->SetMarkerStyle(8);
-        gacceptance->SetLineColor(12);
+  gacceptance->SetLineColor(12);
-        gacceptance->SetLineWidth(1);
+  gacceptance->SetLineWidth(1);
-        if(xmin < xmax)
+  if(xmin < xmax)
-                (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
+    (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
-        (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
+  (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
-        gacceptance->Draw("ACP");
+  gacceptance->Draw("ACP");
-        //cacc->Print("acceptance.pdf");
+  //cacc->Print("acceptance.pdf");
-        //delete gacceptance;
+  //delete gacceptance;
-        //delete cacc;
+  //delete cacc;
+}
 //-----------------------------------------------------------------------------
 void PrintGlassStat(CDetector *detector)
+{
-        int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
+  int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
-        int glass_out = (detector->GetHGlass())->GetEntries();
+  int glass_out = (detector->GetHGlass())->GetEntries();
-        printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
+  printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
+}
 //-----------------------------------------------------------------------------
 void PrintGuideHead()
+{
-        //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
+  //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
-        printf("Acceptance\n");
+  printf("Acceptance\n");
+}
 //-----------------------------------------------------------------------------
 void PrintGuideStat(double acceptance)
+{
-        /*int n_active = (detector->GetHActive())->GetEntries();
+  /*int n_active = (detector->GetHActive())->GetEntries();
         int n_enter = (detector->GetLG())->GetEnteranceHits();
         double izkoristek = n_active/(double)n_enter;
         int fates[7]; (detector->GetLG())->GetVFate(fates);
 //      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
         printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
                    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
         for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);
         printf("%5.1lf\n", 100.0*izkoristek);*/
-        //double n_active = (detector->GetHActive())->GetEntries();
+  //double n_active = (detector->GetHActive())->GetEntries();
-        //double r_acc = 100.0 * n_active / n_in;
+  //double r_acc = 100.0 * n_active / n_in;
-        double r_acc = 100.0 * acceptance;
+  double r_acc = 100.0 * acceptance;
-        printf("%7.3lf\n", r_acc);
+  printf("%7.3lf\n", r_acc);
 }
 //-----------------------------------------------------------------------------
 //-----------------------------------------------------------------------------
 // en zarek
+double Single(CDetector *detector, DetectorParameters& parameters,
-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)
+{
-  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
-        theta = Pi()*theta/180.0;
+  theta = Pi()*theta/180.0;
-        if(theta < 1e-6) theta = 1e-6;
+  if(theta < 1e-6) theta = 1e-6;
-        phi = phi*Pi()/180.0;
+  phi = phi*Pi()/180.0;
-        if(phi < 1e-6) phi = 1e-6;
+  if(phi < 1e-6) phi = 1e-6;
-        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
+  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
-        //detector->Init();
+  //detector->Init();
-        if(show_3d) detector->Draw(draw_width);
+  if(show_3d) detector->Draw(draw_width);
-        CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
+  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));
+      TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
+      TMath::Sin(theta)*TMath::Cos(phi));
+  // 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 = new CRay;
+  CRay *ray1 = new CRay;
-        detector->Propagate(*ray0, ray1, show_3d);
+  detector->Propagate(*ray0, ray1, show_3d);
+  CRay *incidentPolarization = new CRay;
+  incidentPolarization->SetColor(kGreen);
+  TVector3 drawPosition = ray0->GetR();
+  drawPosition.SetX(drawPosition.X() - 5);
+  incidentPolarization->Set(drawPosition, polarization);
+  incidentPolarization->DrawS(drawPosition.X(), 1);
+  TVector3 outPolarization = ray1->GetP();
+  drawPosition = ray1->GetR();
+  drawPosition.SetX(drawPosition.X() + 5);
+  CRay* rayPol = new CRay(drawPosition, outPolarization);
+  rayPol->SetColor(kBlack);
+  rayPol->DrawS(drawPosition.X(), 1);
-        delete ray0;
+  delete ray0;
-        delete ray1;
+  delete ray1;
-        return (detector->GetHActive())->GetEntries() / (double)(1);
+  return (detector->GetHActive())->GetEntries() / (double)(1);
+}
 //-----------------------------------------------------------------------------
 // zarki, razporejeni v mrezi
-double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
+double Grid(CDetector *detector, DetectorParameters& parameters,
+int NN = 10, double theta = 0.0)
 {
-//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
+  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
-        theta = Pi()*theta/180.0;
+  theta = Pi()*theta/180.0;
-        if(theta < 1e-6) theta = 1e-6;
+  if(theta < 1e-6) theta = 1e-6;
-        //detector->Init();
+  //detector->Init();
-        if(show_3d) detector->Draw(draw_width);
+  if(show_3d) detector->Draw(draw_width);
-        const double b = parameters.getB();
+  const double b = parameters.getB();
-        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
+  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
-        for(int i = 0; i < NN; i++) {
+  for(int i = 0; i < NN; i++) {
-                for(int j = 0; j < NN; j++) {
+      for(int j = 0; j < NN; j++) {
-                CRay *ray0 = new CRay(CENTER.x() - offset,
+          CRay *ray0 = new CRay(CENTER.x() - offset,
 .99*(i*b/NN + b/(2.0*NN) - b/2.0),
 .99*(j*b/NN + b/(2.0*NN) - b/2.0),
-                                          TMath::Cos(theta),
+              TMath::Cos(theta),
 .0,
-                                          TMath::Sin(theta));
+              TMath::Sin(theta));
+          TVector3 polarization(0, 1, 0);
+          polarization.RotateY(-theta);
+          ray0->setPolarization(polarization);
-                CRay *ray1 = new CRay;
+          CRay *ray1 = new CRay;
-                if(i == (NN/2))
+          if(i == (NN/2))
-                        detector->Propagate(*ray0, ray1, show_3d);
+            detector->Propagate(*ray0, ray1, show_3d);
-                else
+          else
-                        detector->Propagate(*ray0, ray1, 0);
+            detector->Propagate(*ray0, ray1, 0);
-                delete ray0;
+          delete ray0;
-                delete ray1;
+          delete ray1;
-                }
+      }
-        }
+  }
-        double acceptance = 0.0;
+  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));
+  acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
-        return acceptance;
+  return acceptance;
+}
 //-----------------------------------------------------------------------------
 // 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)
+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);
+  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
-        theta = theta*3.14159265358979312/180.0;
+  theta = theta*3.14159265358979312/180.0;
-        if(theta < MARGIN) theta = MARGIN;
+  if(theta < MARGIN) theta = MARGIN;
-        phi = phi*3.14159265358979312/180.0;
+  phi = phi*3.14159265358979312/180.0;
-        if(phi < MARGIN) phi = MARGIN;
+  if(phi < MARGIN) phi = MARGIN;
-        TDatime now;
+  TDatime now;
-        TRandom rand;
+  TRandom rand;
-        rand.SetSeed(now.Get());
+  rand.SetSeed(now.Get());
-        //detector->Init();
+  //detector->Init();
-        if(show_3d) detector->Draw(draw_width);
+  if(show_3d) detector->Draw(draw_width);
-        double SiPM = parameters.getA();
+  double SiPM = parameters.getA();
-        double    M = parameters.getM();
+  double    M = parameters.getM();
-        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
+  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
-        for(int i = 0; i < NN; i++) {
+  for(int i = 0; i < NN; i++) {
-                double rx = CENTER.x() - offset;
+      double rx = CENTER.x() - offset;
-                double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/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 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
-                double rl = TMath::Cos(theta);
+      double rl = TMath::Cos(theta);
-                double rm = TMath::Sin(theta)*TMath::Sin(phi);
+      double rm = TMath::Sin(theta)*TMath::Sin(phi);
-                double rn = TMath::Sin(theta)*TMath::Cos(phi);
+      double rn = TMath::Sin(theta)*TMath::Cos(phi);
-                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
+      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
+      TVector3 polarization(0, 0, 1);
+      polarization.RotateY(-theta);
+      polarization.RotateZ(phi);
+      ray0->setPolarization(polarization);
-                CRay *ray1 = new CRay;
+      CRay *ray1 = new CRay;
-                if(i < show_rays)
+      if(i < show_rays)
-                        detector->Propagate(*ray0, ray1, show_3d);
+        detector->Propagate(*ray0, ray1, show_3d);
-                else
+      else
-                        detector->Propagate(*ray0, ray1, 0);
+        detector->Propagate(*ray0, ray1, 0);
-                //delete ray0;
+      delete ray1;
-                //delete ray1;
+      delete ray0;
-        }
+  }
-        double acceptance = 0.0;
+  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;
+  acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
-        return acceptance;
+  return acceptance;
+}
 //-----------------------------------------------------------------------------
 // zarki, izotropno porazdeljeni znotraj kota theta
 // = nakljucni vstopni polozaj in kot
 // 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,
-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)
+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);
+  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
-        double pi = 3.14159265358979312;
+  double pi = 3.14159265358979312;
-        theta = theta*3.14159265358979312/180.0;
+  theta = theta*3.14159265358979312/180.0;
-        if(theta < MARGIN) theta = MARGIN;
+  if(theta < MARGIN) theta = MARGIN;
-        TDatime now;
+  TDatime now;
-        TRandom rand;
+  TRandom rand;
-        rand.SetSeed(now.Get());
+  rand.SetSeed(now.Get());
-        double rx, ry, rz, rl, rm, rn;
+  double rx, ry, rz, rl, rm, rn;
-        double rphi, rtheta;
+  double rphi, rtheta;
-        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
+  //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
-        double theta_min_rad = TMath::Cos(theta);
+  double theta_min_rad = TMath::Cos(theta);
-        double theta_max_rad = TMath::Cos(thetaMin);
+  double theta_max_rad = TMath::Cos(thetaMin);
-          TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
+  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
-          hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
+  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
-          hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
+  hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
-          htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
+  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
-          htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
+  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
-          hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
+  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
-          hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
+  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
-          hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
+  hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
-          hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
+  hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
-          hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
+  hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
-          hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
+  hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
-          hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
+  hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
-          hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
+  hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
-        //detector->Init();
+  //detector->Init();
-        if (show_3d) detector->Draw(draw_width);
+  if (show_3d) detector->Draw(draw_width);
-        double SiPM = parameters.getA();
+  double SiPM = parameters.getA();
-        double M = parameters.getM();
+  double M = parameters.getM();
-        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
+  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
-        for(int i = 0; i < NN; i++) {
+  for(int i = 0; i < NN; i++) {
-                rx = CENTER.x() - offset;
+      rx = CENTER.x() - offset;
-                ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
-                rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
-                rphi = 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);
+      //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, 1.0)));
-                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
+      rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
-                rl = TMath::Cos(rtheta);
+      rl = TMath::Cos(rtheta);
-                rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
+      rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
-                rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
+      rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
-                if(show_rand) {
+      if(show_rand) {
-                        htheta->Fill(rtheta);
+          htheta->Fill(rtheta);
-                        hcostheta->Fill( TMath::Cos(rtheta) );
+          hcostheta->Fill( TMath::Cos(rtheta) );
-                        hphi->Fill(rphi);
+          hphi->Fill(rphi);
-                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
+          hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
-                }
+      }
-                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
+      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
-                CRay *ray1 = new CRay;
+      CRay *ray1 = new CRay;
-                if(i < show_rays) {
+      if(i < show_rays) {
-                        detector->Propagate(*ray0, ray1, show_3d);
+          detector->Propagate(*ray0, ray1, show_3d);
-                        }
+      }
-                else {
+      else {
-                        detector->Propagate(*ray0, ray1, 0);
+          detector->Propagate(*ray0, ray1, 0);
-                        }
+      }
-          delete ray0;
+      delete ray0;
-          delete ray1;
+      delete ray1;
-        }
+  }
-        if(show_rand) {
+  if(show_rand) {
-                TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
+      TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
-                if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
+      if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
-                else c2rand->Clear();
+      else c2rand->Clear();
-                c2rand->Divide(2,3);
+      c2rand->Divide(2,3);
-                c2rand->cd(1); hphi->Draw();
+      c2rand->cd(1); hphi->Draw();
-                c2rand->cd(3); htheta->Draw();
+      c2rand->cd(3); htheta->Draw();
-                c2rand->cd(5); hcostheta->Draw();
+      c2rand->cd(5); hcostheta->Draw();
-                c2rand->cd(2); hl->Draw();
+      c2rand->cd(2); hl->Draw();
-                c2rand->cd(4); hm->Draw();
+      c2rand->cd(4); hm->Draw();
-                c2rand->cd(6); hn->Draw();
+      c2rand->cd(6); hn->Draw();
-        }
+  }
-        double acceptance = 0.0;
+  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;
+  return acceptance;
+}
 // Beamtest distribution
 // with fixed theta and phi in interval phiMin, phiMax
+double beamtest(CDetector *detector, DetectorParameters& parameters,
-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;
+  double pi = 3.14159265358979312;
-        theta *= pi/180.0;
+  theta *= pi/180.0;
-        if(theta < MARGIN) theta = MARGIN;
+  if(theta < MARGIN) theta = MARGIN;
-        phiMin *= pi/180.0;
+  phiMin *= pi/180.0;
-        phiMax *= pi/180.0;
+  phiMax *= pi/180.0;
-        TDatime now;
+  TDatime now;
-        TRandom rand;
+  TRandom rand;
-        rand.SetSeed(now.Get());
+  rand.SetSeed(now.Get());
-        double rx, ry, rz, rl, rm, rn;
+  double rx, ry, rz, rl, rm, rn;
-        double rphi;
+  double rphi;
+  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);
-        //double rtheta;
+  //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_min_rad = TMath::Power(TMath::Cos(theta), 2);
+      //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)));
-        //double theta_min_rad = TMath::Cos(theta);
+      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);
+      // inital polarizaton
+      TVector3 polarization(0, 0, 1);
+      polarization.RotateY(-theta);
+      polarization.RotateX(rphi);
+      ray0->setPolarization(polarization);
+      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);
-        TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
-        hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
-        hphi = new TH1F("hphi", "hphi", 100, -pi, pi);
+      c2rand->cd(1); hphi->Draw();
-        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;
+      c2rand->cd(3); htheta->Draw();
-        hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
+      c2rand->cd(5); hcostheta->Draw();
-        hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
+      c2rand->cd(2); hl->Draw();
-        hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
+      c2rand->cd(4); hm->Draw();
-        hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
-        hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
+      c2rand->cd(6); hn->Draw();
+  }
-        //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;
+  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;
+  return acceptance;
+}

Subversion Repositories f9daq

(root)/lightguide/trunk/src/raySimulator.cpp – Rev 71 → 73