Subversion Repositories f9daq

Rev

Rev 72 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 72 Rev 73
Line 16... Line 16...
16 
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false);
 16 
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false);
17 
// Print the detector parameters
 17 
// Print the detector parameters
18 
void getParameters();
 18 
void getParameters();
19 
 
 19 
 
20 
//void SetLGType(int in = 1, int side = 1, int out = 0)
 20 
//void SetLGType(int in = 1, int side = 1, int out = 0)
21 
        //{detector->SetLGType(in, side, out);}
 21 
//{detector->SetLGType(in, side, out);}
22 
 
 22 
 
23 
void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46)
 23 
void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46)
24 
        { parameters.setGuide(SiPM0, b0, d0);
 24 
{ parameters.setGuide(SiPM0, b0, d0);
25 
          parameters.setIndices(n1, n2, n3); }
 25 
parameters.setIndices(n1, n2, n3); }
26 
       
 26 
 
27 
void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
 27 
void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
28 
  { parameters.setGap(x_gap0, y_gap0, z_gap0); }
 28 
{ parameters.setGap(x_gap0, y_gap0, z_gap0); }
29 
 
 29 
 
30 
void setGlass(double glassOn, double glassD)
 30 
void setGlass(double glassOn, double glassD)
31 
  { parameters.setGlass(glassOn, glassD); };
 31 
{ parameters.setGlass(glassOn, glassD); };
32 
 
 32 
 
33 
void setPlate(int plateOn = 1, double plateWidth = 1)
 33 
void setPlate(int plateOn = 1, double plateWidth = 1)
34 
  { parameters.setPlate(plateOn, plateWidth); }
 34 
{ parameters.setPlate(plateOn, plateWidth); }
35 
 
 35 
 
36 
void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
 36 
void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
37 
 
 37 
 
38 
// Refractive Indices
 38 
// Refractive Indices
39 
// n1 - around light guide - air
 39 
// n1 - around light guide - air
40 
// n2 - light guide (plate) material - k9 glass
 40 
// n2 - light guide (plate) material - k9 glass
41 
// n3 - the material at the exit - optical grease, epoxy, air, etc.
 41 
// n3 - the material at the exit - optical grease, epoxy, air, etc.
42 
void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);}
 42 
void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);}
43 
 
 43 
 
44 
int save_ary = 0;
 44 
int save_ary = 0;
45 
//------------------------------------------------------------------------------------------
 45 
//------------------------------------------------------------------------------------------
46 
void Help()
 46 
void Help()
47 
{
 47 
{
48 
        printf("void SetCenter(x, y, z)\n");
 48 
  printf("void SetCenter(x, y, z)\n");
49 
        printf("void SetLGType(in, side, out)\n");
 49 
  printf("void SetLGType(in, side, out)\n");
50 
        printf("SURF_DUMMY 0, ");
 50 
  printf("SURF_DUMMY 0, ");
51 
        printf("SURF_REFRA 1, ");
 51 
  printf("SURF_REFRA 1, ");
52 
        printf("SURF_REFLE 2, ");
 52 
  printf("SURF_REFLE 2, ");
53 
        printf("SURF_TOTAL 3, ");
 53 
  printf("SURF_TOTAL 3, ");
54 
        printf("SURF_IMPER 4\n");
 54 
  printf("SURF_IMPER 4\n");
55 
        printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n");
 55 
  printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n");
56 
        printf("void SetR(R0)\n");
 56 
  printf("void SetR(R0)\n");
57 
        printf("void SetGlass(glass_on0, glass_d0)\n");
 57 
  printf("void SetGlass(glass_on0, glass_d0)\n");
58 
        printf("void SetGap(x_gap0, y_gap0, z_gap0)\n");
 58 
  printf("void SetGap(x_gap0, y_gap0, z_gap0)\n");
59 
       
 59 
 
60 
        printf("LGG(NN, theta)\n");
 60 
  printf("LGG(NN, theta)\n");
61 
        printf("LGR(NN, theta, phi)\n");
 61 
  printf("LGR(NN, theta, phi)\n");
62 
        printf("LGI(NN, theta)\n");
 62 
  printf("LGI(NN, theta)\n");
63 
       
 63 
 
64 
        printf("LGI_gap(NN, min, max, steps, theta)\n");
 64 
  printf("LGI_gap(NN, min, max, steps, theta)\n");
65 
        printf("LGR_gap(NN, min, max, steps, theta, phi)\n");
 65 
  printf("LGR_gap(NN, min, max, steps, theta, phi)\n");
66 
        printf("LGR_th(NN, min, max, steps, phi)\n");
 66 
  printf("LGR_th(NN, min, max, steps, phi)\n");
67 
        printf("LGG_th(NN, min, max, steps)\n");
 67 
  printf("LGG_th(NN, min, max, steps)\n");
68 
       
 68 
 
69 
        printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)");  
 69 
  printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)");
70 
}
 70 
}
71 
//------------------------------------------------------------------------------------------
 71 
//------------------------------------------------------------------------------------------
72 
int show_in_steps = 0;
 72 
int show_in_steps = 0;
73 
void SetShowInSteps(int in = 1)
 73 
void SetShowInSteps(int in = 1)
74 
{                                                                                
 74 
{                                                                                
75 
        show_in_steps = in;
 75 
  show_in_steps = in;
76 
}
 76 
}
77 
//------------------------------------------------------------------------------------------
 77 
//------------------------------------------------------------------------------------------
78 
int gs_set = 0;
 78 
int gs_set = 0;
79 
void SetGS()
 79 
