Subversion Repositories f9daq

// 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++);

          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");

          TH2F* histoPlate = detector->GetHPlate();

          histoPlate->Draw("COLZ");

          cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");

          TH2F* histoGlass = (detector->GetHGlass());

          histoGlass->Draw("COLZ");

          cdata->cd(cpc++);

          TH2F* histoActive = (detector->GetHActive());

          histoActive->Draw("COLZ");

          cdata->cd(cpc++);

          TH2F* histoLaser = (detector->GetHLaser());

          histoLaser->Draw("COLZ");

          cdata->cd(cpc++);

          TH2F* histoDetector = (detector->GetHDetector());

          histoDetector->Draw("COLZ");

          cdata->cd(cpc++);

          TH1F* fate=((detector->GetLG())->GetHFate());

          fate->Draw();

          cdata->cd(cpc++);

          TH1F* reflections = ((detector->GetLG())->GetHNOdb_all());

          reflections->Draw();

          cdata->cd(cpc++);

          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");

      }

  }

}

//-----------------------------------------------------------------------------

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 y-polarization == vertical

  TVector3 k(ray0->GetK());

  TVector3 polarization(k.Orthogonal());

  polarization.Unit();

  polarization.Rotate(phi, k);

  if (polarization.Dot(ray0->GetK()) > 1e-2) printf("ERROR: pol not perep\n");

  ray0->setPolarization(polarization);

  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();

  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 simulated = 0;

  double offset = (parameters._plateOn ? parameters._plateWidth : 0);

  const int GRID = 100;

  for(int i = 0; i < GRID; i++) {

      for(int j = 0; j < GRID; j++) {

          for(int jk = 0; jk<NN; jk++){

              CRay *ray0 = new CRay(CENTER.x() - offset,

                  (i*b/GRID + b/(2.0*GRID) - b/2.0),

                  (j*b/GRID + b/(2.0*GRID) - b/2.0),

                  TMath::Cos(theta),

                  0.0,

                  TMath::Sin(theta));

              TVector3 k(ray0->GetK());

              TVector3 polarization(k.Orthogonal());

              polarization.Unit();

              ray0->setPolarization(polarization);

              CRay ray1;

              if(i == (GRID/2) and jk == 0)

                detector->Propagate(*ray0, ray1, show_3d);

              else

                detector->Propagate(*ray0, ray1, 0);

              delete ray0;

              simulated++;

              //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() / 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)

{

  //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;

  TRandom2 rand;

  rand.SetSeed(now.Get());

  //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++) {

      double rx = CENTER.x() - offset;

      double ry = rand.Uniform(-b/2.0, b/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(k.Orthogonal());

      polarization.Unit();

      polarization.Rotate(phi, k);

      ray0->setPolarization(polarization);

      CRay ray1;

      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 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;

  thetaMin = thetaMin*pi/180.0;

  if(thetaMin < MARGIN) thetaMin = MARGIN;

  if(thetaMax < MARGIN) thetaMax = MARGIN;

  TDatime now;

  TRandom2 rand;

  rand.SetSeed(now.Get());

  //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);

  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);

  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++) {

      double rx = CENTER.x() - offset;

      double ry = rand.Uniform(-b/2.0, b/2.0);

      double rz = rand.Uniform(-b/2.0, b/2.0);

      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)));

      double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));

      double rl = TMath::Cos(theta);

      double rm = TMath::Sin(theta)*TMath::Sin(phi);

      double rn = TMath::Sin(theta)*TMath::Cos(phi);

      if(show_rand) {

          htheta->Fill(theta);

          hcostheta->Fill( TMath::Cos(theta) );

          hphi->Fill(phi);

          hl->Fill(rl);

          hm->Fill(rm);

          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;

      if(i < show_rays) {

          detector->Propagate(*ray0, ray1, show_3d);

      }

      else {

          detector->Propagate(*ray0, ray1, 0);

      }

      //delete ray1;

      delete ray0;

  }

  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;

  TRandom2 rand;

  rand.SetSeed(now.Get());

  double rx, ry, rz, rl, rm, rn;

  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;

  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);

      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(phi);

      rn = TMath::Sin(theta)*TMath::Cos(phi);

      if(show_rand) {

          //htheta->Fill(rtheta);

          //hcostheta->Fill( TMath::Cos(rtheta) );

          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(k.Orthogonal());

      polarization.Unit();

      polarization.Rotate(phi, k);

      ray0->setPolarization(polarization);

      CRay ray1;

      if(i < show_rays) {

          detector->Propagate(*ray0, ray1, show_3d);

      }

      else {

          detector->Propagate(*ray0, ray1, 0);

      }

      //delete ray1;

      delete ray0;

  }

  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;

}