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

       if(!only2d) {
           RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
           cdata->Divide(3,3);
           int cpc = 1;
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
+          TH2F* generated = detector->GetGenerated();
+          int nGenerated = generated->GetEntries();
+          int minimum = generated->GetBinContent(generated->GetMinimumBin());
+          int maximum = generated->GetBinContent(generated->GetMaximumBin());
+          generated->GetZaxis()->SetRangeUser(0, 1.05*maximum);
+          double variation = (maximum-minimum)/(double)nGenerated;
+          printf("Statistical variation (max-min)/all = %f perc. \n", variation*100);
+          generated->SetTitle("Generated");
+          generated->Draw("colz");
-          cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
+          cdata->cd(cpc++); //((detector->GetLG())->GetHIn())->Draw("COLZ");
+          TH2F* histoPlate = detector->GetHPlate();
+          histoPlate->Draw("COLZ");
           cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
+          TH2F* histoGlass = (detector->GetHGlass());
-          (detector->GetHGlass())->Draw("COLZ");
+          histoGlass->Draw("COLZ");
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
+          TH2F* histoActive = (detector->GetHActive());
+          histoActive->Draw("COLZ");
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
+          TH2F* histoLaser = (detector->GetHLaser());
+          histoLaser->Draw("COLZ");
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
+          TH2F* histoDetector = (detector->GetHDetector());
+          histoDetector->Draw("COLZ");
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
+          TH1F* fate=((detector->GetLG())->GetHFate());
+          fate->Draw();
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
+          TH1F* reflections = ((detector->GetLG())->GetHNOdb_all());
+          reflections->Draw();
+          cdata->cd(cpc++);
-          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
+          TH1F *reflectedExit = ((detector->GetLG())->GetHNOdb_exit());
+          reflectedExit->Draw();
       }else{
           RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
           cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
           cdata->SaveAs("cdata.gif");
+      }
 //-----------------------------------------------------------------------------
 //-----------------------------------------------------------------------------
 // en zarek
 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)
+{
   theta = Pi()*theta/180.0;
   if(theta < 1e-6) theta = 1e-6;
   phi = phi*Pi()/180.0;
   if(phi < 1e-6) phi = 1e-6;
   if(show_3d) detector->Draw(draw_width);
   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));
-  // Set z-polarization == vertical
+  // Set y-polarization == vertical
+  TVector3 k(ray0->GetK());
-  TVector3 polarization(0,0,1);
+  TVector3 polarization(k.Orthogonal());
-  polarization.RotateY(-theta);
+  polarization.Unit();
-  polarization.RotateZ(phi);
+  polarization.Rotate(phi, k);
-  if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n");
+  if (polarization.Dot(ray0->GetK()) > 1e-2) printf("ERROR: pol not perep\n");
   ray0->setPolarization(polarization);
-  CRay *ray1 = new CRay;
+  CRay ray1;
   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();
+  TVector3 outPolarization = ray1.GetP();
-  drawPosition = ray1->GetR();
+  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 ray1;
+  //delete ray1;
   return (detector->GetHActive())->GetEntries() / (double)(1);