void SetGS()
80 
{                                                                                
 80 
{                                                                                
81 
        const UInt_t Number = 2;
 81 
  const UInt_t Number = 2;
82 
        Double_t Red[Number]   = { 1.00, 0.00};
 82 
  Double_t Red[Number]   = { 1.00, 0.00};
83 
        Double_t Green[Number] = { 1.00, 0.00};
 83 
  Double_t Green[Number] = { 1.00, 0.00};
84 
        Double_t Blue[Number]  = { 1.00, 0.00};
 84 
  Double_t Blue[Number]  = { 1.00, 0.00};
85 
        Double_t Stops[Number] = { 0.00, 1.00};
 85 
  Double_t Stops[Number] = { 0.00, 1.00};
86 
                                                                               
 - 
 
- 
 
 86 
 
87 
        Int_t nb=50;
 87 
  Int_t nb=50;
88 
        if(!gs_set) {
 88 
  if(!gs_set) {
89 
                TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
 89 
      TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb);
90 
                gs_set = 1;
 90 
      gs_set = 1;
91 
        }
 91 
  }
92 
}
 92 
}
93 
//------------------------------------------------------------------------------------------
 93 
//------------------------------------------------------------------------------------------
94 
void Draw1(CDetector *detector, int gs = 0, double range = 0.0)
 94 
void Draw1(CDetector *detector, int gs = 0, double range = 0.0)
95 
{
 95 
{
96 
        if(gs) SetGS();
 96 
  if(gs) SetGS();
97 
        TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500);
 97 
  TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500);
98 
        gStyle->SetOptTitle(0);
 98 
  gStyle->SetOptTitle(0);
99 
        gStyle->SetOptStat(0);
 99 
  gStyle->SetOptStat(0);
100 
        cdata->cd(0);
 100 
  cdata->cd(0);
101 
                if(range > 0.1) {
 101 
  if(range > 0.1) {
102 
                        ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range);
 102 
      ((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range);
103 
                        ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range);
 103 
      ((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range);
104 
                }
 104 
  }
105 
                (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]");
 105 
  (detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]");
106 
                (detector->GetHLaser())->Draw("COLZ");
 106 
  (detector->GetHLaser())->Draw("COLZ");
107 
        cdata->SaveAs("2dscan.gif");
 107 
  cdata->SaveAs("2dscan.gif");
108 
        //gStyle->SetOptTitle(1);
 108 
  //gStyle->SetOptTitle(1);
109 
        //gStyle->SetOptStat(1);
 109 
  //gStyle->SetOptStat(1);
110 
}
 110 
}
111 
 
 111 
 
112 
//------------------------------------------------------------------------------------------
 112 
//------------------------------------------------------------------------------------------
113 
 
 - 
 
114 
 
 - 
 
115 
TVector3 p_pol0(0.0, 0.0, 1.0);
 - 
 
116 
void SetPol(double x, double y, double z)
 - 
 
117 
{p_pol0.SetXYZ(x,y,z);}
 - 
 
118 
 
 113 
 
119 
// Test function
 114 
// Test function
120 
// creates an optical boundary surface
 115 
// creates an optical boundary surface
121 
// shows the propagation of light ray
 116 
// shows the propagation of light ray
122 
// and polarization state of the incident ray (green)
 117 
// and polarization state of the incident ray (green)
123 
// and surface normal (blue)
 118 
// and surface normal (blue)
124 
void PolTest(double theta = 0.0)
 119 
void PolTest(double theta = 0.0)
125 
{
 120 
{
126 
  int p_type = 1;
 121 
  int p_type = SURF_REFRA;
127 
  double p_n1 = parameters.getN1();
 122 
  double p_n1 = parameters.getN1();
128 
  double p_n2 = parameters.getN2();
 123 
  double p_n2 = parameters.getN2();
129 
        theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
 124 
  theta = 3.141593*theta/180.0;
130 
        TVector3 p_pol;
 125 
  if(theta < 1e-6) theta = 1e-6;
131 
       
 126 
 
132 
        Init();
 127 
  Init();
133 
       
 128 
 
134 
        double cx = 0.0;
 129 
  double cx = 0.0;
135 
        TVector3 vodnik_edge[4];
 130 
  TVector3 vodnik_edge[4];
136 
        double t = 3.0;
 131 
  double t = 3.0;
137 
        vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t);
 132 
  vodnik_edge[0].SetXYZ(cx, t,-t);
138 
        vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t);
 133 
  vodnik_edge[1].SetXYZ(cx, t, t);
139 
        /*
 - 
 
140 
#define SURF_DUMMY 0
 - 
 
141 
#define SURF_REFRA 1
 - 
 
142 
#define SURF_REFLE 2
 - 
 
143 
#define SURF_TOTAL 3
 134 
  vodnik_edge[2].SetXYZ(cx,-t, t);
144 
#define SURF_IMPER 4
 135 
  vodnik_edge[3].SetXYZ(cx,-t,-t);
145 
        */
 136 
 
146 
        CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
 137 
  CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
- 
 
 138 
  surf->FlipN();
147 
        surf->SetFresnel(1);
 139 
  surf->SetFresnel(1);
148 
        surf->Draw();
 140 
  surf->Draw();
149 
       
 141 
 
150 
        CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
 142 
  CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
151 
        ray->SetColor(kBlack);
 143 
  ray->SetColor(kRed);
152 
        //p_pol = rotatey(p_pol0, -theta);
 144 
  //p_pol = rotatey(p_pol0, -theta);
- 
 
 145 
  TVector3 sE(0,1,0);
- 
 
 146 
  TVector3 pE(0,0,1);
- 
 
 147 
  TVector3 p_pol = pE;
153 
        p_pol = p_pol0; p_pol.RotateY(-theta);
 148 
  p_pol.RotateY(-theta);
