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


// own include
//#include "include/raySimulator.h"
#include "include/RTUtil.h"
#include "include/guide.h"
 
// general include
#include "TROOT.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TView3D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TBenchmark.h"
#include "TPolyMarker.h"
#include "TGraph.h"
#include "TF1.h"
 
//=================================================================================
TCanvas *c3dview;
TCanvas *cacc;
int show_3d = 0, show_data = 0;
int draw_width = 2;
 
// 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);}
 
 
//void SetLGType(int in = 1, int side = 1, int out = 0)
//{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);}
 
//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 = 10.0)
{      
  RTSetStyle(gStyle);
  gStyle->SetOptStat(10);
 
  if(show_3d) {
      double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
      double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
 
      c3dview = (TCanvas*)gROOT->FindObject("c3dview");
      if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
      c3dview->SetFillColor(0);}
      else c3dview->Clear();
 
      TView3D *view = new TView3D(1, rsys0, rsys1);
      view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
      //view->SetPerspective();
      view->SetParallel();
      view->FrontView();
      view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
      //view->ToggleRulers();
  }
}
//-----------------------------------------------------------------------------
void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
{
  if(show_data) {
      char sbuff[256];
      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.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
 
      if(!only2d) {
          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
          cdata->Divide(3,3);
          int cpc = 1;
          cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
          cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
          cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
          (detector->GetHGlass())->Draw("COLZ");
 
          cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
          cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
          cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
 
          cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
          cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
      }else{
          RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
          cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
          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)
{
 
  cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
 
  (cacc->cd(0))->SetRightMargin(0.05);
  (cacc->cd(0))->SetLeftMargin(0.13);
  (cacc->cd(0))->SetTopMargin(0.08);
 
  TGraph *gacceptance = new TGraph(n, x, y);
  gacceptance->SetTitle(htitle);
  gacceptance->SetMarkerStyle(8);
  gacceptance->SetLineColor(12);
  gacceptance->SetLineWidth(1);
  if(xmin < xmax)
    (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
  (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
 
  gacceptance->Draw("ACP");
  //cacc->Print("acceptance.pdf");
  //delete gacceptance;
  //delete cacc;
}
//-----------------------------------------------------------------------------
void PrintGlassStat(CDetector *detector)
{
  int glass_in = ((detector->GetLG())->GetHOut())->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);
}
//-----------------------------------------------------------------------------
void PrintGuideHead()
{
  //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
  printf("Acceptance\n");
}
//-----------------------------------------------------------------------------
void PrintGuideStat(double acceptance)
{
  /*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 r_acc = 100.0 * n_active / n_in;
  double r_acc = 100.0 * acceptance;
  printf("%7.3lf\n", r_acc);
}
//-----------------------------------------------------------------------------
 
//-----------------------------------------------------------------------------
// 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)
{
 
  theta = Pi()*theta/180.0;
  if(theta < 1e-6) theta = 1e-6;
  phi = phi*Pi()/180.0;
  if(phi < 1e-6) phi = 1e-6;
 
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 
  //detector->Init();
  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
  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;
 
  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 ray1;
 
  return (detector->GetHActive())->GetEntries() / (double)(1);
}
//-----------------------------------------------------------------------------
// zarki, razporejeni v mrezi
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);
  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 offset = (parameters._plateOn ? parameters._plateWidth : 0);
 
  for(int i = 0; i < NN; i++) {
      for(int j = 0; j < NN; j++) {
          CRay *ray0 = new CRay(CENTER.x() - offset,
              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
              TMath::Cos(theta),
              0.0,
              TMath::Sin(theta));
          TVector3 polarization(0, 1, 0);
          polarization.RotateY(-theta);
          ray0->setPolarization(polarization);
          CRay *ray1 = new CRay;
          if(i == (NN/2))
            detector->Propagate(*ray0, ray1, show_3d);
          else
            detector->Propagate(*ray0, ray1, 0);
          delete ray0;
          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));
  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)
{
  //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;
  rand.SetSeed(now.Get());
 
  //detector->Init();
  if(show_3d) detector->Draw(draw_width);
 
  double SiPM = parameters.getA();
  double    M = parameters.getM();
  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 rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/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 polarization(0, 0, 1);
      polarization.RotateY(-theta);
      polarization.RotateZ(phi);
      ray0->setPolarization(polarization);
      CRay *ray1 = new CRay;
 
      if(i < show_rays)
        detector->Propagate(*ray0, ray1, show_3d);
      else
        detector->Propagate(*ray0, ray1, 0);
      delete ray1;
      delete ray0;
  }
 
  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;
  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,
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);
  double pi = 3.14159265358979312;
  theta = theta*3.14159265358979312/180.0;
  if(theta < MARGIN) theta = MARGIN;
 
  TDatime now;
  TRandom 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_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);
  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);
 
 
  //detector->Init();
  if (show_3d) detector->Draw(draw_width);
 
  double SiPM = parameters.getA();
  double M = parameters.getM();
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 
  for(int i = 0; i < NN; i++) {
 
      rx = CENTER.x() - offset;
      ry = 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);
      //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(rtheta);
      rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
      rn = TMath::Sin(rtheta)*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;
  /*
        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;
}
 
// 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)
{
  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;
  rand.SetSeed(now.Get());
  double rx, ry, rz, rl, rm, rn;
  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);
 
  //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);
      // 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);
 
      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;
  /*
        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;
}
 

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

Subversion Repositories f9daq

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