+}
 //-----------------------------------------------------------------------------
 // zarki, razporejeni v mrezi
 double Grid(CDetector *detector, DetectorParameters& parameters,
-int NN = 10, double theta = 0.0)
+    int NN = 10, double theta = 0.0)
+{
   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
   theta = Pi()*theta/180.0;
   if(theta < 1e-6) theta = 1e-6;
   //detector->Init();
   if(show_3d) detector->Draw(draw_width);
   const double b = parameters.getB();
+  double simulated = 0;
   double offset = (parameters._plateOn ? parameters._plateWidth : 0);
+  const int GRID = 100;
-  for(int i = 0; i < NN; i++) {
+  for(int i = 0; i < GRID; i++) {
-      for(int j = 0; j < NN; j++) {
+      for(int j = 0; j < GRID; j++) {
+          for(int jk = 0; jk<NN; jk++){
-          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),
+                  (i*b/GRID + b/(2.0*GRID) - b/2.0),
-.99*(j*b/NN + b/(2.0*NN) - b/2.0),
+                  (j*b/GRID + b/(2.0*GRID) - b/2.0),
-              TMath::Cos(theta),
+                  TMath::Cos(theta),
 .0,
-              TMath::Sin(theta));
+                  TMath::Sin(theta));
+              TVector3 k(ray0->GetK());
-          TVector3 polarization(0, 1, 0);
+              TVector3 polarization(k.Orthogonal());
-          polarization.RotateY(-theta);
+              polarization.Unit();
-          ray0->setPolarization(polarization);
+              ray0->setPolarization(polarization);
-          CRay *ray1 = new CRay;
+              CRay ray1;
-          if(i == (NN/2))
+              if(i == (GRID/2) and jk == 0)
-            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;
+              simulated++;
-          delete ray1;
+              //delete ray1;
+          }
+      }
+  }
   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() / simulated;
   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)