154 
        //printf("p_pol = "); printv(p_pol);
 149 
  //printf("p_pol = "); printv(p_pol);
155 
        ray->SetPolarization(p_pol);
 150 
  ray->setPolarization(p_pol);
156 
       
 151 
 
157 
        ray->DrawS(cx, -5.0);
 152 
  ray->DrawS(cx, -5.0);
158 
       
 153 
 
159 
        CRay *out = new CRay; out->SetColor(kRed);
 154 
  CRay *out = new CRay;
- 
 
 155 
  out->SetColor(kBlack);
160 
        TVector3 *inters = new TVector3;
 156 
  TVector3 *inters = new TVector3;
161 
        surf->PropagateRay(*ray, out, inters);
 157 
  surf->PropagateRay(*ray, out, inters);
- 
 
 158 
  printf(" n1 = %f, n2 = %f\n", p_n1, p_n2);
162 
        //if(fate == 1) out->DrawS(cx, 5.0);
 159 
  //if(fate == 1) out->DrawS(cx, 5.0);
163 
  out->DrawS(cx, 5.0);
 160 
  out->DrawS(cx, 5.0);
164 
 
 161 
 
165 
        CRay *incidentPolarization = new CRay;
 162 
  CRay *incidentPolarization = new CRay;
166 
        incidentPolarization->SetColor(kGreen);
 163 
  incidentPolarization->SetColor(kGreen);
167 
        incidentPolarization->Set(ray->GetR(), p_pol);
 164 
  incidentPolarization->Set(ray->GetR(), p_pol);
168 
        incidentPolarization->DrawS(cx, 1.0);
 165 
  incidentPolarization->DrawS(cx, 1.0);
- 
 
 166 
  printf(" GREEN: polarization vector\n");
169 
       
 167 
 
170 
        CRay *surfaceNormal = new CRay;
 168 
  CRay *surfaceNormal = new CRay;
171 
        surfaceNormal->SetColor(kBlue);
 169 
  surfaceNormal->SetColor(kBlue);
172 
        surfaceNormal->Set(ray->GetR(), surf->GetN());
 170 
  surfaceNormal->Set(ray->GetR(), surf->GetN());
173 
        surfaceNormal->DrawS(cx, 1.0);
 171 
  surfaceNormal->DrawS(cx, 1.0);
- 
 
 172 
  printf(" BLUE: surface normal vector\n");
174 
}
 173 
}
175 
 
 174 
 
176 
void ptt()
 175 
void ptt()
177 
{
 176 
{
178 
        for(double th = 0.0; th < 91.0; th += 5.0) {
 177 
  for(double th = 0.0; th < 91.0; th += 5.0) {
179 
                printf("%lf ", th);
 178 
      printf("%lf ", th);
180 
                PolTest(th);
 179 
      PolTest(th);
181 
        }
 180 
  }
182 
}
 181 
}
183 
//------------------------------------------------------------------------------------------
 182 
//------------------------------------------------------------------------------------------
184 
// Propagate single ray and show the statistics
 183 
// Propagate single ray and show the statistics
185 
void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
 184 
void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0)
186 
{
 185 
{
187 
  CDetector *detector = new CDetector(CENTER, parameters);
 186 
  CDetector *detector = new CDetector(CENTER, parameters);
188 
        Init();
 187 
  Init();
189 
        double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
 188 
  double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
190 
        PrintGuideHead();
 189 
  PrintGuideHead();
191 
        PrintGuideStat(izkoristek);
 190 
  PrintGuideStat(izkoristek);
192 
        DrawData(detector, parameters, theta, izkoristek);
 191 
  DrawData(detector, parameters, theta, izkoristek);
193 
}
 192 
}
194 
//------------------------------------------------------------------------------------------
 193 
//------------------------------------------------------------------------------------------
195 
// Propagate NN rays generated as grid and show the statistics
 194 
// Propagate NN rays generated as grid and show the statistics
196 
void LGG(int NN=10, double theta=0.0, bool coupling=false)
 195 
void LGG(int NN=10, double theta=0.0, bool coupling=false)
197 
{
 196 
{
198 
  parameters.setCoupling(coupling);
 197 
  parameters.setCoupling(coupling);
199 
 CDetector *detector = new CDetector(CENTER, parameters);
 198 
  CDetector *detector = new CDetector(CENTER, parameters);
200 
        Init();
 199 
  Init();
201 
        double izkoristek = Grid(detector, parameters, NN, theta);
 200 
  double izkoristek = Grid(detector, parameters, NN, theta);
202 
        //printf("izkoristek = %.3lf\n", izkoristek);
 201 
  //printf("izkoristek = %.3lf\n", izkoristek);
203 
        PrintGuideHead();
 202 
  PrintGuideHead();
204 
        PrintGuideStat(izkoristek);
 203 
  PrintGuideStat(izkoristek);
205 
        DrawData(detector, parameters, theta, izkoristek);
 204 
  DrawData(detector, parameters, theta, izkoristek);
206 
}
 205 
}
207 
//------------------------------------------------------------------------------------------
 206 
//------------------------------------------------------------------------------------------
208 
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
 207 
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
209 
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
 208 
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
210 
{
 209 
{
211 
  CDetector *detector = new CDetector(CENTER, parameters);
 210 
  CDetector *detector = new CDetector(CENTER, parameters);
212 
        Init();
 211 
  Init();
213 
        double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
 212 
  double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
214 
        //printf("izkoristek = %.3lf\n", izkoristek);
 213 
  //printf("izkoristek = %.3lf\n", izkoristek);
215 
        PrintGuideHead(); PrintGuideStat(izkoristek);
 214 
  PrintGuideHead(); PrintGuideStat(izkoristek);
216 
        DrawData(detector, parameters, theta, izkoristek);
 215 
  DrawData(detector, parameters, theta, izkoristek);
217 
}
 216 
}
218 
//------------------------------------------------------------------------------------------
 217 
