FairRoot/PandaRoot
PndRichGeo.cxx
Go to the documentation of this file.
1 #include "PndRichGeo.h"
2 #include "FairGeoNode.h"
3 #include "TRandom.h"
4 
6 
7 // ----- Default constructor -------------------------------------------
9  : FairGeoSet()
10 {
11  // Constructor
12  // fName has to be the name used in the geometry for all volumes.
13  // If there is a mismatch the geometry cannot be build.
14  fName="rich";
15  maxSectors=0;
16  maxModules=10;
17 
18  fRichOffset = TVector3( 0, 0, 700-35-10/*+50*/ ); //-70
19  fAlBoxSize = TVector3( 600, 600, 100 );
20  fAlBoxWallThickness = 0.05;
21  fAerogelSize = TVector3( 0*590+1*290, 0*590+1*120, 4 );
22  fAerogelOffset = TVector3( 0, 0, 1 );
23  fnOpt = std::vector<Double_t>(1,1.05);
24  fAerogelLayers = std::vector<Double_t>(1,1);
25  fAngleExtansionInner = 1.0;
26  fAngleExtansionOuter = 1.0;
27  fMirrorCurvature = 20;
28  fAngleOfMirrorPosition = 50;
29  fMirrorThickness = 0.1;
30  fMirrorLength = 290;
31  fPhDetLength = 290;
32  fPhDetThickness = 3;
33  fBeamPipeHoleX = 10;
34  fBeamPipeHoleY = 10;
35 
36  fMirrorType = 0;
37  fFlatMirrorZ = std::vector<Double_t>(2);
38  fFlatMirrorY = std::vector<Double_t>(2);
39  fFlatMirrorZGlob = std::vector<Double_t>(2);
40  fFlatMirrorYGlob = std::vector<Double_t>(2);
41  fPhDetZ = std::vector<Double_t>(2);
42  fPhDetY = std::vector<Double_t>(2);
43 
44  fSenseLevel = 4;
45  fSensorsPerDevice = 1<<(fSenseLevel-1); // number of sensores per device in one direction
46  // Quantum efficiency
47  fPhDetDev = 0; // pde_dpc3200_22.dat
48  fPhDetDev = 1; // pde_h12700.dat
49  std::string workdir(getenv( "VMCWORKDIR" ));
50  std::string effFileName = workdir;
51  Double_t keff;
52  if (fPhDetDev==0)
53  {
54  effFileName += "/detectors/rich/pde_dpc3200_22.dat";
55  fPhDetSizeX = 3.26;
56  fPhDetSizeY = 3.26;
57  fPhDetGapX = 0.01;
58  fPhDetGapY = 0.01;
59  keff = 1.7; // measured difference MC-EXP
60  }
61  if (fPhDetDev==1)
62  {
63  effFileName += "/detectors/rich/pde_h12700.dat";
64  fPhDetSizeX = 5.2;
65  fPhDetSizeY = 5.2;
66  fPhDetGapX = 0.01;
67  fPhDetGapY = 0.01;
68  keff = 1;
69  }
70  std::ifstream from( effFileName.c_str() );
71  Double_t wli, pdei;
72  from >> wli >> pdei;
73  while( !from.eof() ) {
74  fWlPhoton.push_back(wli); // nm
75  fPDE.push_back(pdei/100.0/keff); // %/100
76  from >> wli >> pdei;
77  };
78  fPhDetEff = new TGraph(fWlPhoton.size(),fWlPhoton.data(),fPDE.data());
79 
80 }
81 
82 void PndRichGeo::init(size_t ver0) {
83  size_t nRefInd = (ver0%100000)/1000;
84  size_t nlayers = (ver0%1000)/100;
85  size_t ver = ver0%100;
86  nlayers = nlayers ? nlayers : 3;
87 
88  double ka1 = angleExtansionInner();
89  double ka2 = angleExtansionOuter();
90  double beta = mirrorCurvature();
91  double alpha = angleOfMirrorPosition();
92  // Mirror
93  double za = richOffset().Z() /*+ alBoxWallThickness()*/ + aerogelOffset().Z();
94  double ya = aerogelSize().Y()/2;
95  double wa = aerogelSize().Z();
96  double yp1 = ya + wa;
97  double zp1 = wa;
98  alpha = alpha < 45 - beta ? alpha : 45 -beta;
99  alpha *= M_PI/180;
100  beta *= M_PI/180;
101  double thetaCh = acos(1/nOpt()[0]);
102  double theta = atan(ya/za);
103  double alpha1 = 2*(alpha+beta) - ka1*thetaCh;
104  double alpha2 = 2*(alpha-beta) + ka2*thetaCh + theta;
105  double alpham = (alpha1+alpha2)/2;
106  double zm1 = zp1 + yp1/tan(alpha1);
107  double ym1 = 0;
108  double zm2 = (zm1+ya*tan(alpha))/(1-tan(theta+ka2*thetaCh)*tan(alpha));
109  double ym2 = (zm2-zm1)/tan(alpha);
110  double zp2 = ((ym2-yp1+zm2*tan(alpha2))*tan(alpham)+zp1)/(1+tan(alpha2)*tan(alpham));
111  double yp2 = yp1 + (zp2-zp1)/tan(alpham);
112  if (yp2<ym2) {
113  yp2 = ym2;
114  zp2 = zm2;
115  alpham = atan((zp2-zp1)/(yp2-yp1));
116  }
117  //double phDetWidth = sqrt((zp1-zp2)*(zp1-zp2)+(yp1-yp2)*(yp1-yp2)); //[R.K. 01/2017] unused variable?
118  double zmc = (zm1+zm2)/2;
119  double ymc = (ym1+ym2)/2;
120  double wm = sqrt((zm1-zm2)*(zm1-zm2)+(ym1-ym2)*(ym1-ym2));
121  double rm = wm/2/sin(beta);
122  double hm = rm*cos(beta);
123  //double zm0 = zmc-hm*cos(alpha)-alBoxSize().Z()/2+aerogelOffset().Z(); //[R.K. 01/2017] unused variable?
124  //double ym0 = ymc+hm*sin(alpha); //[R.K. 01/2017] unused variable?
125  double theta1 = 360-(alpha+beta)*180/M_PI;
126  double theta2 = 360-(alpha-beta)*180/M_PI;
127  //double theta3 = 360-theta2; //[R.K. 01/2017] unused variable?
128 
129  // photosensor pixel sizes
130  // http://www.digitalphotoncounting.com/wp-content/uploads/PDPC_leaflet_A4_2015_10.pdf
131  fdX = 0.4075/2; //0.38016;
132  fdY = 0.4075/2; //0.32;
133  fdZ = 0;
134  fiXmax = 2*(int)(fMirrorLength/2/fdX);
135  fiYmax = 0; // see further
136 
138 
139  fMirrorRadius = rm;
140  fMirrorAxis = TVector3( 0, ymc+hm*sin(alpha), za + zmc-hm*cos(alpha) );
141  theta1 -= theta1<180 ? 0 : 360;
142  theta2 -= theta2<180 ? 0 : 360;
143  fMirrorThetaMin = theta1*M_PI/180;
144  fMirrorThetaMax = theta2*M_PI/180;
145 
146  UInt_t nm;
147  Double_t yShift, zShift, dn;
148 
149  switch (nlayers)
150  {
151  case 1:
152  fnOpt.resize(1);
153  fnOpt[0] = 1.05;
154  fAerogelLayers.resize(1);
155  fAerogelLayers[0] = 1;
156  break;
157  case 2:
158  fnOpt.resize(2);
159  dn = nRefInd ? (nRefInd-1)*0.00005 : 0.000817682;
160  fnOpt[0] = 1.05 - dn;
161  fnOpt[1] = 1.05 + dn;
162  fAerogelLayers.resize(2);
163  fAerogelLayers[0] = 0.5;
164  fAerogelLayers[1] = 0.5;
165  break;
166  case 3:
167  fnOpt.resize(3);
168  dn = nRefInd ? (nRefInd-1)*0.00005 : 0.00108652;
169  fnOpt[0] = 1.05 - dn;
170  fnOpt[1] = 1.05;
171  fnOpt[2] = 1.05 + dn;
172  fAerogelLayers.resize(3);
173  fAerogelLayers[0] = 0.333333;
174  fAerogelLayers[1] = 0.333334;
175  fAerogelLayers[2] = 0.333333;
176  break;
177  case 4:
178  fnOpt.resize(4);
179  dn = nRefInd ? (nRefInd-1)*0.00005 : 0.00122976;
180  fnOpt[0] = 1.05 - dn;
181  fnOpt[1] = 1.05;
182  fnOpt[2] = 1.05;
183  fnOpt[3] = 1.05 + dn;
184  fAerogelLayers.resize(4);
185  fAerogelLayers[0] = 0.25;
186  fAerogelLayers[1] = 0.25;
187  fAerogelLayers[2] = 0.25;
188  fAerogelLayers[3] = 0.25;
189  break;
190  case 5:
191  fnOpt.resize(5);
192  fnOpt[0] = 1.05;
193  fnOpt[1] = 1.05;
194  fnOpt[2] = 1.05;
195  fnOpt[3] = 1.05;
196  fnOpt[4] = 1.05;
197  fAerogelLayers.resize(5);
198  fAerogelLayers[0] = 0.2;
199  fAerogelLayers[1] = 0.2;
200  fAerogelLayers[2] = 0.2;
201  fAerogelLayers[3] = 0.2;
202  fAerogelLayers[4] = 0.2;
203  break;
204  default:
205  break;
206  }
207  yShift = 0;
208  zShift = richOffset().Z() + aerogelOffset().Z();
209  switch (ver)
210  {
211  case 1: // round mirror
212  fPhDetY[0] = 60;
213  fPhDetY[1] = 82.5719;
214  fPhDetZ[0] = 5.00000000000000;
215  fPhDetZ[1] = 58.957;
216  fMirrorRadius = 130.35;
217  fMirrorAxis = TVector3( 0, 100.366, -70.1726 ); //z relative to aerogel entrence
218  fMirrorAxisGlob = TVector3( 0, 100.366+yShift, -70.1726+zShift ); //z relative to aerogel entrence
219  fMirrorThetaMin = -0.878805;
220  fMirrorThetaMax = -0.13694;
221  break;
222  case 11: // flat mirror (one part)
223  case 21: // flat mirror (one part) + side mirrors
224  fFlatMirrorY.resize(2);
225  fFlatMirrorZ.resize(2);
226  fFlatMirrorY[0] = 0;
227  fFlatMirrorY[1] = 113.511542159065;
228  fFlatMirrorZ[0] = 47.6318222110121;
229  fFlatMirrorZ[1] = 129.893234066359;
230  fPhDetY[0] = 60;
231  fPhDetY[1] = fFlatMirrorY[1];
232  fPhDetZ[0] = 5.00000000000000;
233  fPhDetZ[1] = fFlatMirrorZ[1];
234  //
235  nm = fFlatMirrorZ.size();
236  fFlatMirrorYGlob.resize(nm);
237  fFlatMirrorZGlob.resize(nm);
238  yShift = 0;
239  zShift = richOffset().Z() + aerogelOffset().Z();
240  for(UInt_t i=0; i<nm; i++) {
241  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
242  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
243  }
244  break;
245  case 12: // flat mirror (two part)
246  case 22: // flat mirror (two part) + side mirrors
247  fFlatMirrorY.resize(3);
248  fFlatMirrorZ.resize(3);
249  fFlatMirrorY[0] = 0;
250  fFlatMirrorY[1] = 30.7192103643191;
251  fFlatMirrorY[2] = 91.7457580720718;
252  fFlatMirrorZ[0] = 21.0968997324958;
253  fFlatMirrorZ[1] = 53.4669621286154;
254  fFlatMirrorZ[2] = 79.0929943996688;
255  fPhDetY[0] = 60;
256  fPhDetY[1] = fFlatMirrorY[2];
257  fPhDetZ[0] = 5.00000000000000;
258  fPhDetZ[1] = fFlatMirrorZ[2];
259  //
260  nm = fFlatMirrorZ.size();
261  fFlatMirrorYGlob.resize(nm);
262  fFlatMirrorZGlob.resize(nm);
263  yShift = 0;
264  zShift = richOffset().Z() + aerogelOffset().Z();
265  for(UInt_t i=0; i<nm; i++) {
266  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
267  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
268  }
269  break;
270  case 13: // flat mirror (three part)
271  case 23: // flat mirror (three part) + side mirrors
272  fFlatMirrorY.resize(4);
273  fFlatMirrorZ.resize(4);
274  // with refraction on the surface
275  fFlatMirrorY[0] = 0;
276  fFlatMirrorY[1] = 15.7768998527303;
277  fFlatMirrorY[2] = 53.8249581943035;
278  fFlatMirrorY[3] = 86.6039907197112;
279  fFlatMirrorZ[0] = 19.8059481359439;
280  fFlatMirrorZ[1] = 36.7817321826854;
281  fFlatMirrorZ[2] = 63.6406349592571;
282  fFlatMirrorZ[3] = 67.0923693467737;
283  fPhDetY[0] = 60;
284  fPhDetY[1] = fFlatMirrorY[3];
285  fPhDetZ[0] = 5.00000000000000;
286  fPhDetZ[1] = fFlatMirrorZ[3];
287  //
288  nm = fFlatMirrorZ.size();
289  fFlatMirrorYGlob.resize(nm);
290  fFlatMirrorZGlob.resize(nm);
291  yShift = 0;
292  zShift = richOffset().Z() + aerogelOffset().Z();
293  for(UInt_t i=0; i<nm; i++) {
294  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
295  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
296  }
297  break;
298  case 14: // flat mirror (four part)
299  case 24: // flat mirror (four part) + side mirrors
300  fFlatMirrorY.resize(5);
301  fFlatMirrorZ.resize(5);
302  // with refraction on the surface
303  fFlatMirrorY[0] = 0;
304  fFlatMirrorY[1] = 8.03026468817330;
305  fFlatMirrorY[2] = 30.7475489960267;
306  fFlatMirrorY[3] = 59.1955006448811;
307  fFlatMirrorY[4] = 84.0239842069457;
308  fFlatMirrorZ[0] = 18.5825566438019;
309  fFlatMirrorZ[1] = 27.3977855428941;
310  fFlatMirrorZ[2] = 47.7739516802010;
311  fFlatMirrorZ[3] = 60.7121551030865;
312  fFlatMirrorZ[4] = 61.0707645809510;
313  fPhDetY[0] = 60;
314  fPhDetY[1] = fFlatMirrorY[4];
315  fPhDetZ[0] = 5.00000000000000;
316  fPhDetZ[1] = fFlatMirrorZ[4];
317  //
318  nm = fFlatMirrorZ.size();
319  fFlatMirrorYGlob.resize(nm);
320  fFlatMirrorZGlob.resize(nm);
321  yShift = 0;
322  zShift = richOffset().Z() + aerogelOffset().Z();
323  for(UInt_t i=0; i<nm; i++) {
324  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
325  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
326  }
327  break;
328  case 15: // flat mirror (five part)
329  case 25: // flat mirror (five part) + side mirrors
330  fFlatMirrorY.resize(6);
331  fFlatMirrorZ.resize(6);
332  // with refraction on the surface
333  fFlatMirrorY[0] = 0;
334  fFlatMirrorY[1] = 6.43722379520563;
335  fFlatMirrorY[2] = 24.4256033426640;
336  fFlatMirrorY[3] = 49.2742119933454;
337  fFlatMirrorY[4] = 65.6045650294544;
338  fFlatMirrorY[5] = 82.9932569195603;
339  fFlatMirrorZ[0] = 19.7751084893358;
340  fFlatMirrorZ[1] = 26.7049745668535;
341  fFlatMirrorZ[2] = 43.2466340126002;
342  fFlatMirrorZ[3] = 57.0589176253026;
343  fFlatMirrorZ[4] = 60.4234882450729;
344  fFlatMirrorZ[5] = 58.6650992017072;
345  fPhDetY[0] = 60;
346  fPhDetY[1] = fFlatMirrorY[5];
347  fPhDetZ[0] = 5.00000000000000;
348  fPhDetZ[1] = fFlatMirrorZ[5];
349  //
350  nm = fFlatMirrorZ.size();
351  fFlatMirrorYGlob.resize(nm);
352  fFlatMirrorZGlob.resize(nm);
353  yShift = 0;
354  zShift = richOffset().Z() + aerogelOffset().Z();
355  for(UInt_t i=0; i<nm; i++) {
356  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
357  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
358  }
359  break;
360  case 16: // flat mirror (six part)
361  case 26: // flat mirror (six part) + side mirrors
362  fFlatMirrorY.resize(7);
363  fFlatMirrorZ.resize(7);
364  // with refraction on the surface
365  fFlatMirrorY[0] = 0;
366  fFlatMirrorY[1] = 4.32030152236946;
367  fFlatMirrorY[2] = 16.3735214678972;
368  fFlatMirrorY[3] = 34.0234063082931;
369  fFlatMirrorY[4] = 53.4151371528727;
370  fFlatMirrorY[5] = 67.4943304893323;
371  fFlatMirrorY[6] = 82.0821597606174;
372  fFlatMirrorZ[0] = 20.0304540223155;
373  fFlatMirrorZ[1] = 24.6621130797753;
374  fFlatMirrorZ[2] = 36.3146724796001;
375  fFlatMirrorZ[3] = 48.7721721707139;
376  fFlatMirrorZ[4] = 56.7889550463438;
377  fFlatMirrorZ[5] = 58.5790965356983;
378  fFlatMirrorZ[6] = 56.5386444942197;
379  fPhDetY[0] = 60;
380  fPhDetY[1] = fFlatMirrorY[6];
381  fPhDetZ[0] = 5.00000000000000;
382  fPhDetZ[1] = fFlatMirrorZ[6];
383  //
384  nm = fFlatMirrorZ.size();
385  fFlatMirrorYGlob.resize(nm);
386  fFlatMirrorZGlob.resize(nm);
387  yShift = 0;
388  zShift = richOffset().Z() + aerogelOffset().Z();
389  for(UInt_t i=0; i<nm; i++) {
390  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
391  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
392  }
393  break;
394  case 17: // flat mirror (seven part)
395  case 27: // flat mirror (seven part) + side mirrors
396  fFlatMirrorY.resize(8);
397  fFlatMirrorZ.resize(8);
398  // with refraction on the surface
399  fFlatMirrorY[0] = 0;
400  fFlatMirrorY[1] = 2.07504791151454;
401  fFlatMirrorY[2] = 8.03938318086690;
402  fFlatMirrorY[3] = 22.6202807022538;
403  fFlatMirrorY[4] = 38.2473277661409;
404  fFlatMirrorY[5] = 55.1167285140505;
405  fFlatMirrorY[6] = 68.1475547317666;
406  fFlatMirrorY[7] = 81.4987977109529;
407  fFlatMirrorZ[0] = 19.0632831569999;
408  fFlatMirrorZ[1] = 21.3232597676276;
409  fFlatMirrorZ[2] = 27.5080822022865;
410  fFlatMirrorZ[3] = 40.5611035370563;
411  fFlatMirrorZ[4] = 49.9253177088730;
412  fFlatMirrorZ[5] = 56.0112344905128;
413  fFlatMirrorZ[6] = 57.2586668011183;
414  fFlatMirrorZ[7] = 55.1771069627915;
415  fPhDetY[0] = 60;
416  fPhDetY[1] = fFlatMirrorY[7];
417  fPhDetZ[0] = 5.00000000000000;
418  fPhDetZ[1] = fFlatMirrorZ[7];
419  //
420  nm = fFlatMirrorZ.size();
421  fFlatMirrorYGlob.resize(nm);
422  fFlatMirrorZGlob.resize(nm);
423  yShift = 0;
424  zShift = richOffset().Z() + aerogelOffset().Z();
425  for(UInt_t i=0; i<nm; i++) {
426  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
427  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
428  }
429  break;
430  case 18: // flat mirror (eight part)
431  case 28: // flat mirror (eight part) + side mirrors
432  fFlatMirrorY.resize(9);
433  fFlatMirrorZ.resize(9);
434  // with refraction on the surface
435  fFlatMirrorY[0] = 0;
436  fFlatMirrorY[1] = 1.53560904998618;
437  fFlatMirrorY[2] = 5.89212190378038;
438  fFlatMirrorY[3] = 17.5565807772290;
439  fFlatMirrorY[4] = 30.0104598999171;
440  fFlatMirrorY[5] = 43.1204560600741;
441  fFlatMirrorY[6] = 57.0242748178388;
442  fFlatMirrorY[7] = 69.0047342707571;
443  fFlatMirrorY[8] = 81.1329884527176;
444  fFlatMirrorZ[0] = 19.5409594890924;
445  fFlatMirrorZ[1] = 21.2004105360782;
446  fFlatMirrorZ[2] = 25.7419689325749;
447  fFlatMirrorZ[3] = 36.6887140347995;
448  fFlatMirrorZ[4] = 45.2863228266215;
449  fFlatMirrorZ[5] = 51.5860160769717;
450  fFlatMirrorZ[6] = 55.7435949020503;
451  fFlatMirrorZ[7] = 56.4461575595980;
452  fFlatMirrorZ[8] = 54.3233266479462;
453  fPhDetY[0] = 60;
454  fPhDetY[1] = fFlatMirrorY[8];
455  fPhDetZ[0] = 5.00000000000000;
456  fPhDetZ[1] = fFlatMirrorZ[8];
457  //
458  nm = fFlatMirrorZ.size();
459  fFlatMirrorYGlob.resize(nm);
460  fFlatMirrorZGlob.resize(nm);
461  yShift = 0;
462  zShift = richOffset().Z() + aerogelOffset().Z();
463  for(UInt_t i=0; i<nm; i++) {
464  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
465  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
466  }
467  break;
468  case 19: // flat mirror (nine part)
469  case 29: // flat mirror (nine part) + side mirrors
470  fFlatMirrorY.resize(10);
471  fFlatMirrorZ.resize(10);
472  // with refraction on the surface
473  fFlatMirrorY[0] = 0;
474  fFlatMirrorY[1] = 1.49645897228357;
475  fFlatMirrorY[2] = 5.64778607612293;
476  fFlatMirrorY[3] = 16.5322752536511;
477  fFlatMirrorY[4] = 27.8183214052547;
478  fFlatMirrorY[5] = 40.2374704149686;
479  fFlatMirrorY[6] = 52.9871945876487;
480  fFlatMirrorY[7] = 62.1756927165643;
481  fFlatMirrorY[8] = 71.5135677331208;
482  fFlatMirrorY[9] = 80.8735573790680;
483  fFlatMirrorZ[0] = 20.3538808252916;
484  fFlatMirrorZ[1] = 21.9497914345741;
485  fFlatMirrorZ[2] = 26.2257519619150;
486  fFlatMirrorZ[3] = 36.3783289992003;
487  fFlatMirrorZ[4] = 44.3190814791673;
488  fFlatMirrorZ[5] = 50.6718607175374;
489  fFlatMirrorZ[6] = 54.9437977944303;
490  fFlatMirrorZ[7] = 56.2339910167212;
491  fFlatMirrorZ[8] = 55.8374506290110;
492  fFlatMirrorZ[9] = 53.7178276378618;
493  fPhDetY[0] = 60;
494  fPhDetY[1] = fFlatMirrorY[9];
495  fPhDetZ[0] = 5.00000000000000;
496  fPhDetZ[1] = fFlatMirrorZ[9];
497  //
498  nm = fFlatMirrorZ.size();
499  fFlatMirrorYGlob.resize(nm);
500  fFlatMirrorZGlob.resize(nm);
501  yShift = 0;
502  zShift = richOffset().Z() + aerogelOffset().Z();
503  for(UInt_t i=0; i<nm; i++) {
504  fFlatMirrorYGlob[i] = fFlatMirrorY[i] + yShift;
505  fFlatMirrorZGlob[i] = fFlatMirrorZ[i] + zShift;
506  }
507  break;
508  default:
509  break;
510  }
511  fPhDetAngle = std::atan((fPhDetY[1]-fPhDetY[0])/(fPhDetZ[1]-fPhDetZ[0]));
512  fPhDetP0U = TVector3(0,fPhDetY[0] + yShift,fPhDetZ[0] + zShift);
513  fPhDetNxU = TVector3(1,0,0);
514  fPhDetNyU = TVector3(0,fPhDetY[1]-fPhDetY[0],fPhDetZ[1]-fPhDetZ[0]).Unit();
515  fPhDetNzU = fPhDetNxU.Cross(fPhDetNyU).Unit();
516  fPhDetP0D = TVector3(0,-fPhDetY[0] - yShift,fPhDetZ[0] + zShift);
517  fPhDetNxD = TVector3(1,0,0);
518  fPhDetNyD = TVector3(0,fPhDetY[1]-fPhDetY[0],-(fPhDetZ[1]-fPhDetZ[0])).Unit();
519  fPhDetNzD = fPhDetNxD.Cross(fPhDetNyD).Unit();
520  fiYmax = (int)(2*std::sqrt((fPhDetY[1]-fPhDetY[0])*(fPhDetY[1]-fPhDetY[0])+
521  (fPhDetZ[1]-fPhDetZ[0])*(fPhDetZ[1]-fPhDetZ[0]))/fdY);
522  // number of photo-devices
523  Double_t phDetWidth = std::sqrt((fPhDetY[1]-fPhDetY[0])*(fPhDetY[1]-fPhDetY[0])+
524  (fPhDetZ[1]-fPhDetZ[0])*(fPhDetZ[1]-fPhDetZ[0]));
525  fPhDetNumX = 2*(int)(fMirrorLength/2/(fPhDetSizeX + fPhDetGapX)+1);
526  fPhDetNumY = 2*(int)(phDetWidth/(fPhDetSizeY + fPhDetGapY)+1);
529 }
530 
532 {
533  return (( wl >= fWlPhoton.front() )&&( wl <= fWlPhoton.back() )) ? fPhDetEff->Eval(wl) : 0 ;
534 }
535 
537 {
538  //photodetector: cell size
539  if (pos.Y()>=0) {
540  TVector3 dP = pos - fPhDetP0U;
541  return TVector3(dP*fPhDetNxU,dP*fPhDetNyU,dP*fPhDetNzU);
542  }
543  else {
544  TVector3 dP = pos - fPhDetP0D;
545  return TVector3(dP*fPhDetNxD,dP*fPhDetNyD,dP*fPhDetNzD);
546  }
547 }
548 
550 {
551  if (pos.Y()>=0)
552  return fPhDetP0U + pos.X()*fPhDetNxU + pos.Y()*fPhDetNyU + pos.Z()*fPhDetNzU;
553  else
554  return fPhDetP0D + pos.X()*fPhDetNxD + pos.Y()*fPhDetNyD + pos.Z()*fPhDetNzD;
555 }
556 
557 TVector3 PndRichGeo::PositionDiscretization(TVector3 pos, bool cell)
558 {
559  // full variant
560  TVector3 posl = PhDetPositionLocal(pos);
561  Double_t xl = posl.X();
562  Double_t yl = posl.Y();
563  // to local coordinate system of the device
565  UInt_t Ix = xl/wx + fPhDetNumX/2;
566  Double_t xlc = wx*(Ix + 0.5 - fPhDetNumX/2);
568  UInt_t Iy = yl/wy + fPhDetNumY/2;
569  Double_t ylc = wy*(Iy + 0.5 - fPhDetNumY/2);
570  xl -= xlc;
571  yl -= ylc;
572  //
573  Double_t xc[4];
574  Double_t yc[4];
575  Double_t dxc[4];
576  Double_t dyc[4];
577  Double_t xcl3[2];
578  Double_t ycl3[2];
579  Double_t xcl2;
580  Double_t ycl2;
581  Double_t dx;
582  Double_t dy;
583  bool gep;
584  // dpc3200-22
585  // http://www.digitalphotoncounting.com/wp-content/uploads/PDPC_leaflet_A4_2015_10.pdf
586  if (fPhDetDev==0)
587  {
588  // pixel center x,y, pixel half widths wx,wy
589  xc[0] = 0.23; yc[0] = 0.195; dxc[0] = 0.16; dyc[0] = 0.19;
590  xc[1] = 0.56; yc[1] = 0.585; dxc[1] = 0.16; dyc[1] = 0.19;
591  xc[2] = 1.02; yc[2] = 0.975; dxc[2] = 0.16; dyc[2] = 0.19;
592  xc[3] = 1.35; yc[3] = 1.365; dxc[3] = 0.16; dyc[3] = 0.19;
593  // die center x,y
594  xcl3[0] = 0.395; ycl3[0] = 0.39; // {(xc[0]+xc[1])/2,(xc[2]+xc[3])/2}
595  xcl3[1] = 1.185; ycl3[1] = 1.17; // {(yc[0]+yc[1])/2,(yc[2]+yc[3])/2}
596  // quarter center x,y
597  xcl2 = 0.79; // (xc[1]+xc[2])/2
598  ycl2 = 0.78; // (yc[1]+yc[2])/2
599  //
600  dx = 0.395;
601  dy = 0.39;
602  // geometrical efficiency of pixel (cell filling)
603  gep = gRandom->Uniform()<=(cell?0.74:1);
604  }
605  // h12700
606  // https://www.hamamatsu.com/resources/pdf/etd/H12700_TPMH1348E.pdf
607  if (fPhDetDev==1)
608  {
609  // pixel center x,y, pixel half widths wx,wy
610  xc[0] = 0.3; yc[0] = 0.3; dxc[0] = 0.3; dyc[0] = 0.3;
611  xc[1] = 0.9; yc[1] = 0.9; dxc[1] = 0.3; dyc[1] = 0.3;
612  xc[2] = 1.5; yc[2] = 1.5; dxc[2] = 0.3; dyc[2] = 0.3;
613  xc[3] = 2.1125; yc[3] = 2.1125; dxc[3] = 0.3125; dyc[3] = 0.3125;
614  // die center x,y
615  xcl3[0] = 0.6; ycl3[0] = 0.6; // {(xc[0]+xc[1])/2,(xc[2]+xc[3])/2}
616  xcl3[1] = 1.8125; ycl3[1] = 1.8125; // {(yc[0]+yc[1])/2,(yc[2]+yc[3])/2}
617  // quarter center x,y
618  xcl2 = 1.2125; // ~(xc[1]+xc[2])/2
619  ycl2 = 1.2125; // ~(yc[1]+yc[2])/2
620  //
621  dx = 0.6;
622  dy = 0.6;
623  // geometrical efficiency of pixel (cell filling)
624  gep = true;
625  }
626  Int_t sx = xl>0?1:-1;
627  Int_t sy = yl>0?1:-1;
628  UInt_t ix = (sx*xl)/dx;
629  UInt_t iy = (sy*yl)/dy;
630  ix = ix>3?3:ix;
631  iy = iy>3?3:iy;
632  //pixel level = 4
633  Double_t hx = sx*xc[ix];
634  Double_t hy = sy*yc[iy];
635  if ((std::fabs(xl-hx)<dxc[ix])&&
636  (std::fabs(yl-hy)<dyc[iy])&&
637  gep)
638  {
639  //tile level = 1
640  if (fSenseLevel==1) {
641  hx = 0;
642  hy = 0;
643  fSensorIndexX = Ix;
644  fSensorIndexY = Iy;
645  }
646  //qurter level = 2
647  if (fSenseLevel==2) {
648  hx = sx*xcl2;
649  hy = sy*ycl2;
650  fSensorIndexX = Ix*fSensorsPerDevice + (sx<0?0:1);
651  fSensorIndexY = Iy*fSensorsPerDevice + (sy<0?0:1);
652  }
653  //die level = 3
654  if (fSenseLevel==3) {
655  UInt_t ixx = ix/2;
656  UInt_t iyy = iy/2;
657  hx = sx*xcl3[ixx];
658  hy = sy*ycl3[iyy];
659  fSensorIndexX = Ix*fSensorsPerDevice + (sx<0?1:2) + sx*ixx;
660  fSensorIndexY = Iy*fSensorsPerDevice + (sy<0?1:2) + sy*iyy;
661  }
662  //pixel level = 4
663  if (fSenseLevel==4) {
664  fSensorIndexX = Ix*fSensorsPerDevice + (sx<0?3:4) + sx*ix;
665  fSensorIndexY = Iy*fSensorsPerDevice + (sy<0?3:4) + sy*iy;
666  }
668  fSensorPosition = PhDetPositionGlobal(TVector3(hx+xlc,hy+ylc,posl.Z()));
669  return fSensorPosition; // photon hits a phdet
670  }
671  else
672  {
673  fSensorIndexX = -1;
674  fSensorIndexY = -1;
675  fSensorIndex = -1;
676  fSensorPosition = TVector3(0,0,0);
677  return fSensorPosition; // photon conversion out of sensitive region
678  }
679 }
680 
681 TVector3 PndRichGeo::PixelPosition(UInt_t ix, UInt_t iy)
682 {
683  UInt_t Ix = ix/fSensorsPerDevice; // device intex
684  UInt_t Iy = iy/fSensorsPerDevice; // device index
685  UInt_t ixl = ix - Ix*fSensorsPerDevice; // local index
686  UInt_t iyl = iy - Iy*fSensorsPerDevice; // local index
687  UInt_t nsh = fSensorsPerDevice/2;
688  Int_t sx = ixl<nsh?-1:1;
689  Int_t sy = iyl<nsh?-1:1;
690  UInt_t ixm = ixl<nsh?nsh-ixl-1:ixl-nsh;
691  UInt_t iym = iyl<nsh?nsh-iyl-1:iyl-nsh;
692  //
693  Double_t hx = (Ix+0.5-fPhDetNumX/2)*(fPhDetSizeX + fPhDetGapX); // div. center position (x)
694  Double_t hy = (Iy+0.5-fPhDetNumY/2)*(fPhDetSizeY + fPhDetGapY); // div. center position (y)
695  //
696  Double_t xc[4];
697  Double_t yc[4];
698  //Double_t dxc[4]; //[R.K.04/2017] unused
699  //Double_t dyc[4]; //[R.K.04/2017] unused
700  Double_t xcl3[2];
701  Double_t ycl3[2];
702  Double_t xcl2;
703  Double_t ycl2;
704  // dpc3200-22
705  // http://www.digitalphotoncounting.com/wp-content/uploads/PDPC_leaflet_A4_2015_10.pdf
706  if (fPhDetDev==0)
707  {
708  // pixel center x,y, pixel half widths wx,wy
709  xc[0] = 0.23; yc[0] = 0.195; //dxc[0] = 0.16; dyc[0] = 0.19; //[R.K.04/2017] unused
710  xc[1] = 0.56; yc[1] = 0.585; //dxc[1] = 0.16; dyc[1] = 0.19; //[R.K.04/2017] unused
711  xc[2] = 1.02; yc[2] = 0.975; //dxc[2] = 0.16; dyc[2] = 0.19; //[R.K.04/2017] unused
712  xc[3] = 1.35; yc[3] = 1.365; //dxc[3] = 0.16; dyc[3] = 0.19; //[R.K.04/2017] unused
713  // die center x,y
714  xcl3[0] = 0.395; ycl3[0] = 0.39; // {(xc[0]+xc[1])/2,(xc[2]+xc[3])/2}
715  xcl3[1] = 1.185; ycl3[1] = 1.17; // {(yc[0]+yc[1])/2,(yc[2]+yc[3])/2}
716  // quarter center x,y
717  xcl2 = 0.79; // (xc[1]+xc[2])/2
718  ycl2 = 0.78; // (yc[1]+yc[2])/2
719  }
720  // h12700
721  // https://www.hamamatsu.com/resources/pdf/etd/H12700_TPMH1348E.pdf
722  if (fPhDetDev==1)
723  {
724  // pixel center x,y, pixel half widths wx,wy
725  xc[0] = 0.3; yc[0] = 0.3; //dxc[0] = 0.3; dyc[0] = 0.3; //[R.K.04/2017] unused
726  xc[1] = 0.9; yc[1] = 0.9; //dxc[1] = 0.3; dyc[1] = 0.3; //[R.K.04/2017] unused
727  xc[2] = 1.5; yc[2] = 1.5; //dxc[2] = 0.3; dyc[2] = 0.3; //[R.K.04/2017] unused
728  xc[3] = 2.1125; yc[3] = 2.1125; //dxc[3] = 0.3125; dyc[3] = 0.3125; //[R.K.04/2017] unused
729  // die center x,y
730  xcl3[0] = 0.6; ycl3[0] = 0.6; // {(xc[0]+xc[1])/2,(xc[2]+xc[3])/2}
731  xcl3[1] = 1.8125; ycl3[1] = 1.8125; // {(yc[0]+yc[1])/2,(yc[2]+yc[3])/2}
732  // quarter center x,y
733  xcl2 = 1.2125; // ~(xc[1]+xc[2])/2
734  ycl2 = 1.2125; // ~(yc[1]+yc[2])/2
735  }
736  //tile level = 1
737  if (fSenseLevel==1) {
738  hx += 0;
739  hy += 0;
740  }
741  //qurter level = 2
742  if (fSenseLevel==2) {
743  hx += sx*xcl2;
744  hy += sy*ycl2;
745  }
746  //die level = 3
747  if (fSenseLevel==3) {
748  hx += sx*xcl3[ixm];
749  hy += sy*ycl3[iym];
750  }
751  //pixel level = 4
752  if (fSenseLevel==4) {
753  hx += sx*xc[ixm];
754  hy += sy*yc[iym];
755  }
756  fSensorIndexX = ix;
757  fSensorIndexY = iy;
759  return PhDetPositionGlobal(TVector3(hx,hy,0));
760 }
761 
763  Double_t dX,
764  Double_t dY,
765  Double_t dZ)
766 {
767  // simple variant
768  Double_t dX_ = dX>0 ? dX : fdX;
769  Double_t dY_ = dY>0 ? dY : fdY;
770  Double_t dZ_ = dZ>0 ? dZ : fdZ;
771  Double_t x = pos.X();
772  Double_t y = pos.Y();
773  Double_t z = pos.Z();
774  if (x&&dX) x = ((int)(x/dX_) + x/std::fabs(x)/2)*dX_;
775  if (y&&dY) y = ((int)(y/dY_) + y/std::fabs(y)/2)*dY_;
776  if (z&&dZ) z = ((int)(z/dZ_) + z/std::fabs(z)/2)*dZ_;
777  return TVector3(x,y,z);
778 }
779 
780 UInt_t PndRichGeo::IndexX(TVector3 pos)
781 {
782  if (pos!=fSensorPosition) PositionDiscretization(pos,false);
783  return fSensorIndexX;
784  //return (int)((pos.X()+fdX*fiXmax/2)/fdX);
785 }
786 
787 UInt_t PndRichGeo::IndexY(TVector3 pos)
788 {
789  if (pos!=fSensorPosition) PositionDiscretization(pos,false);
790  return fSensorIndexY;
791  //return (int)((pos.Y()+fdY*fiXmax/2)/fdY);
792 }
793 
794 TVector3 PndRichGeo::PixelPositionLocal(UInt_t ix, UInt_t iy)
795 {
796  Double_t x = ix*fdX-fiXmax*fdX/2+0.5*fdX;
797  Double_t y = iy*fdY-fiYmax*fdY/2+0.5*fdY;
798  return TVector3(x,y,0);
799 }
800 
801 TVector3 PndRichGeo::PixelPositionGlobal(UInt_t ix, UInt_t iy)
802 {
804 }
805 
806 // -------------------------------------------------------------------------
807 
808 const char* PndRichGeo::getModuleName(Int_t m)
809 {
815  sprintf(modName,"rich0%i",m+1);
816  return modName;
817 }
818 
819 const char* PndRichGeo::getEleName(Int_t m)
820 {
822  sprintf(eleName,"rich0%i",m+1);
823  return eleName;
824 }
Int_t fSensorIndexY
Definition: PndRichGeo.h:64
UInt_t fSensorsPerDevice
Definition: PndRichGeo.h:62
Double_t mirrorCurvature()
Definition: PndRichGeo.h:144
UInt_t fiYmax
Definition: PndRichGeo.h:83
std::vector< Double_t > fnOpt
refraction index of the aerogel
Definition: PndRichGeo.h:18
char eleName[20]
Definition: PndRichGeo.h:87
TVector3 pos
friend F32vec4 acos(const F32vec4 &a)
Definition: P4_F32vec4.h:113
UInt_t fPhDetNumX
Definition: PndRichGeo.h:73
Double_t phDetQEff(Double_t wl)
Definition: PndRichGeo.cxx:531
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
std::vector< Double_t > fFlatMirrorY
No idea (SS)
Definition: PndRichGeo.h:50
double fMirrorRadius
Definition: PndRichGeo.h:44
double dy
double fAerogelEntrancePositionZ
Definition: PndRichGeo.h:41
UInt_t IndexY(TVector3 pos)
Definition: PndRichGeo.cxx:787
TVector3 aerogelOffset()
Definition: PndRichGeo.h:129
UInt_t fiXmax
Definition: PndRichGeo.h:83
TVector3 fMirrorAxis
Definition: PndRichGeo.h:45
TVector3 PhDetPositionLocal(TVector3 pos)
Definition: PndRichGeo.cxx:536
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
TVector3 aerogelSize()
Definition: PndRichGeo.h:126
std::vector< Double_t > fPhDetZ
No idea (SS)
Definition: PndRichGeo.h:55
Double_t fPhDetGapX
Definition: PndRichGeo.h:71
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TVector3 PositionDiscretization(TVector3 pos, bool cell=true)
Definition: PndRichGeo.cxx:557
Double_t angleOfMirrorPosition()
Definition: PndRichGeo.h:147
const char * getEleName(Int_t)
Definition: PndRichGeo.cxx:819
UInt_t fPhDetPixelNumX
Definition: PndRichGeo.h:75
TGraph * fPhDetEff
Definition: PndRichGeo.h:59
UInt_t fPhDetDev
Definition: PndRichGeo.h:78
double fMirrorThetaMax
Definition: PndRichGeo.h:43
TVector3 fPhDetNyU
Definition: PndRichGeo.h:80
std::vector< Double_t > fFlatMirrorZGlob
No idea (SS)
Definition: PndRichGeo.h:51
std::vector< Double_t > fFlatMirrorYGlob
No idea (SS)
Definition: PndRichGeo.h:52
TVector3 fPhDetNxU
Definition: PndRichGeo.h:80
TVector3 fPhDetNzU
Definition: PndRichGeo.h:80
TVector3 PhDetPositionGlobal(TVector3 pos)
Definition: PndRichGeo.cxx:549
std::vector< Double_t > fFlatMirrorZ
No idea (SS)
Definition: PndRichGeo.h:49
Double_t
UInt_t IndexX(TVector3 pos)
Definition: PndRichGeo.cxx:780
std::vector< Double_t > fWlPhoton
Definition: PndRichGeo.h:57
Int_t fSensorIndex
Definition: PndRichGeo.h:65
Double_t fPhDetGapY
Definition: PndRichGeo.h:72
void init(size_t ver=0)
Definition: PndRichGeo.cxx:82
Double_t z
Double_t angleExtansionOuter()
Definition: PndRichGeo.h:141
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
std::vector< Double_t > fAerogelLayers
No idea (SS)
Definition: PndRichGeo.h:19
TVector3 PixelPositionLocal(UInt_t ix, UInt_t iy)
Definition: PndRichGeo.cxx:794
TVector3 PixelPositionGlobal(UInt_t ix, UInt_t iy)
Definition: PndRichGeo.cxx:801
Double_t fdZ
Definition: PndRichGeo.h:82
double dx
UInt_t fSenseLevel
Definition: PndRichGeo.h:61
UInt_t fPhDetNumY
Definition: PndRichGeo.h:74
Double_t fPhDetAngle
Definition: PndRichGeo.h:79
Double_t angleExtansionInner()
Definition: PndRichGeo.h:138
TVector3 LocalPositionDiscretization(TVector3 pos, Double_t dX=-1, Double_t dY=-1, Double_t dZ=-1)
Definition: PndRichGeo.cxx:762
Double_t x
Double_t fPhDetSizeX
Definition: PndRichGeo.h:69
std::vector< Double_t > fPhDetY
No idea (SS)
Definition: PndRichGeo.h:56
TVector3 richOffset()
Definition: PndRichGeo.h:117
TVector3 fPhDetNzD
Definition: PndRichGeo.h:81
TVector3 fMirrorAxisGlob
Definition: PndRichGeo.h:46
double fMirrorLength
Mirror length [cm].
Definition: PndRichGeo.h:30
ClassImp(PndAnaContFact)
TVector3 fPhDetNyD
Definition: PndRichGeo.h:81
char modName[20]
Definition: PndRichGeo.h:86
Double_t y
double alpha
Definition: f_Init.h:9
const char * getModuleName(Int_t)
Definition: PndRichGeo.cxx:808
std::vector< Double_t > nOpt()
Definition: PndRichGeo.h:132
UInt_t fPhDetPixelNumY
Definition: PndRichGeo.h:76
Double_t fdX
Definition: PndRichGeo.h:82
Double_t fdY
Definition: PndRichGeo.h:82
TVector3 fSensorPosition
Definition: PndRichGeo.h:66
TVector3 fPhDetP0D
Definition: PndRichGeo.h:81
Double_t fPhDetSizeY
Definition: PndRichGeo.h:70
Int_t fSensorIndexX
Definition: PndRichGeo.h:63
TVector3 PixelPosition(UInt_t ix, UInt_t iy)
Definition: PndRichGeo.cxx:681
TVector3 fPhDetP0U
Definition: PndRichGeo.h:80
double fMirrorThetaMin
Definition: PndRichGeo.h:42
TVector3 fPhDetNxD
Definition: PndRichGeo.h:81