+    int NN, double theta, double phi, int show_rays)
+{
   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
   theta = theta*3.14159265358979312/180.0;
   if(theta < MARGIN) theta = MARGIN;
   phi = phi*3.14159265358979312/180.0;
   if(phi < MARGIN) phi = MARGIN;
   TDatime now;
-  TRandom rand;
+  TRandom2 rand;
   rand.SetSeed(now.Get());
   //detector->Init();
   if(show_3d) detector->Draw(draw_width);
-  double SiPM = parameters.getA();
-  double    M = parameters.getM();
+  double b = parameters.getB();
   double offset = (parameters._plateOn ? parameters._plateWidth : 0);
   for(int i = 0; i < NN; i++) {
       double rx = CENTER.x() - offset;
-      double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      double ry = rand.Uniform(-b/2.0, b/2.0);
-      double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      double rz = rand.Uniform(-b/2.0, b/2.0);
       double rl = TMath::Cos(theta);
       double rm = TMath::Sin(theta)*TMath::Sin(phi);
       double rn = TMath::Sin(theta)*TMath::Cos(phi);
       CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
+      TVector3 k(ray0->GetK());
-      TVector3 polarization(0, 0, 1);
+      TVector3 polarization(k.Orthogonal());
-      polarization.RotateY(-theta);
+      polarization.Unit();
-      polarization.RotateZ(phi);
+      polarization.Rotate(phi, k);
       ray0->setPolarization(polarization);
-      CRay *ray1 = new CRay;
+      CRay ray1;
       if(i < show_rays)
         detector->Propagate(*ray0, ray1, show_3d);
       else
         detector->Propagate(*ray0, ray1, 0);
-      delete ray1;
+      //delete ray1;
       delete ray0;
+  }
   double acceptance = 0.0;
   /*
 // NN - number of rays to be simulated with angles theta distributed uniformly
 //      in cos theta and phi uniformly:
 //      theta [0, theta]
 //      phi [0,360]
 double RandIso(CDetector *detector, DetectorParameters& parameters,
-int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
+    int NN = 1e3, double thetaMin=0.0, double thetaMax = 30.0, int show_rays = 30, int show_rand = 0)
+{
   //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
   double pi = 3.14159265358979312;
+  thetaMax = thetaMax*pi/180.0;
-  theta = theta*3.14159265358979312/180.0;
+  thetaMin = thetaMin*pi/180.0;
+  if(thetaMin < MARGIN) thetaMin = MARGIN;
-  if(theta < MARGIN) theta = MARGIN;
+  if(thetaMax < MARGIN) thetaMax = MARGIN;
   TDatime now;
-  TRandom rand;
+  TRandom2 rand;
   rand.SetSeed(now.Get());
-  double rx, ry, rz, rl, rm, rn;
-  double rphi, rtheta;
   //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
-  double theta_min_rad = TMath::Cos(theta);
+  double theta_min_rad = TMath::Cos(thetaMax);
   double theta_max_rad = TMath::Cos(thetaMin);
   TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
   hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
   hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
   //detector->Init();
   if (show_3d) detector->Draw(draw_width);
-  double SiPM = parameters.getA();
-  double M = parameters.getM();
+  double b = parameters.getB();
   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;
+      double rx = CENTER.x() - offset;
-      ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      double ry = rand.Uniform(-b/2.0, b/2.0);
-      rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
+      double rz = rand.Uniform(-b/2.0, b/2.0);
-      rphi = rand.Uniform(0.0, 2.0*pi);
+      double phi = rand.Uniform(0.0, 2.0*pi);
       //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 = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
-      rl = TMath::Cos(rtheta);
+      double rl = TMath::Cos(theta);
-      rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
+      double rm = TMath::Sin(theta)*TMath::Sin(phi);
-      rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
+      double rn = TMath::Sin(theta)*TMath::Cos(phi);
       if(show_rand) {
-          htheta->Fill(rtheta);
+          htheta->Fill(theta);
-          hcostheta->Fill( TMath::Cos(rtheta) );
+          hcostheta->Fill( TMath::Cos(theta) );
-          hphi->Fill(rphi);
+          hphi->Fill(phi);
+          hl->Fill(rl);
+          hm->Fill(rm);
-          hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
+          hn->Fill(rn);
+      }
       CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
+      TVector3 k(ray0->GetK());
+      TVector3 polarization(k.Orthogonal());
+      polarization.Unit();
+      polarization.Rotate(phi, k);
+      ray0->setPolarization(polarization);
-      CRay *ray1 = new CRay;
+      CRay ray1;
       if(i < show_rays) {
           detector->Propagate(*ray0, ray1, show_3d);
+      }
       else {
           detector->Propagate(*ray0, ray1, 0);
+      }
+      //delete ray1;
       delete ray0;
-      delete ray1;
+  }
   if(show_rand) {
       TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
       if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
+}
 // Beamtest distribution
 // with fixed theta and phi in interval phiMin, phiMax
 double beamtest(CDetector *detector, DetectorParameters& parameters,
-int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
+    int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
+{
   double pi = 3.14159265358979312;
   theta *= pi/180.0;
   if(theta < MARGIN) theta = MARGIN;
   phiMin *= pi/180.0;
   phiMax *= pi/180.0;
   TDatime now;
-  TRandom rand;
+  TRandom2 rand;
   rand.SetSeed(now.Get());
   double rx, ry, rz, rl, rm, rn;
-  double rphi;
+  double phi;
   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;
       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);
+      phi = 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);
+      rm = TMath::Sin(theta)*TMath::Sin(phi);
-      rn = TMath::Sin(theta)*TMath::Cos(rphi);
+      rn = TMath::Sin(theta)*TMath::Cos(phi);
       if(show_rand) {
           //htheta->Fill(rtheta);
           //hcostheta->Fill( TMath::Cos(rtheta) );
-          hphi->Fill(rphi);
+          hphi->Fill(phi);
           hl->Fill(rl);
           hm->Fill(rm);
           hn->Fill(rn);
+      }
       CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
       // inital polarizaton
+      TVector3 k(ray0->GetK());
-      TVector3 polarization(0, 0, 1);
+      TVector3 polarization(k.Orthogonal());
-      polarization.RotateY(-theta);
+      polarization.Unit();
-      polarization.RotateX(rphi);
+      polarization.Rotate(phi, k);
       ray0->setPolarization(polarization);
-      CRay *ray1 = new CRay;
+      CRay ray1;
       if(i < show_rays) {
           detector->Propagate(*ray0, ray1, show_3d);
+      }
       else {
           detector->Propagate(*ray0, ray1, 0);
+      }
+      //delete ray1;
       delete ray0;
-      delete ray1;
+  }
   if(show_rand) {
       TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
       if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);

Subversion Repositories f9daq

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