//------------------------------------------------------------------------------------------
219 
// Propagate NN rays isotropically generated in solid angle theta and show the statistics
 218 
// Propagate NN rays isotropically generated in solid angle theta and show the statistics
220 
void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
 219 
void LGI(int NN = 1e4, double theta = 30.0, int nnrays = 30, int showr = 0)
221 
{
 220 
{
222 
  Init();
 221 
  Init();
223 
  CDetector *detector = new CDetector(CENTER, parameters);
 222 
  CDetector *detector = new CDetector(CENTER, parameters);
224 
  //CDetector detector = new CDetector();
 223 
  //CDetector detector = new CDetector();
225 
        double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
 224 
  double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
226 
        //printf("izkoristek = %.3lf\n", izkoristek);
 225 
  //printf("izkoristek = %.3lf\n", izkoristek);
227 
        PrintGuideHead();
 226 
  PrintGuideHead();
228 
        PrintGuideStat(izkoristek);
 227 
  PrintGuideStat(izkoristek);
229 
        DrawData(detector, parameters, theta, izkoristek);
 228 
  DrawData(detector, parameters, theta, izkoristek);
230 
        //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
 229 
  //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
231 
        //canvasDetector->cd();
 230 
  //canvasDetector->cd();
232 
        //detector->Draw();
 231 
  //detector->Draw();
233 
}
 232 
}
234 
 
 233 
 
235 
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
 234 
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
236 
{
 235 
{
237 
  Init();
 236 
  Init();
238 
  parameters.setCoupling(coupling);
 237 
  parameters.setCoupling(coupling);
239 
  CDetector *detector = new CDetector(CENTER, parameters);
 238 
  CDetector *detector = new CDetector(CENTER, parameters);
240 
         
 239 
 
241 
         const double theta = 18.5;
 240 
  const double theta = 18.5;
242 
         double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr);
 241 
  double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr);
243 
         
 242 
 
244 
         PrintGuideHead();
 243 
  PrintGuideHead();
245 
         PrintGuideStat(izkoristek);
 244 
  PrintGuideStat(izkoristek);
246 
         DrawData(detector, parameters, 18.5, izkoristek);
 245 
  DrawData(detector, parameters, 18.5, izkoristek);
247 
         
 246 
 
248 
         //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
 247 
  //TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
249 
         //canvasDetector->cd();
 248 
  //canvasDetector->cd();
250 
         //detector->Draw();
 249 
  //detector->Draw();
251 
}
 250 
}
252 
       
 251 
 
253 
       
 252 
 
254 
//------------------------------------------------------------------------------------------
 253 
//------------------------------------------------------------------------------------------
255 
void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
 254 
void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0)
256 
{
 255 
{
257 
int tmp3d = show_3d, tmpdata = show_data;
 256 
  int tmp3d = show_3d, tmpdata = show_data;
258 
double step[steps], acc[steps];
 257 
  double step[steps], acc[steps];
259 
 
 258 
 
260 
        show_3d = 0; show_data = 0;
 259 
  show_3d = 0; show_data = 0;
261 
       
 260 
 
262 
        //printf(" x_gap | acceptance \n");
 261 
  //printf(" x_gap | acceptance \n");
263 
        PrintGuideHead();
 262 
  PrintGuideHead();
264 
        for(int i=0; i<=steps; i++) {
 263 
  for(int i=0; i<=steps; i++) {
265 
                step[i] = min + i*(max - min)/(double)steps;
 264 
      step[i] = min + i*(max - min)/(double)steps;
266 
                parameters.setGap(step[i], 0.0, 0.0);
 265 
      parameters.setGap(step[i], 0.0, 0.0);
267 
                CDetector *detector = new CDetector(CENTER, parameters);
 266 
      CDetector *detector = new CDetector(CENTER, parameters);
268 
                Init();
 267 
      Init();
269 
                acc[i] = RandIso(detector, parameters, NN, 0, theta);
 268 
      acc[i] = RandIso(detector, parameters, NN, 0, theta);
270 
                //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
 269 
      //printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
271 
                PrintGuideStat(acc[i]);
 270 
      PrintGuideStat(acc[i]);
272 
        }
 271 
  }
273 
       
 272 
 
274 
        //char sbuff[256];
 273 
  //char sbuff[256];
275 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
 274 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
276 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
275 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
277 
                //      (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
 276 
  //    (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
278 
       
 277 
 
279 
        DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
 278 
  DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max);
280 
       
 279 
 
281 
        show_3d = tmp3d; show_data = tmpdata;
 280 
  show_3d = tmp3d; show_data = tmpdata;
282 
}
 281 
}
283 
//------------------------------------------------------------------------------------------
 282 
//------------------------------------------------------------------------------------------
284 
 
 283 
 
