CalcFiberLoopCorr.C
Go to the documentation of this file.
1 #include <TRandom3.h>
2 #include <TH1D.h>
3 #include <TH3D.h>
4 #include <TMath.h>
5 #include <TFile.h>
6 
7 #include <cmath>
8 #include <iostream>
9 #include <iomanip>
10 
11 using namespace std;
12 using namespace TMath;
13 
14 bool inTorus(double x, double y, double z, double r, double R)
15 {
16  if (x > 0)
17  {
18  //double det = pow( x*x + y*y + R*R - r*r, 2) - 4.0*R*R*(x*x + y*y);
19  double det = pow( R - sqrt(x*x + y*y),2) + z*z - r*r;
20  if (det <= 0.0) return true;
21  else return false;
22  }
23  else
24  {
25  double det1 = pow(y - R,2) + z*z - r*r;
26  double det2 = pow(y + R,2) + z*z - r*r;
27  if (det1 <= 0.0 || det2 <= 0.0) return true;
28  else return false;
29  }
30 }
31 
33 {
34  //ring diameter = 61.5 mm
35  //groove 1.6 mm deep
36  //fiber radius = 0.7 mm
37 
38  double ringD = 6.15 - 0.16;
39  double r = 0.07;
40 
41  double R = 0.5*ringD+r;
42 
43  TRandom3 gen;
44  gen.SetSeed(0);
45 
46  double xmin = -2*R;
47  //double xmin = 0;
48  double xmax = (R+r);
49 
50  double ymin = (-R-r);
51  double ymax = (R+r);
52 
53  double zmin = (-r);
54  double zmax = (r);
55 
56  double Abox = (ymax-ymin)*(zmax-zmin);
57  double Afiber = 2.0*Pi()*r*r;
58 
59  TFile* outfile = new TFile("FiberLoopCorr.root", "recreate");
60  TH1D* histo = new TH1D("hFiberLoopCorr", ";Distance from start of the bend;Correction Factor", 100, xmin, xmax);
61  histo->SetDirectory(outfile);
62  histo->SetStats(false);
63  histo->GetXaxis()->CenterTitle();
64  histo->GetYaxis()->CenterTitle();
65 
66  double m = R+r;
67 
68  //TH3D* hBox = new TH3D("hBox", "", 50, -m, m, 50, -m, m, 50, -m, m);
69 
70  int ntot = 10000000;
71  /*
72  for (int iter = 0; iter < ntot; ++iter)
73  {
74  double x = gen.Uniform(-m, m);
75  double y = gen.Uniform(-m, m);
76  double z = gen.Uniform(-m, m);
77  if (inTorus(x, y, z, r, R)) hBox->Fill(x, y, z);
78  }
79  */
80 
81  for (int iX = 1; iX <= histo->GetXaxis()->GetNbins(); ++iX)
82  {
83  cout << "iX = " << iX << endl;
84  double successes = 0.0;
85  double x = histo->GetBinCenter(iX);
86  for (int iter = 0; iter < ntot; ++iter)
87  {
88  //if (iter % 2000 == 0) cout << "iter = " << iter << endl;
89  double y = gen.Uniform(ymin, ymax);
90  double z = gen.Uniform(zmin, zmax);
91  //assume retaining ring blocks about half the light...
92  if (inTorus(x, y, z, r, R))
93  {
94  //hBox->Fill(x, y, z);
95  successes += 1.0;
96  }
97  }
98  histo->SetBinContent(iX, (Abox/Afiber)*(successes/ntot));
99  }
100 
101  outfile->cd();
102  histo->Write();
103  histo->Draw();
104 
105  //TCanvas* cBox = new TCanvas();
106  //hBox->Draw();
107 
108 }
109 
110 
std::map< std::string, double > xmax
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Double_t ymax
Definition: plot.C:25
bool inTorus(double x, double y, double z, double r, double R)
#define R(x)
TH2D * histo
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
Float_t iX
TRandom3 r(0)
void CalcFiberLoopCorr()
Double_t ymin
Definition: plot.C:24
FILE * outfile
Definition: dump_event.C:13