285 
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
 284 
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
286 
{
 285 
{
287 
  //int tmp3d = show_3d, tmpdata = show_data;
 286 
  //int tmp3d = show_3d, tmpdata = show_data;
288 
  double show_rays = 10;
 287 
  double show_rays = 10;
289 
  double step[steps], acc[steps];
 288 
  double step[steps], acc[steps];
290 
 
 289 
 
291 
        //show_3d = 0; show_data = 0;
 290 
  //show_3d = 0; show_data = 0;
292 
       
 291 
 
293 
        //printf(" theta | acceptance \n");
 292 
  //printf(" theta | acceptance \n");
294 
        PrintGuideHead();
 293 
  PrintGuideHead();
295 
        for(int i=0; i<=steps; i++) {
 294 
  for(int i=0; i<=steps; i++) {
296 
                Init();
 295 
      Init();
297 
                step[i] = min + i*(max - min)/(double)steps;
 296 
      step[i] = min + i*(max - min)/(double)steps;
298 
                CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
 297 
      CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
299 
                acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
 298 
      acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
300 
                //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
 299 
      //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
301 
                PrintGuideStat(acc[i]);
 300 
      PrintGuideStat(acc[i]);
302 
                delete detector;
 301 
      delete detector;
303 
        }
 302 
  }
304 
       
 303 
 
305 
        //char sbuff[256];
 304 
  //char sbuff[256];
306 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
 305 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
307 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
306 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
308 
                //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
 307 
  //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
309 
       
 308 
 
310 
        DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
 309 
  DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
311 
       
 310 
 
312 
        //show_3d = tmp3d; show_data = tmpdata;
 311 
  //show_3d = tmp3d; show_data = tmpdata;
313 
}
 312 
}
314 
 
 313 
 
315 
 
 314 
 
316 
//------------------------------------------------------------------------------------------
 315 
//------------------------------------------------------------------------------------------
317 
//void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
 316 
//void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
318 
void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
 317 
void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
319 
{
 318 
{
320 
  //int tmp3d = show_3d, tmpdata = show_data;
 319 
  //int tmp3d = show_3d, tmpdata = show_data;
321 
  double step[steps], acc[steps];
 320 
  double step[steps], acc[steps];
322 
 
 321 
 
323 
        //show_3d = 0; show_data = 0;
 322 
  //show_3d = 0; show_data = 0;
324 
       
 323 
 
325 
        //printf(" theta | acceptance \n");
 324 
  //printf(" theta | acceptance \n");
326 
        PrintGuideHead();
 325 
  PrintGuideHead();
327 
        for(int i=0; i<=steps; i++) {
 326 
  for(int i=0; i<=steps; i++) {
328 
                Init();
 327 
      Init();
329 
                step[i] = min + i*(max - min)/(double)steps;
 328 
      step[i] = min + i*(max - min)/(double)steps;
330 
                CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
 329 
      CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
331 
                //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
 330 
      //acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
332 
                acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
 331 
      acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
333 
                //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
 332 
      //printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
334 
                PrintGuideStat(acc[i]);
 333 
      PrintGuideStat(acc[i]);
335 
                delete detector;
 334 
      delete detector;
336 
        }
 335 
  }
337 
       
 336 
 
338 
        //char sbuff[256];
 337 
  //char sbuff[256];
339 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
 338 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
340 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
339 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
341 
                //      (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
 340 
  //    (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
342 
       
 341 
 
343 
        DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
 342 
  DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
344 
       
 343 
 
345 
        //show_3d = tmp3d; show_data = tmpdata;
 344 
  //show_3d = tmp3d; show_data = tmpdata;
346 
}
 345 
}
347 
//------------------------------------------------------------------------------------------
 346 
//------------------------------------------------------------------------------------------
348 
 
 347 
 
349 
void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
 348 
void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
350 
{
 349 
{
351 
  int tmp3d = show_3d, tmpdata = show_data;
 350 
  int tmp3d = show_3d, tmpdata = show_data;
352 
  show_3d = 1; show_data = 0;
 351 
  show_3d = 1; show_data = 0;
353 
  //PrintGuideHead();
 352 
  //PrintGuideHead();
354 
  TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
 353 
  TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
355 
  double min = 0.0;
 354 
  double min = 0.0;
356 
  printf("\nWait, this takes a while ");
 355 
  printf("\nWait, this takes a while ");
357 
        for(int i=0; i<=steps; i++) {
 356 
  for(int i=0; i<=steps; i++) {
358 
          double theta = min + i*(maxTheta - min) / (double)steps;
 357 
      double theta = min + i*(maxTheta - min) / (double)steps;
359 
          for (int j=0; j<=steps; j++) {
 358 
      for (int j=0; j<=steps; j++) {
360 
            Init();
 359 
          Init();
361 
                  double phi = min + j*(maxPhi - min) / (double)steps;
 360 
          double phi = min + j*(maxPhi - min) / (double)steps;
362 
                  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
 361 
          CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
363 
      double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
 362 
          double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
364 
      //PrintGuideStat(acc);
 363 
          //PrintGuideStat(acc);
365 
      hAcceptance->Fill(i, j, acc);
 364 
          hAcceptance->Fill(i, j, acc);
366 
      //printf("Acc: %f ", acc);
 365 
          //printf("Acc: %f ", acc);
367 
      delete detector;
 366 
          delete detector;
368 
      }
 367 
      }
369 
    printf(".");
 368 
      printf(".");
370 
  }
 369 
  }
371 
  printf("\n");
 370 
  printf("\n");
372 
  //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
 371 
  //DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
373 
  show_3d = tmp3d; show_data = tmpdata;
 372 
  show_3d = tmp3d; show_data = tmpdata;
374 
  TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
 373 
  TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
375 
  cp->cd();
 374 
  cp->cd();
376 
  hAcceptance->Draw("COLZ");
 375 
  hAcceptance->Draw("COLZ");
377 
}
 376 
}
378 
 
 377 
 
379 
// a vs. d
 378 
// a vs. d
380 
void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0)
 379 
void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0)
381 
{
 380 
{
382 
        //char sbuff[256];
 381 
  //char sbuff[256];
383 
 
 382 
 
384 
        //show_3d = 0; // don't show simulations
383 
  //show_3d = 0; // don't show simulations
385 
        //show_data = 1;
 384 
  //show_data = 1;
386 
       
 385 
 
387 
        //double d  = detector->GetD();
 386 
  //double d  = detector->GetD();
388 
        const double b = parameters.getB(); // upper side of LG
 387 
  const double b = parameters.getB(); // upper side of LG
389 
        //const double SiPM = 3.0; // the length of the detector itself
 388 
  //const double SiPM = 3.0; // the length of the detector itself
390 
        //double reflectivity = detector->guide->getR();
 389 
  //double reflectivity = detector->guide->getR();
391 
       
 390 
 
392 
        TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
 391 
  TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);
393 
       
 392 
 
394 
        // Use the Fresnel eq. instead of fixed reflectivity 96%
 393 
  // Use the Fresnel eq. instead of fixed reflectivity 96%
395 
        //detector->SetFresnel(1);
 394 
  //detector->SetFresnel(1);
396 
       
 395 
 
397 
        //printf("   d |   a  | Acceptance\n");
 396 
  //printf("   d |   a  | Acceptance\n");
398 
        getParameters();
 397 
  getParameters();
399 
        printf("Wait, this takes a while ");
 398 
  printf("Wait, this takes a while ");
400 
        for(int i=0; i<steps; i++) {
 399 
  for(int i=0; i<steps; i++) {
401 
                const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
 400 
      const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
402 
                for(int j=0; j<steps; j++) {
 401 
      for(int j=0; j<steps; j++) {
403 
                  Init();
 402 
          Init();
404 
                  const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
 403 
          const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
405 
                  //const double M = b/a;
 404 
          //const double M = b/a;
406 
      parameters.setGuide(a, b, d);
 405 
          parameters.setGuide(a, b, d);
407 
                  CDetector *detector = new CDetector(CENTER, parameters);
 406 
          CDetector *detector = new CDetector(CENTER, parameters);
408 
                  //detector->guide->setLG(y, M, x);
 407 
          //detector->guide->setLG(y, M, x);
409 
                  //Init(); exclude simulation
 408 
          //Init(); exclude simulation
410 
                  double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
 409 
          double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
411 
                  //double acceptance = Grid(NN, theta);
 410 
          //double acceptance = Grid(NN, theta);
412 
                  //double acceptance = -1.0;
 411 
          //double acceptance = -1.0;
413 
                  hAcceptance->Fill(d, a, acceptance);
 412 
          hAcceptance->Fill(d, a, acceptance);
414 
                  //printf("%.2lf | %.2lf | ", d, a);
 413 
          //printf("%.2lf | %.2lf | ", d, a);
415 
                  //PrintGuideStat(acceptance);
 414 
          //PrintGuideStat(acceptance);
416 
                  delete detector; //works fine, 50x50 grid takes ~4MB of RAM
 415 
          delete detector; //works fine, 50x50 grid takes ~4MB of RAM
417 
                  }
 416 
      }
418 
                printf(".");
 417 
      printf(".");
419 
        }
 418 
  }
420 
        printf("\n");
 419 
  printf("\n");
421 
       
 420 
 
422 
 
 421 
 
423 
       
 422 
 
424 
        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
 423 
  TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
425 
        cp->cd();
 424 
  cp->cd();
426 
        //TVirtualPad *pacc = cp->cd(0);
 425 
  //TVirtualPad *pacc = cp->cd(0);
427 
       
 426 
 
428 
        /*
 427 
  /*
429 
        pacc->SetRightMargin(0.10);
 428 
      pacc->SetRightMargin(0.10);
430 
        pacc->SetLeftMargin(0.10);
 429 
      pacc->SetLeftMargin(0.10);
431 
        pacc->SetTopMargin(0.10);
 430 
      pacc->SetTopMargin(0.10);
432 
        pacc->SetBottomMargin(0.10);
 431 
      pacc->SetBottomMargin(0.10);
433 
        */
 432 
   */
434 
       
 433 
 
435 
        TFile *file = new TFile("acceptance.root","RECREATE");
 434 
  TFile *file = new TFile("acceptance.root","RECREATE");
436 
        hAcceptance->Write();
 435 
  hAcceptance->Write();
437 
        file->Close();
 436 
  file->Close();
438 
        //delete file;
 437 
  //delete file;
439 
       
 438 
 
440 
        //hAcceptance->SetContour(100);
 439 
  //hAcceptance->SetContour(100);
441 
        //gStyle->SetPalette(1,0);
 440 
  //gStyle->SetPalette(1,0);
442 
        hAcceptance->SetTitle(";d [mm];a [mm]");
 441 
  hAcceptance->SetTitle(";d [mm];a [mm]");
443 
        hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
 442 
  hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
444 
        hAcceptance->Draw("COLZ");
 443 
  hAcceptance->Draw("COLZ");
445 
        char filename[128];
 444 
  char filename[128];
446 
        sprintf(filename,"LGI_ad%d.C", steps);
 445 
  sprintf(filename,"LGI_ad%d.C", steps);
447 
        cp->SaveAs(filename);
 446 
  cp->SaveAs(filename);
448 
       
 447 
 
449 
}
 448 
}
450 
 
 449 
 
451 
//------------------------------------------------------------------------------------------
 450 
//------------------------------------------------------------------------------------------
452 
// collection efficiency vs length
 451 
// collection efficiency vs length
453 
// a changes, b rests the same, d changes
 452 
// a changes, b rests the same, d changes
Line 457... Line 456...
457 
// d vs a vs acceptance
 456 
// d vs a vs acceptance
458 
void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
 457 
void LGI_d2(int NN = 1e4, double min = 0.5, double max = 3.0, const int steps = 25, double theta = 30.0)
459 
{
 458 
{
460 
  TVector3 gapGuideSiPM(0.3, 0, 0);
 459 
  TVector3 gapGuideSiPM(0.3, 0, 0);
461 
  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
 460 
  CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
462 
        int tmp3d = show_3d, tmpdata = show_data;
 461 
  int tmp3d = show_3d, tmpdata = show_data;
463 
        double step[steps], acc[steps];
 462 
  double step[steps], acc[steps];
464 
        double sipm_d, M0;
 463 
  double sipm_d, M0;
465 
        char sbuff[256];
 464 
  char sbuff[256];
- 
 
 465 
 
- 
 
 466 
  show_3d = 1; show_data = 1;
- 
 
 467 
 
- 
 
 468 
  M0 = parameters.getM();
- 
 
 469 
 
- 
 
 470 
  PrintGuideHead();
- 
 
 471 
  for(int i=0; i<=steps; i++) {
- 
 
 472 
      step[i] = min + i*(max - min)/(double)steps;
- 
 
 473 
      sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
- 
 
 474 
      parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
- 
 
 475 
      //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
- 
 
 476 
      printf("%.1lf | %.2lf | ", step[i], sipm_d);
- 
 
 477 
      Init();
- 
 
 478 
      acc[i] = RandIso(detector, parameters, NN, 0, theta);
- 
 
 479 
      PrintGuideStat(acc[i]);
- 
 
 480 
      sprintf(sbuff, "d%2d.gif", i);
- 
 
 481 
      //c3dview->SaveAs(sbuff);
- 
 
 482 
  }
- 
 
 483 
 
- 
 
 484 
 
- 
 
 485 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
- 
 
 486 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
- 
 
 487 
  //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
- 
 
 488 
 
- 
 
 489 
  //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
- 
 
 490 
 
- 
 
 491 
  TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
- 
 
 492 
  TVirtualPad *pacc = cp->cd(0);
- 
 
 493 
 
- 
 
 494 
  pacc->SetRightMargin(0.05);
- 
 
 495 
  pacc->SetLeftMargin(0.13);
- 
 
 496 
  pacc->SetTopMargin(0.08);
- 
 
 497 
 
- 
 
 498 
  TGraph *gacceptance = new TGraph(steps+1, step, acc);
- 
 
 499 
  gacceptance->SetTitle("; d [mm];Collection efficiency");
- 
 
 500 
  gacceptance->SetMarkerStyle(8);
- 
 
 501 
  gacceptance->SetLineColor(12);
- 
 
 502 
  gacceptance->SetLineWidth(1);
- 
 
 503 
  (gacceptance->GetXaxis())->SetRangeUser(min, max);
- 
 
 504 
  //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
- 
 
 505 
 
- 
 
 506 
  gacceptance->Draw("ACP");
- 
 
 507 
 
466 
 
 508 
 
467 
        show_3d = 1; show_data = 1;
 - 
 
468 
       
 - 
 
469 
        M0 = parameters.getM();
 - 
 
470 
       
 - 
 
471 
        PrintGuideHead();
 - 
 
472 
        for(int i=0; i<=steps; i++) {
 - 
 
473 
                step[i] = min + i*(max - min)/(double)steps;           
 - 
 
474 
                sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
 - 
 
475 
                parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
 - 
 
476 
                //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
 - 
 
477 
                printf("%.1lf | %.2lf | ", step[i], sipm_d);
 - 
 
478 
                Init();
 - 
 
479 
                acc[i] = RandIso(detector, parameters, NN, 0, theta);
 - 
 
480 
                PrintGuideStat(acc[i]);
 - 
 
481 
                sprintf(sbuff, "d%2d.gif", i);
 - 
 
482 
                //c3dview->SaveAs(sbuff);
 - 
 
483 
        }
 - 
 
484 
       
 - 
 
485 
       
 - 
 
486 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
 - 
 
487 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
- 
 
488 
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
 - 
 
489 
       
 - 
 
490 
        //DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max);
 - 
 
491 
       
 - 
 
492 
        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
 - 
 
493 
        TVirtualPad *pacc = cp->cd(0);
 - 
 
494 
       
 - 
 
495 
        pacc->SetRightMargin(0.05);
 - 
 
496 
        pacc->SetLeftMargin(0.13);
 - 
 
497 
        pacc->SetTopMargin(0.08);
 - 
 
498 
       
 - 
 
499 
        TGraph *gacceptance = new TGraph(steps+1, step, acc);
 - 
 
500 
        gacceptance->SetTitle("; d [mm];Collection efficiency");
 - 
 
501 
        gacceptance->SetMarkerStyle(8);
 - 
 
502 
        gacceptance->SetLineColor(12);
 - 
 
503 
        gacceptance->SetLineWidth(1);
 - 
 
504 
        (gacceptance->GetXaxis())->SetRangeUser(min, max);
 - 
 
505 
        //(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
 - 
 
506 
       
 - 
 
507 
        gacceptance->Draw("ACP");
 - 
 
508 
       
 - 
 
509 
       
 - 
 
510 
        show_3d = tmp3d; show_data = tmpdata;
 509 
  show_3d = tmp3d; show_data = tmpdata;
511 
}
 510 
}
512 
//------------------------------------------------------------------------------------------
 511 
//------------------------------------------------------------------------------------------
513 
// Acceptance vs. the length
 512 
// Acceptance vs. the length
514 
// The magnification ratio M rests the same all the time
 513 
// The magnification ratio M rests the same all the time
515 
void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
 514 
void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0)
516 
{
 515 
{
517 
        int tmp3d = show_3d, tmpdata = show_data;
 516 
  int tmp3d = show_3d, tmpdata = show_data;
518 
        double step[steps], acc[steps];
 517 
  double step[steps], acc[steps];
519 
        const double a = parameters.getA();
 518 
  const double a = parameters.getA();
520 
        const double magnif = parameters.getM();
 519 
  const double magnif = parameters.getM();
521 
 
 520 
 
522 
        //show_3d = show_in_steps; show_data = 1;
 521 
  //show_3d = show_in_steps; show_data = 1;
523 
       
 522 
 
524 
        // Use Fresnel equations
 523 
  // Use Fresnel equations
525 
        //detector->SetFresnel(1);
 524 
  //detector->SetFresnel(1);
526 
       
 525 
 
527 
        // Set glass (n=1.5) at the exit
 526 
  // Set glass (n=1.5) at the exit
528 
        //detector->SetGlass(1,0);
 527 
  //detector->SetGlass(1,0);
529 
       
 528 
 
530 
        PrintGuideHead();
 529 
  PrintGuideHead();
531 
        for(int i=0; i<=steps; i++) {
 530 
  for(int i=0; i<=steps; i++) {
532 
                step[i] = min + i*(max - min)/(double)steps;
 531 
      step[i] = min + i*(max - min)/(double)steps;
533 
                //parameters.setGuide(a, magnif, step[i]);
 532 
      //parameters.setGuide(a, magnif, step[i]);
534 
                parameters.setGuide(a, magnif, step[i]);
 533 
      parameters.setGuide(a, magnif, step[i]);
535 
                CDetector *detector = new CDetector(CENTER, parameters);
 534 
      CDetector *detector = new CDetector(CENTER, parameters);
536 
                //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
 535 
      //printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
537 
                printf("%.1lf | ", step[i]);
 536 
      printf("%.1lf | ", step[i]);
538 
                Init();
 537 
      Init();
539 
                acc[i] = RandIso(detector, parameters, NN, 0, theta);
 538 
      acc[i] = RandIso(detector, parameters, NN, 0, theta);
540 
                PrintGuideStat(acc[i]);
 539 
      PrintGuideStat(acc[i]);
541 
                delete detector;
 540 
      delete detector;
542 
        }
 541 
  }
543 
       
 542 
 
544 
        //char sbuff[256];
 543 
  //char sbuff[256];
545 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
 544 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
546 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
545 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
547 
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
 546 
  //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
548 
       
 547 
 
549 
        DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
 548 
  DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max);
550 
       
 549 
 
551 
        show_3d = tmp3d; show_data = tmpdata;
 550 
  show_3d = tmp3d; show_data = tmpdata;
552 
}
 551 
}
553 
//------------------------------------------------------------------------------------------
 552 
//------------------------------------------------------------------------------------------
554 
// Magnification optimization. a and d are fixed, M and b change in steps
 553 
// Magnification optimization. a and d are fixed, M and b change in steps
555 
void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
 554 
void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0)
556 
{
 555 
{
557 
        int tmp3d = show_3d, tmpdata = show_data;
 556 
  int tmp3d = show_3d, tmpdata = show_data;
558 
        double step[steps], acc[steps];
 557 
  double step[steps], acc[steps];
559 
       
 558 
 
560 
        const double a = parameters.getA();
 559 
  const double a = parameters.getA();
561 
        const double d = parameters.getD();
 560 
  const double d = parameters.getD();
562 
 
 561 
 
563 
        show_3d = 0; show_data = 0;
 562 
  show_3d = 0; show_data = 0;
564 
       
 563 
 
565 
        PrintGuideHead();
 564 
  PrintGuideHead();
566 
        for(int i=0; i<=steps; i++) {
 565 
  for(int i=0; i<=steps; i++) {
567 
                step[i] = min + i*(max - min)/(double)steps;
 566 
      step[i] = min + i*(max - min)/(double)steps;
568 
                //parameters.setGuide(a, step[i], d);
 567 
      //parameters.setGuide(a, step[i], d);
569 
                parameters.setGuide(a, a*step[i], d);
 568 
      parameters.setGuide(a, a*step[i], d);
570 
                CDetector *detector = new CDetector(CENTER, parameters);
 569 
      CDetector *detector = new CDetector(CENTER, parameters);
571 
                Init();
 570 
      Init();
572 
                acc[i] = RandIso(detector, parameters, NN, 0, theta);
 571 
      acc[i] = RandIso(detector, parameters, NN, 0, theta);
573 
                PrintGuideStat(acc[i]);
 572 
      PrintGuideStat(acc[i]);
574 
        }
 573 
  }
575 
       
 574 
 
576 
        //char sbuff[256];
 575 
  //char sbuff[256];
577 
        //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
 576 
  //sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf",
578 
                //      detector->GetSiPM(), detector->GetM(), detector->GetD(),
577 
  //    detector->GetSiPM(), detector->GetM(), detector->GetD(),
579 
                        //                                (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
 578 
  //                              (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta);
580 
       
 579 
 
581 
        DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
 580 
  DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max);
582 
       
 581 
 
583 
        show_3d = tmp3d; show_data = tmpdata;
 582 
  show_3d = tmp3d; show_data = tmpdata;
584 
}
 583 
}
585 
 
 584 
 
586 
//------------------------------------------------------------------------------------------
 585 
//------------------------------------------------------------------------------------------
587 
 
 586 
 
588 
void getParameters()
 587 
void getParameters()
589 
{
 588 
{
590 
  printf("LIGHT GUIDE\n"
 589 
  printf("LIGHT GUIDE\n"
591 
    "   b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
 590 
      "   b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
592 
  printf("MATERIAL REFRACITVE INDICES\n"
 591 
  printf("MATERIAL REFRACITVE INDICES\n"
593 
    "   n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
 592 
      "   n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
594 
  printf("PLATE\n"
 593 
  printf("PLATE\n"
595 
    "   ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
 594 
      "   ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
596 
  return;
 595 
  return;
597 
}
 596 
}
598 
 
 597