FairRoot/PandaRoot
createRootGeometry_DIRC_fsEVdroplens_MCPs.C
Go to the documentation of this file.
1 // macro that creates the Barrel DIRC geometry
2 // to be able to run this macro please be sure you have the following configuration in gconfig/g4Config.C file:
3 // TG4RunConfiguration* runConfiguration
4 // = new TG4RunConfiguration("geomRoot", "QGSP_BERT_EMV+optical", "stepLimiter+specialCuts+specialControls");
5 
6 // input parameters are:
7 // fFocusingSystem = 0 - no focusing is used
8 // fFocusingSystem = 1 - lens
9 // fFocusingSystem = 2 - forward mirror is used
10 // fFocusingSystem = 3 - new lens (no air gap) is used
11 // fprizm = kFALSE - no prism
12 // fprizm = kTRUE - prism
13 
14 // allowed configurations:
15 // fFocusingSystem = 0 && fprizm = kFALSE output file is in geometry/dirc_l0_p0.root
16 // fFocusingSystem = 0 && fprizm = kTRUE output file is in geometry/dirc_l0_p1.root
17 // fFocusingSystem = 1 && fprizm = kFALSE output file is in geometry/dirc_l1_p0.root
18 // fFocusingSystem = 2 && fprizm = kFALSE output file is in geometry/dirc_l2_p0.root
19 // fFocusingSystem = 2 && fprizm = kTRUE output file is in geometry/dirc_l2_p1.root
20 // corresponding .root files with geometry please fone in the geometry directory
21 
22 // B E F O R E R U N N I N G T H E P R O G R A M
23 // C H E C K T H E N A M E O F T H E O U T P U T F I L E ! ! !
24 const Double_t pi = 3.1415926535;
25 
26 void getXrange(double dphi, double radiusCornerOut, double radiusCornerIn, double MCPactiveArea, int cas1, int cas2, double ylow, double step, double &x1, double &x2);
27 double getX(TVector3 corner1, TVector3 corner2, double yy);
28 int findSectorIn(double y, double dphi, double radius, double EVdrop, double hthick);
29 int findSectorOut(double y, double dphi_rad, double radiusCornerOut);
30 
31 int createRootGeometry_DIRC_fsEVdroplens_MCPs(Int_t fFocusingSystem = 0, Bool_t fprizm = kFALSE){
32 
33  //const Double_t pi = 3.1415926535;
34 
35  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
36 
37  TString vmcWorkdir = getenv("VMCWORKDIR");
38 
39  // Load this libraries
40  gSystem->Load("libGeoBase");
41  gSystem->Load("libParBase");
42  gSystem->Load("libBase");
43  gSystem->Load("libPndData");
44  gSystem->Load("libPassive");
45 
46  // variable to achieve the DIRC basic parameters
47  PndGeoDrc* fGeo = new PndGeoDrc();
48 
49  // units = cm
50 
51 
52  Double_t eps = 0.01; // epsilon
53  Double_t mirr_hthick = 0.02;
54  Double_t PDbaseLayer = 5.; // [cm] thickness of the carbon at the back of the EV
55 
56  Double_t radius = fGeo->radius(); // 50. radius in middle of the barbox (x and y)
57  Double_t hthick = fGeo->barHalfThick(); // 1.7/2. half thickness of the bars
58  Double_t barnum = fGeo->barNum(); // 6 number of bars per barbox
59  Double_t bbnum = fGeo->BBoxNum(); //16. total number of sides = barboxes
60  Double_t bbGap = fGeo->BBoxGap(); //1.5 gap btw the neighboring barboxes (at the middle height)
61  Double_t pipehAngle = fGeo->PipehAngle(); //3.6 [degrees] half of the angular space needed for the target pipe
62  Double_t barWin_hthick = 0.1; // [cm] thickness of the 'glas' window at the readout end of the barbox;
63 
64  Double_t bbox_zdown = fGeo->barBoxZDown(); // 130. bar box z downstream
65  Double_t bbox_zup = fGeo->barBoxZUp(); // -120. bar box z upstream
66  Double_t bbox_hlen = 0.5*(bbox_zdown - bbox_zup); // 125. bar box half length
67  Double_t bbox_shift = bbox_zup + bbox_hlen; // 5. bar box shift
68  Double_t bargap = fGeo->barGap(); // 0.01 half gap between bars
69  Double_t boxgap = fGeo->boxGap(); // 0.1 gap between bars and the bar box
70  Double_t boxthick = fGeo->boxThick(); // 0.05 thickness of the bar box
71  Double_t gluehthick = 0.01;//[cm] // glue layer, connects two halfs into one long radiator bar
72  Double_t len = 0.; // length of the lenses block. see further
73  Double_t fSlabEnd = 0.; // [cm] position at which bar in front of EV ends
74  Double_t fdz_mirr1 = 0.;
75  Double_t fdz_mirr2 = 0.;
76  Double_t fdz_lens3 = 0.;
77  Double_t fdz_lens2 = 0.;
78  Double_t fdz_lens1 = 0.;
79 
80  Double_t sob_len = fGeo->EVlen(); // 30. in current version
81  Double_t sob_shift = -bbox_hlen + bbox_shift - sob_len; // -150.
82  Double_t sob_angleB = fGeo->EVbackAngle(); //90. [degrees] angle of the EV (usually it is 90)
83  Double_t sob_angle = 45.;//fGeo->EVangle(); //80. [degrees] opening angle of the EV ############################
84  Double_t sob_Rout = radius + hthick + sob_len*tan(sob_angle/180.*pi)/cos(pi/bbnum/2.);
85 
86  Double_t bbAngle = ( 180. - 2.*pipehAngle - bbGap/radius/pi*180.*(bbnum/2.-1.) )/(bbnum/2.); // ~20 degrees
87  Double_t bbX = radius*bbAngle/180.*pi; // 17.45 cm
88  cout<<"bbAngle = "<<bbAngle<<", bbX = "<<bbX<<endl;
89 
90 
91  Double_t phi0 = (180.-2.*pipehAngle)/bbnum + pipehAngle;
92  Double_t dphi = (180.-2.*pipehAngle)/bbnum*2.;
93  Double_t dphi_rad = dphi/180.*pi;
94 
95  cout<<"dphi = "<<dphi<<", phi0 = "<<phi0<<endl;
96  //----------------------------------------------------------------------
97 
98  //parameters for MCPs:
99  Double_t MCPsize = 5.9; //[cm] width=height of an MCP
100  Double_t MCPactiveArea = 5.3; //[cm] width=height of the active area of an MCP
101  Double_t MCPgap = 0.1; //[cm] gap between MCPs
102  Double_t PixelSize = 0.65; //[cm] size of a pixel
103  Int_t Npix = 8; // number of mcps in a row/column
104  const Int_t NpixTotal = Npix*Npix; // total number of pixels for 1 mcp
105  Double_t PixelGap = (MCPactiveArea - Double_t(Npix)*PixelSize)/(Double_t(Npix)-1.);
106  cout<<"pixel gap = "<< PixelGap<<endl;
107  Double_t PDwindowThick = 0.1; //[cm]
108  Double_t PhCathodeThick = 0.1; //[cm]
109  Double_t PDsensitiveThick = 0.1; //[cm]
110  Double_t PDgreaseLayer = 0.1; //[cm]
111 
112  // parameters for the EV
113  Double_t EVdrop = fGeo->EVdrop(); //1. [cm] drop of the EV - inner radius ##############
114  Double_t EVoffset = fGeo->EVoffset(); //0. [cm] offset of the EV - outer radius
115  Double_t EVgreaseLayer = 0.1; //[cm] grease layer thickness btw the bar window and the EV
116 
117  // prizm:
118  Double_t phlength = fGeo->PrismhLength(); //4.5 [cm] half length of the prizm
119  Double_t pdrop = fGeo->PrismDrop(); //0.5 [cm] drop of the prizm - inner side
120  Double_t pangle = fGeo->PrismAngle(); //30. [degrees] angle of the edge
121  Double_t poffset = fGeo->PrismOffset(); //1.[cm] prizm offset - outer side
122  Double_t pheight = 2.*phlength * tan(pangle/180.*pi); // [cm] half heigth of the prism
123 
124  Double_t sob_Rprizm = radius + hthick + poffset + pheight + EVoffset + (sob_len-2.*phlength)*tan(60./180.*pi);
125 
126  cout<<"prizm length = "<<phlength*2.<<endl;
127  cout<<"prizm height = "<<2.*phlength * tan(pangle/180.*pi)<<endl;
128 
129 
130  Double_t radiusCornerIn = (radius-hthick-EVdrop)/cos(dphi/2./180.*pi);
131  Double_t radiusCornerOut = sob_Rout;
132  Double_t radiusMiddleSmall = radius-hthick-EVdrop;
133  //----------------------------------------------------------
134 
135  // Rotations:
136  TGeoRotation rot1;
137  rot1.RotateZ(90.);
138 
139  TString fGeoFile= Form("../../geometry/dirc_l%d_p%d_mirrorGap_MCPs1_updated.root",fFocusingSystem, fprizm);
140  TFile* fi = new TFile(fGeoFile,"RECREATE");
141  cout<<"Output file = "<<fGeoFile<<endl;
142 
143  Double_t par[10];
144  Double_t aa;
145  Double_t z, density, radl, absl, w;
146  Int_t nel, numed, nz;
147 
148  new TGeoManager("Drc", "Drc");
149 
150  // MATERIALS, MIXTURES AND TRACKING MEDIA
151 // Mixture: air
152  nel = 3;
153  density = 0.001205;
154  TGeoMixture *air = new TGeoMixture("air", nel,density);
155  aa = 14.010000; z = 7.000000; w = 0.780000; // N
156  air->DefineElement(0,aa,z,w);
157  aa = 16.000000; z = 8.000000; w = 0.210000; // O
158  air->DefineElement(1,aa,z,w);
159  aa = 39.950000; z = 18.000000; w = 0.010000; // AR
160  air->DefineElement(2,aa,z,w);
161  //air->SetIndex(0);
162 // Medium: air
163  numed = 1; // medium number
164  par[0] = 0.000000; // isvol
165  par[1] = 1.000000; // ifield
166  par[2] = 30.000000; // fieldm
167  par[3] = -1.000000; // tmaxfd
168  par[4] = -1.000000; // stemax
169  par[5] = -1.000000; // deemax
170  par[6] = 0.001000; // epsil
171  par[7] = -1.000000; // stmin
172  TGeoMedium *air_m = new TGeoMedium("air", numed,air, par);
173 
174 // Mixture: DIRCairNoSens
175  nel = 3;
176  density = 0.001205;
177  TGeoMixture *DIRCairNoSens = new TGeoMixture("DIRCairNoSens", nel,density);
178  aa = 14.010000; z = 7.000000; w = 0.780000; // N
179  DIRCairNoSens->DefineElement(0,aa,z,w);
180  aa = 16.000000; z = 8.000000; w = 0.210000; // O
181  DIRCairNoSens->DefineElement(1,aa,z,w);
182  aa = 39.950000; z = 18.000000; w = 0.010000; // AR
183  DIRCairNoSens->DefineElement(2,aa,z,w);
184  //DIRCairNoSens->SetIndex(1);
185 // Medium: DIRCairNoSens
186  numed = 2; // medium number
187  par[0] = 0.000000; // isvol
188  par[1] = 1.000000; // ifield
189  par[2] = 30.000000; // fieldm
190  par[3] = -1.000000; // tmaxfd
191  par[4] = -1.000000; // stemax
192  par[5] = -1.000000; // deemax
193  par[6] = 0.001000; // epsil
194  par[7] = -1.000000; // stmin
195  TGeoMedium *DIRCairNoSens_m = new TGeoMedium("DIRCairNoSens", numed, DIRCairNoSens, par);
196 
197  // Mixture: Epotek301_2 (optical glue)
198  nel = 2;
199  density = 2.200000;
200  TGeoMixture* Epotek301_2 = new TGeoMixture("Epotek301_2", nel,density);
201  aa = 28.090000; z = 14.000000; w = 0.467475; // SI
202  Epotek301_2->DefineElement(0,aa,z,w);
203  aa = 15.999400; z = 8.000000; w = 0.532525; // O
204  Epotek301_2->DefineElement(1,aa,z,w);
205  //Epotek301_2->SetIndex(0);
206 // Medium: Epotek
207  numed = 5; // medium number
208  par[0] = 1.000000; // isvol
209  par[1] = 1.000000; // ifield
210  par[2] = 20.000000; // fieldm
211  par[3] = -1.000000; // tmaxfd
212  par[4] = -1.000000; // stemax
213  par[5] = -1.000000; // deemax
214  par[6] = 0.001000; // epsil
215  par[7] = -1.000000; // stmin
216  TGeoMedium *Epotek301_2_m = new TGeoMedium("Epotek301_2", numed, Epotek301_2, par);
217 
218  // Mixture: OpticalGrease
219  nel = 2;
220  density = 2.200000;
221  TGeoMixture* OpticalGrease = new TGeoMixture("OpticalGrease", nel,density);
222  aa = 28.090000; z = 14.000000; w = 0.467475; // SI
223  OpticalGrease->DefineElement(0,aa,z,w);
224  aa = 15.999400; z = 8.000000; w = 0.532525; // O
225  OpticalGrease->DefineElement(1,aa,z,w);
226  //OpticalGrease->SetIndex(0);
227 // Medium: OpticalGrease
228  numed = 6; // medium number
229  par[0] = 1.000000; // isvol
230  par[1] = 1.000000; // ifield
231  par[2] = 20.000000; // fieldm
232  par[3] = -1.000000; // tmaxfd
233  par[4] = -1.000000; // stemax
234  par[5] = -1.000000; // deemax
235  par[6] = 0.001000; // epsil
236  par[7] = -1.000000; // stmin
237  TGeoMedium *OpticalGrease_m = new TGeoMedium("OpticalGrease", numed, OpticalGrease, par);
238 
239 // Material: DIRCcarbonFiber
240  aa = 12.011000;
241  z = 6.000000;
242  density = 2.500000;
243  radl = 16.999068;
244  absl = 32.080525;
245  TGeoMaterial *DIRCcarbonFiber = new TGeoMaterial("DIRCcarbonFiber", aa,z,density,radl,absl);
246  //DIRCcarbonFiber->SetIndex(6);
247 // Medium: DIRCcarbonFiber
248  numed = 3; // medium number
249  par[0] = 0.000000; // isvol
250  par[1] = 0.000000; // ifield
251  par[2] = 20.000000; // fieldm
252  par[3] = -1.000000; // tmaxfd
253  par[4] = -1.000000; // stemax
254  par[5] = -1.000000; // deemax
255  par[6] = 0.001000; // epsil
256  par[7] = -1.000000; // stmin
257  TGeoMedium *DIRCcarbonFiber_m = new TGeoMedium("DIRCcarbonFiber", numed,DIRCcarbonFiber, par);
258 
259 // Mixture: FusedSil
260  nel = 2;
261  density = 2.200000;
262  TGeoMixture* FusedSil = new TGeoMixture("FusedSil", nel,density);
263  aa = 28.090000; z = 14.000000; w = 0.467475; // SI
264  FusedSil->DefineElement(0,aa,z,w);
265  aa = 15.999400; z = 8.000000; w = 0.532525; // O
266  FusedSil->DefineElement(1,aa,z,w);
267  //FusedSil->SetIndex(0);
268 // Medium: FusedSil
269  numed = 4; // medium number
270  par[0] = 1.000000; // isvol
271  par[1] = 1.000000; // ifield
272  par[2] = 20.000000; // fieldm
273  par[3] = -1.000000; // tmaxfd
274  par[4] = -1.000000; // stemax
275  par[5] = -1.000000; // deemax
276  par[6] = 0.001000; // epsil
277  par[7] = -1.000000; // stmin
278  TGeoMedium *FusedSil_m = new TGeoMedium("FusedSil", numed,FusedSil, par);
279 
280 // Mixture: Mirror
281  nel = 2;
282  density = 2.200000;
283  TGeoMixture* Mirror = new TGeoMixture("Mirror", nel,density);
284  aa = 28.090000; z = 14.000000; w = 0.467475; // SI
285  Mirror->DefineElement(0,aa,z,w);
286  aa = 15.999400; z = 8.000000; w = 0.532525; // O
287  Mirror->DefineElement(1,aa,z,w);
288  //Mirror->SetIndex(4);
289 // Medium: Mirror
290  numed = 5; // medium number
291  par[0] = 0.000000; // isvol
292  par[1] = 0.000000; // ifield
293  par[2] = 20.000000; // fieldm
294  par[3] = -1.000000; // tmaxfd
295  par[4] = -1.000000; // stemax
296  par[5] = -1.000000; // deemax
297  par[6] = 0.001000; // epsil
298  par[7] = -1.000000; // stmin
299  TGeoMedium *Mirror_m = new TGeoMedium("Mirror", numed,Mirror, par);
300 
301 // Mixture: Marcol82
302  nel = 2;
303  density = 0.850000;
304  TGeoMixture *Marcol82 = new TGeoMixture("Marcol82", nel,density);
305  aa = 1.007940; z = 1.000000; w = 0.148605; // H
306  Marcol82->DefineElement(0,aa,z,w);
307  aa = 12.010700; z = 6.000000; w = 0.851395; // C
308  Marcol82->DefineElement(1,aa,z,w);
309  //Marcol82->SetIndex(5);
310 // Medium: Marcol82
311  numed = 6; // medium number
312  par[0] = 0.000000; // isvol
313  par[1] = 1.000000; // ifield
314  par[2] = 30.000000; // fieldm
315  par[3] = -1.000000; // tmaxfd
316  par[4] = -1.000000; // stemax
317  par[5] = -1.000000; // deemax
318  par[6] = 0.001000; // epsil
319  par[7] = -1.000000; // stmin
320  TGeoMedium *Marcol82_m = new TGeoMedium("Marcol82", numed,Marcol82, par);
321 
322 // Mixture: Marcol82-7
323  nel = 2;
324  density = 0.850000;
325  TGeoMixture *Marcol82_7 = new TGeoMixture("Marcol82-7", nel,density);
326  aa = 1.007940; z = 1.000000; w = 0.148605; // H
327  Marcol82_7->DefineElement(0,aa,z,w);
328  aa = 12.010700; z = 6.000000; w = 0.851395; // C
329  Marcol82_7->DefineElement(1,aa,z,w);
330  //Marcol82_7->SetIndex(5);
331 // Medium: Marcol82_7
332  numed = 7; // medium number
333  par[0] = 0.000000; // isvol
334  par[1] = 1.000000; // ifield
335  par[2] = 30.000000; // fieldm
336  par[3] = -1.000000; // tmaxfd
337  par[4] = -1.000000; // stemax
338  par[5] = -1.000000; // deemax
339  par[6] = 0.001000; // epsil
340  par[7] = -1.000000; // stmin
341  TGeoMedium *Marcol82_7_m = new TGeoMedium("Marcol82-7", numed,Marcol82_7, par);
342 
343 
344 // Mixture: NLAK33A
345  nel = 2;
346  density = 4.220000;
347  TGeoMixture *NLAK33A = new TGeoMixture("NLAK33A", nel, density);
348  aa = 28.090000; z = 14.000000; w = 0.467475; // SI
349  NLAK33A->DefineElement(0,aa,z,w);
350  aa = 15.999400; z = 8.000000; w = 0.532525; // O
351  NLAK33A->DefineElement(1,aa,z,w);
352  //NLAK33A->SetIndex(2);
353 // Medium: NLAK33A
354  numed = 8; // medium number
355  par[0] = 1.000000; // isvol
356  par[1] = 1.000000; // ifield
357  par[2] = 20.000000; // fieldm
358  par[3] = -1.000000; // tmaxfd
359  par[4] = -1.000000; // stemax
360  par[5] = -1.000000; // deemax
361  par[6] = 0.001000; // epsil
362  par[7] = -1.000000; // stmin
363  TGeoMedium *NLAK33A_m = new TGeoMedium("NLAK33A", numed,NLAK33A, par);
365 
366 
367 
368  // focusing systems:
369  // 0 no ficusing:
370  if(fFocusingSystem == 0){
371  Double_t flen = 0.;
372  fSlabEnd = -bbox_hlen + bbox_shift + flen;
373  cout<<"bar ends at = "<<fSlabEnd<<endl;
374 
375  //fAtBarEnd = "DrcBar";
376  }
377 // 1 lens
378  if (fFocusingSystem == 1){ // L E N S E S (no airgaps, thin nlak33)
379  // some notes to lens operations (revision 9649) is in
380  // ~carsten/work/documents/software/PandaRoot/lens_definitions_2.pdf
381 
382  Double_t r = 75.18; // first lens radius (cm)
383  Double_t alpha = TMath::ASin(bbX/12./r); // lside -> bbX, assuming that there are 6 bars in the bar box??
384  Double_t a = r - r*TMath::Cos(alpha);
385  Double_t b = a + 0.5; // box dimension .6 instead of .5 due to strong curvature
386 
387  cout<<" DIRC a,b = "<<a<<" "<<b<<endl;
388 
389 
390  Double_t r2 = 75.18; // radius second lens (cm)
391  Double_t b2 = 0.2; // + a2;
392 
393  Double_t r3 = 17.95; // third lens radius (cm)
394  Double_t alpha3 = TMath::ASin(bbX/12./r3); // lside -> bbX, assuming that there are 6 bars in the bar box??
395  Double_t a3 = r3 - r3*TMath::Cos(alpha3);
396  Double_t b3 = a3 + 0.0; // box dimension .6 instead of .5 due to strong curvature
397 
398  // lens1 (b) + lens2 (0.6) + lens3 (b3) + gap (0.5)
399 
400  len = b + 0.2 + b3 + 0.5;//+ a2; // dimension of the box containing both lenses
401 
402  cout<<"DIRC len= "<<len<<endl;
403 
404  fSlabEnd = -bbox_hlen + bbox_shift + len; // used in processHits
405  cout<<"bar ends at = "<<fSlabEnd<<endl;
406 
407  // Lenses
408 
409  // Lens 1
410  Double_t t = -r +b/2;
411 
412  TGeoSphere* logicSphere= new TGeoSphere("S",0.,r, 0. ,180.,0.,360.);
413  TGeoBBox* lBox = new TGeoBBox("B", (bbX/barnum)/2-bargap, hthick, b/2.);
414  TGeoTranslation *tr1 = new TGeoTranslation("tr1", 0.,0., t);
415  tr1->RegisterYourself();
416  TGeoCompositeShape *cs = new TGeoCompositeShape("cs","S*(B:tr1)");
417  TGeoVolume *lens1 = new TGeoVolume("DrcLENS1Sensor",cs, gGeoManager->GetMedium("FusedSil"));
418  lens1->SetLineColor(kRed-8);
419  lens1->SetTransparency(40);
420 
421  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r + b
422  // with -0.01 one can make a gap visible (.1mm) for orientation
423  //barContainer->AddNode(lens1, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r+b) /*-0.01*/, new TGeoRotation (0)));
424 /* barContainer->AddNode(lens1, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r+b), new TGeoRotation (0)));
425 */ fdz_lens1 = -(bbox_hlen)+len -(-r+b);
426 
427 
428  //Lens 2
429  Double_t t2 = -r2;// +b2/2 r2 is the reference point (concave lens)
430  TGeoSphere* logicSphere2 = new TGeoSphere("S2",0 ,r2, 0. ,180.,0.,360.);
431  TGeoBBox* lBox2 = new TGeoBBox("B2", (bbX/barnum)/2.-bargap, hthick, b2/2.); // lside -> bbX
432  TGeoTranslation *tr2 = new TGeoTranslation("tr2", 0.,0., t2);
433  tr2->RegisterYourself();
434  TGeoCompositeShape *cs2 = new TGeoCompositeShape("cs2","(B2:tr2)-S2");
435  TGeoVolume *lens2 = new TGeoVolume("DrcLENS2Sensor",cs2, gGeoManager->GetMedium("NLAK33A"));
436  lens2->SetLineColor(kRed+2);
437  lens2->SetTransparency(40);
438 
439  // place the lens exactly on lens1
440  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r2
441  // the tip of lens1 is at b
442  // with -0.02 one can make a gap visible (.1mm due to lens1) for orientation
443  //barContainer->AddNode(lens2, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r2) -b /*-0.02*/, new TGeoRotation (0)));
444 /* barContainer->AddNode(lens2, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r2) -b , new TGeoRotation (0)));
445 */ fdz_lens2 = -(bbox_hlen)+len -(-r2) -b;
446 
447  //Lens3 (like lens1, same treatment)
448  Double_t t3 = -r3+b3/2;
449  TGeoSphere* logicSphere3= new TGeoSphere("S3",0.,r3, 0. ,180.,0.,360.);
450  TGeoBBox* lBox3 = new TGeoBBox("B3", (bbX/barnum)/2-bargap, hthick, b3/2.); // lside -> bbX
451  TGeoTranslation *tr3 = new TGeoTranslation("tr3", 0.,0., t3);
452  tr3->RegisterYourself();
453  TGeoCompositeShape *cs3 = new TGeoCompositeShape("cs3","S3*(B3:tr3)");
454  TGeoVolume *lens3 = new TGeoVolume("DrcLENS3Sensor",cs3, gGeoManager->GetMedium("NLAK33A"));
455  lens3->SetLineColor(kRed-6);
456  lens3->SetTransparency(40);
457 
458  // place the lens exactly on lens2 plane side
459  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r3 + b3
460  // with -0.03 one can make a gap visible (.1mm due to lens 1&2) for orientation
461  // b2/2 is the thickness of lens2 in the middle
462  //barContainer->AddNode(lens3, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r3+b3) -b - b2/2 /*-0.03*/, new TGeoRotation (0)));
463 /* barContainer->AddNode(lens3, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r3+b3) -b - b2/2 , new TGeoRotation (0)));
464  //cout<<" DIRC r,r2,r3 = "<<r<<" "<<r2<<" "<<r3<<endl;
465 */ fdz_lens3 = -(bbox_hlen)+len -(-r3+b3) -b - b2/2. ;
466 
467  //fAtBarEnd = "DrcLENS3";
468  } // E N D O F N E W L E N S E S
469 
470  // 2 - mirror:
471  if (fFocusingSystem == 2){ // Mirrors at front
472  // Put some mirrors at the downstream end with a focal plane at the PD.
473 
474  double zpos=120.;
475 
476  // The bar is produced with "bbox_hlen-fabs(len)/2.-mirr_hthick"
477  len = -130. + zpos + 1.;
478  //len = -247.0; // negative to make space at downstream end of bar.
479 
480  Double_t len1 = 1.5; // 1st block
481 
482  Double_t mirror_angle = 0.0;// 0.0 is pointing upstream, -90.0 is pointing to the beam axis
483  Double_t focal_length = 2. * bbox_hlen + sob_len - fabs(len) + len1;
484  Double_t mirror_radius = 2. * focal_length;
485 
486  cout<<" mirror radius: "<<mirror_radius<<endl;
487 
488 
489  // no angle for the moment
490  // block
491  TGeoSphere* logicSphere = new TGeoSphere("S",0.,mirror_radius, 0. ,180.,0.,360.);
492  TGeoBBox* lBox = new TGeoBBox("B", (bbX/barnum)/2-bargap, hthick, fabs(len1)/2.); // lside -> bbX
493 
494  Double_t t = mirror_radius - len1/2.;
495 
496  TGeoTranslation *tr1 = new TGeoTranslation("tr1", 0.,0., t);
497  tr1->RegisterYourself();
498  TGeoCompositeShape *cs = new TGeoCompositeShape("cs","S*(B:tr1)");
499 
500  TGeoVolume *block1 = new TGeoVolume("DrcBlock1",cs, gGeoManager->GetMedium("FusedSil"));
501  block1->SetLineColor(kRed);
502  block1->SetTransparency(40);
503 
504 
505  Double_t shift1 = len1-mirror_radius; // Now the start of block1 is at zero
506  shift1 += bbox_hlen-fabs(len)-2.*mirr_hthick;
507 
508 /*
509  barContainer->AddNode(block1, 1, new TGeoCombiTrans(0., 0., shift1, new TGeoRotation (0)));
510 */ fdz_mirr1 = shift1;
511 
512 
513  Double_t gap = 0;
514  Double_t len2 = fabs(len) - len1 - gap;
515 
516  // no angle for the moment
517  // block
518  TGeoSphere* logicSphere2 = new TGeoSphere("S2",0,mirror_radius, 0. ,180.,0.,360.);
519  TGeoBBox* lBox2 = new TGeoBBox("B2", (bbX/barnum)/2-bargap, hthick, fabs(len2)/2.); // lside -> bbX
520 
521  // make radius part of block
522  Double_t t2 = mirror_radius + len2/2 - 1;
523 
524  TGeoTranslation *tr2 = new TGeoTranslation("tr2", 0.,0., t2);
525  tr2->RegisterYourself();
526  TGeoCompositeShape *cs2 = new TGeoCompositeShape("cs2","(B2:tr2)-S2");
527 
528  TGeoVolume *block2 = new TGeoVolume("DrcBlock2",cs2, gGeoManager->GetMedium("Mirror"));
529  block2->SetLineColor(kGreen);
530  block2->SetTransparency(40);
531 
532 
533  Double_t shift2 = -mirror_radius; // Now the start of the block is at zero
534  shift2 += bbox_hlen-fabs(len)-2*mirr_hthick; // at end of bar
535  shift2 += len1 + gap;
536 
537  fdz_mirr2 = shift2;
538 
539 /*
540  barContainer->AddNode(block2,
541  1,
542  new TGeoCombiTrans(0.,
543  0.,
544  shift2,
545  // bbox_hlen-mirror_radius-2*mirr_hthick-len2+len1+1,
546  //bbox_hlen-mirror_radius-len2+0.1,
547  new TGeoRotation (0)
548  )
549  );
550 */
551 
552 
553  Double_t flen = 0.;
554  fSlabEnd = -bbox_hlen + bbox_shift + flen;
555  cout<<"bar ends at = "<<fSlabEnd<<endl;
556 
557  //fAtBarEnd = "DrcBar";
558  } // E N D O F MIRRORS
559 
560  if(fFocusingSystem == 3){ // cylindrical lens w/o airgap
561 
562  // main parameters:
563 
564  Double_t Hcyl2 = 0.1; //[cm] thickness in the middle of the second (FusedSil) lens
565  Double_t Rcyl = 7.35; // [cm]
566  Double_t Lcyl = (bbX/barnum)/2.-bargap; //cylinder length/2.
567  Double_t Acyl = TMath::ASin(hthick/Rcyl); // angle
568  Double_t Tcyl1 = 0.5*Rcyl*(1. + TMath::Cos(Acyl)); // distance between centers of the cylinder and the box
569 
570  len = Hcyl2 + Rcyl*(1. - TMath::Cos(Acyl));
571  cout<<"DIRC: len = "<<len<<endl;
572 
573  cout<<"length = "<< Lcyl<<", angle = "<<Acyl/3.1415*180.<<", translation = "<<Tcyl1<<endl;
574  cout<<"lens 0.5*width = "<<hthick<<", 0.5*thickness = "<< 0.5*Rcyl*(1.-TMath::Cos(Acyl))<<", 0.5*heigth = "<<Lcyl <<endl;
575 
576  //Lens1
577  TGeoEltu* lCylinder = new TGeoEltu("Cyl",Rcyl,Rcyl,Lcyl);
578  TGeoBBox* lCylBox = new TGeoBBox("CylBox", hthick, hthick, Lcyl);
579  //TGeoBBox* lCylBox = new TGeoBBox("CylBox",36.75, 36.75, Lcyl);
580  //TGeoTranslation *trCyl = new TGeoTranslation("trCyl", 0., 100.4, 0.);
581  TGeoTranslation *trCyl1 = new TGeoTranslation("trCyl", 0., Rcyl*TMath::Cos(Acyl)+hthick, 0.);// !!!
582  trCyl1->RegisterYourself();
583  TGeoCompositeShape *llens1 = new TGeoCompositeShape("llens1","Cyl*(CylBox:trCyl)");
584  TGeoVolume *CylLens1 = new TGeoVolume("DrcLENS1Sensor",llens1, gGeoManager->GetMedium("NLAK33A"));
585  CylLens1->SetLineColor(kRed-8);
586  CylLens1->SetTransparency(40);
587 
588  fdz_lens1 = -(bbox_hlen) + len - (-Rcyl + Rcyl*(1.-TMath::Cos(Acyl)));
589 
590  //Lens2
591  TGeoTranslation *trCyl2 = new TGeoTranslation("trCyl2", 0., Rcyl + Hcyl2 - hthick, 0.);
592  trCyl2->RegisterYourself();
593  TGeoCompositeShape *llens2 = new TGeoCompositeShape("llens2", "(CylBox:trCyl2) - Cyl");
594  TGeoVolume* CylLens2 = new TGeoVolume("DrcLENS2Sensor", llens2, gGeoManager->GetMedium("FusedSil"));
595  CylLens2->SetLineColor(kRed+2);
596  CylLens2->SetTransparency(40);
597 
598  // TGeoRotation rot_lens3;
599  // rot_lens3.RotateZ(-90.);
600  // rot_lens3.RotateX( 90.);
601 
602 
603 
604  fdz_lens2 = -(bbox_hlen) + len + (hthick - Rcyl*(1.-TMath::Cos(Acyl)) - Hcyl2);
605  fSlabEnd = -bbox_hlen + bbox_shift;
606  cout<<"bar ends at = "<<fSlabEnd<<endl;
607 
608  }
609  //---------------------------------------------------------------------------------------------
610 
611 
612  // create top volume:
613  TGeoVolume* cave;
614  TGeoBBox* lTop = new TGeoBBox(500,500,300);
615  cave = new TGeoVolume("DIRC", lTop, air_m);
616  gGeoManager->SetTopVolume(cave);
617 
618  // create pre-top volume:
619  TGeoVolume* vLocalMother;
620  TGeoPcon* shape = new TGeoPcon("BarrelDIRCShape", 0, 360., 4);
621  shape->DefineSection(0, bbox_zdown, 35., 60.);
622  shape->DefineSection(1, bbox_zup, 35., 60.);
623  shape->DefineSection(2, bbox_zup - sob_len, 35., sob_Rout+poffset+pheight+EVoffset+1.);
624  shape->DefineSection(3, bbox_zup - sob_len - PDbaseLayer - EVgreaseLayer, 35., sob_Rout+poffset+pheight+EVoffset+1.);
625  vLocalMother = new TGeoVolume("BarrelDIRC", shape, air_m); //DIRCairNoSens_m); //################
626  cave->AddNode(vLocalMother, 0,0);
627 
628  cout<<"bbox length = "<<2.*(bbox_hlen-0.5*(boxgap+boxthick))<<", prizm length = "<<2.*(phlength+0.5*(boxgap+boxthick))<<endl;
629  cout<<"prizm shift = "<<-bbox_hlen+0.5*(boxgap+boxthick)-phlength<<endl;
630 
631  // create BarBoxes:
632  TGeoBBox* logicbbL;
633  TGeoVolume *bbox;
634  logicbbL = new TGeoBBox("logicbbL", bbX/2.+boxthick, hthick+boxgap+boxthick, bbox_hlen);
635  bbox = new TGeoVolume("DrcBarBox", logicbbL,DIRCcarbonFiber_m);
636  bbox->SetLineColor(30);
637 
638 
639  TGeoBBox* logicbbS;
640  TGeoVolume *abox;
641 
642  logicbbS = new TGeoBBox("logicbbS", bbX/2., hthick+boxgap, bbox_hlen-barWin_hthick);
643  abox = new TGeoVolume("DrcAirBox", logicbbS, DIRCairNoSens_m);
644  bbox->AddNode(abox, 1, new TGeoCombiTrans(0., 0., barWin_hthick, new TGeoRotation(0)));
645  abox->SetLineColor(19); // grau
646 
647  //create windows at the readout end of the bar boxes:
648  TGeoBBox* logicBarWin = new TGeoBBox("logicBarWin", bbX/2., hthick+boxgap, barWin_hthick);
649  TGeoVolume* barwin = new TGeoVolume("DrcBarboxWindow", logicBarWin, FusedSil_m);
650  barwin->SetLineColor(kBlue-4);
651  bbox->AddNode(barwin, 1, new TGeoCombiTrans(0.,0.,-bbox_hlen+barWin_hthick,new TGeoRotation(0)));
652 
653  //create layers of grease at the readout end of the bar boxes, between the windows and the EV:
654  TGeoBBox* logicEVgrease = new TGeoBBox("logicEVgrease", bbX/2., hthick+boxgap, EVgreaseLayer/2.);
655  TGeoVolume* evgrease = new TGeoVolume("DrcEVgrease", logicEVgrease, OpticalGrease_m);
656  evgrease->SetLineColor(kSpring);
657 
658  // put barboxes into right positions:
659  Double_t dx_bbox, dy_bbox, dz_bbox, phi_curr;
660 
661  for(Int_t m = 0; m < bbnum; m ++){
662  phi_curr = (90. - phi0 - dphi*m)/180.*pi;
663  if(m > bbnum/2-1){ phi_curr = (90. - phi0 - dphi*m - 2.*pipehAngle)/180.*pi; }
664  dx_bbox = radius * cos(phi_curr);
665  dy_bbox = radius * sin(phi_curr);
666  dz_bbox = bbox_shift;
667 
668  TGeoRotation rot_bbox;
669  rot_bbox.RotateZ( -phi0 - m*dphi - (TMath::Floor(2.*m/bbnum))*(2.*pipehAngle));
670  vLocalMother->AddNode(bbox, m+1, new TGeoCombiTrans(dx_bbox, dy_bbox, dz_bbox, new TGeoRotation(rot_bbox)));
671  vLocalMother->AddNode(evgrease, m+1, new TGeoCombiTrans(dx_bbox, dy_bbox, bbox_zup-EVgreaseLayer/2.,new TGeoRotation(rot_bbox)));
672  }
673 
674  cout<<"bar width = "<<2.*(((bbX/barnum)/2.)-bargap) << ", bar with gaps = "<<(bbX/barnum)<<endl;
675  cout<<"barboxL width = "<<bbX/2.+boxgap+boxthick<<", barboxS width = "<<bbX/2.+boxgap<<endl;
676 
677  // create logic bar:
678  TGeoBBox* logicBar = new TGeoBBox("logicBar", ((bbX/barnum)/2.)-bargap, hthick, bbox_hlen-fabs(len)/2.-mirr_hthick-barWin_hthick);
679  TGeoVolume* bar;
680  bar = new TGeoVolume("DrcBarSensor",logicBar, FusedSil_m );
681  bar->SetLineColor(kCyan-9);
682  bar->SetTransparency(60);
683 
684  // create logic mirror:
685  TGeoBBox* logicMirror = new TGeoBBox("logicMirror", bbX/barnum/2.-bargap, hthick, mirr_hthick /2.); //$$$$$$$
686  TGeoVolume *mirr = new TGeoVolume("DrcMirr", logicMirror, Mirror_m);
687  mirr->SetLineColor(5);
688 
689  // create glue layer inside the bar (connects two halves):
690  TGeoBBox* logicBarGlue = new TGeoBBox("logicBarGlue", bbX/barnum/2.-bargap, hthick, gluehthick); //$$$$$$$
691  TGeoVolume *barglue = new TGeoVolume("DrcBarGlue", logicBarGlue, Epotek301_2_m);
692  barglue->SetLineColor(kSpring-5);
693  bar->AddNode(barglue, 1, new TGeoCombiTrans(0., 0., 0., new TGeoRotation (0)));
694 
695  Double_t dx, dy, dz_bar, dz_mirr;
696 
697  for(Int_t j=0; j<barnum; j++){
698  dx = - (bbX/2.) + (bbX/barnum)/2. + j * (bbX/barnum);
699  dy = 0.;//len/2. - mirr_hthick;
700  dz_bar = len/2.-mirr_hthick;
701  dz_mirr = bbox_hlen - barWin_hthick - mirr_hthick;
702  if(fFocusingSystem == 1){ // lens
703  abox->AddNode(lens1, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation (0)));
704  abox->AddNode(lens2, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens2, new TGeoRotation (0)));
705  abox->AddNode(lens3, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens3, new TGeoRotation (0)));
706  }
707  if(fFocusingSystem == 2){ // mirror
708  abox->AddNode(block1, 1+j, new TGeoCombiTrans(dx, dy, fdz_mirr1, new TGeoRotation (0)));
709  abox->AddNode(block2, 1+j, new TGeoCombiTrans(dx, dy, fdz_mirr2, new TGeoRotation (0)));
710  }
711  if(fFocusingSystem == 3){
712  TGeoRotation rot_lens3;
713  rot_lens3.RotateZ(-90.);
714  rot_lens3.RotateY( 90.);
715  abox->AddNode(CylLens1, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation(rot_lens3)));
716  abox->AddNode(CylLens2, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation(rot_lens3)));
717  }
718 
719  abox->AddNode(bar, 1+j, new TGeoCombiTrans(dx, dy, dz_bar, new TGeoRotation(0)));
720  if(fFocusingSystem != 2){ // not a forward mirror
721  abox->AddNode(mirr, 1+j, new TGeoCombiTrans(dx, dy, dz_mirr, new TGeoRotation(0))); //$$$$$$$$$$$$$$$$$
722  }
723  }
724 
725  // Expansion volume:
726  Double_t dR;
727  Double_t xEV;
728  Double_t cosFactor1;
729 
730  TGeoPgon* logicEV1, * logicEV2, *logicEV3, * logicEV4;
731  cosFactor1 = cos(pipehAngle/180.*pi)/cos(dphi/180.*pi/2.);
732  dR = (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi) - (radius-hthick);
733  xEV = (dR + sob_len*tan(sob_angle/180.*pi))/ (tan(sob_angle/180.*pi) + tan(sob_angleB/180.*pi));
734  if(sob_angleB == 90.){
735  logicEV1 = new TGeoPgon("logicEV1", 93.6, 172.8, bbnum/2, 2);
736  logicEV1->DefineSection(0, 0., radiusMiddleSmall, sob_Rout);
737  logicEV1->DefineSection(1, sob_len, radiusMiddleSmall, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi));
738  logicEV2 = new TGeoPgon("logicEV2", -86.4, 172.8, bbnum/2, 2);
739  logicEV2->DefineSection(0, 0., radiusMiddleSmall, sob_Rout);
740  logicEV2->DefineSection(1, sob_len, radiusMiddleSmall, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi));
741  logicEV3 = new TGeoPgon("logicEV3", 86.4, 7.2, 1, 2);
742  logicEV3->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1, sob_Rout*cosFactor1);
743  logicEV3->DefineSection(1, sob_len, (radiusMiddleSmall)*cosFactor1, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi)*cosFactor1);
744  logicEV4 = new TGeoPgon("logicEV4", -93.6, 7.2, 1, 2);
745  logicEV4->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1, sob_Rout*cosFactor1);
746  logicEV4->DefineSection(1, sob_len, (radiusMiddleSmall)*cosFactor1, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi)*cosFactor1);
747 
748  }
749  if(sob_angleB != 90.){
750  logicEV1 = new TGeoPgon("logicEV1", 93.6, 172.8, bbnum/2, 3);
751  logicEV1->DefineSection(0, 0., radius-hthick, radius-hthick+eps);
752  logicEV1->DefineSection(1, xEV, radius-hthick, sob_Rout - xEV*tan(sob_angle/180.*pi));
753  logicEV1->DefineSection(2, sob_len, radius-hthick, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi));
754  logicEV2 = new TGeoPgon("logicEV2", -86.4, 172.8, bbnum/2, 3);
755  logicEV2->DefineSection(0, 0., radius-hthick, radius-hthick+eps);
756  logicEV2->DefineSection(1, xEV, radius-hthick, sob_Rout - xEV*tan(sob_angle/180.*pi));
757  logicEV2->DefineSection(2, sob_len, radius-hthick, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi));
758  logicEV3 = new TGeoPgon("logicEV3", 86.4, 7.2, 1, 3);
759  logicEV3->DefineSection(0, 0., (radius-hthick)*cosFactor1, (radius-hthick+eps)*cosFactor1);
760  logicEV3->DefineSection(1, xEV, (radius-hthick)*cosFactor1, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1);
761  logicEV3->DefineSection(2, sob_len, (radius-hthick)*cosFactor1, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi)*cosFactor1);
762  logicEV4 = new TGeoPgon("logicEV4", -93.6, 7.2, 1, 3);
763  logicEV4->DefineSection(0, 0., (radius-hthick)*cosFactor1, (radius-hthick+eps)*cosFactor1);
764  logicEV4->DefineSection(1, xEV, (radius-hthick)*cosFactor1, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1);
765  logicEV4->DefineSection(2, sob_len, (radius-hthick)*cosFactor1, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi)*cosFactor1);
766  }
767 
768  TGeoCompositeShape *logicEV = new TGeoCompositeShape("logicEV","logicEV1 + logicEV2 + logicEV3 + logicEV4");
769  TGeoVolume* baseEV = new TGeoVolume("DrcEVSensor", logicEV, FusedSil_m);//Marcol82_7_m); //####################################
770  baseEV->SetLineColor(kMagenta+2);
771  baseEV->SetTransparency(50);
772  vLocalMother->AddNode(baseEV, 1, new TGeoCombiTrans(0.,0.,sob_shift - EVgreaseLayer, new TGeoRotation(0)));
773 
774  // PhotoDetector Basic material - Carbon. It is placed on the back side of the EV to simulate the support structure of the MCPs, the photons are to hit it in gaps between MCPs:
775  TGeoPgon *logicPDbase1 = new TGeoPgon("logicPDbase1", 93.6, 172.8, bbnum/2, 2);
776  TGeoPgon *logicPDbase2 = new TGeoPgon("logicPDbase2", -86.4, 172.8, bbnum/2, 2);
777  TGeoPgon *logicPDbase3 = new TGeoPgon("logicPDbase3", 86.4, 7.2, 1, 2);
778  TGeoPgon *logicPDbase4 = new TGeoPgon("logicPDbase4", -93.6, 7.2, 1, 2);
779  Double_t rad_delta = (MCPsize-MCPactiveArea)/2./cos(45./180.*pi);
780  if(sob_angleB == 90.){
781  logicPDbase1->DefineSection(0, 0., radiusMiddleSmall-rad_delta, sob_Rout);
782  logicPDbase1->DefineSection(1, PDbaseLayer, radiusMiddleSmall-rad_delta, sob_Rout);
783  logicPDbase2->DefineSection(0, 0., radiusMiddleSmall-rad_delta, sob_Rout);
784  logicPDbase2->DefineSection(1, PDbaseLayer, radiusMiddleSmall-rad_delta, sob_Rout);
785  logicPDbase3->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
786  logicPDbase3->DefineSection(1, PDbaseLayer, (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
787  logicPDbase4->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
788  logicPDbase4->DefineSection(1, PDbaseLayer, (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
789  }
790 /* if(sob_angleB != 90.){
791  logicPD1->DefineSection(0, 0., radius-hthick+eps, radius-hthick+eps+0.1);
792  logicPD1->DefineSection(1, xEV, sob_Rout - xEV*tan(sob_angle/180.*pi), sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1);
793  logicPD2->DefineSection(0, 0., radius-hthick+eps, radius-hthick+eps+0.1);
794  logicPD2->DefineSection(1, xEV, sob_Rout - xEV*tan(sob_angle/180.*pi), sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1);
795  logicPD3->DefineSection(0, 0., (radius-hthick+eps)*cosFactor1, (radius-hthick+eps+0.1)*cosFactor1);
796  logicPD3->DefineSection(1, xEV, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1, (sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1)*cosFactor1);
797  logicPD4->DefineSection(0, 0., (radius-hthick+eps)*cosFactor1, (radius-hthick+eps+0.1)*cosFactor1);
798  logicPD4->DefineSection(1, xEV, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1, (sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1)*cosFactor1);
799  }*/
800 
801  TGeoCompositeShape *logicPDbase = new TGeoCompositeShape("logicPDbase","logicPDbase1 + logicPDbase2 + logicPDbase3 + logicPDbase4");
802  TGeoVolume *pdbase = new TGeoVolume("DrcPDbase", logicPDbase, DIRCcarbonFiber_m);
803  pdbase->SetLineColor(kGreen-6);
804  if(sob_angleB == 90.){
805  vLocalMother->AddNode(pdbase, 1, new TGeoCombiTrans(0., 0., sob_shift-EVgreaseLayer-PDbaseLayer, new TGeoRotation(0)));
806  }
807 /* if(sob_angleB != 90.){
808  vLocalMother->AddNode(pd, 1, new TGeoCombiTrans(0., 0., sob_shift, new TGeoRotation(0)));
809  }
810 */
811  //MCPs as PD
812  TGeoBBox* logicPixel = new TGeoBBox("logicPixel", PixelSize/2., PixelSize/2., PDsensitiveThick/2.);
813  TGeoVolume* pix = new TGeoVolume("DrcPDSensor",logicPixel, FusedSil_m);
814  pix->SetLineColor(kGreen+1);
815 
816  // create pixel holders:
817  TGeoBBox* logicPixelHolder = new TGeoBBox("logicPixelHolder", MCPactiveArea/2., MCPactiveArea/2., PDsensitiveThick/2.);
818  TGeoVolume *pixelholder = new TGeoVolume("DrcPixelHolder", logicPixelHolder, DIRCcarbonFiber_m);
819  pixelholder->SetLineColor(kOrange-8);
820 
821  //put pixels into holders:
822  for(Int_t itr=0; itr<Npix; itr++){
823  Double_t dy_row = (-Npix/2 + itr)*(PixelSize + PixelGap) + 0.5*(PixelSize+PixelGap);
824  for(Int_t jtr=0; jtr<Npix; jtr++){
825  Double_t dx_row = (-Npix/2 + jtr)*(PixelSize + PixelGap) + 0.5*(PixelSize+PixelGap);
826  pixelholder->AddNode(pix, itr*Npix+jtr, new TGeoCombiTrans(dx_row, dy_row, 0., new TGeoRotation(0)));
827  }
828  }
829 
830  // create photo cathodes for each MCP
831  TGeoBBox* logicPhCathode = new TGeoBBox("logicPhCathode", MCPactiveArea/2., MCPactiveArea/2., PhCathodeThick/2.);
832  TGeoVolume *phcathode = new TGeoVolume("DrcPhCathodeSensor", logicPhCathode, FusedSil_m);
833  phcathode->SetLineColor(kGray+1);
834 
835  // create Windows for each MCP
836  TGeoBBox* logicWindow = new TGeoBBox("logicWindow", MCPsize/2., MCPsize/2., PDwindowThick/2.);
837  TGeoVolume *window = new TGeoVolume("DrcPDwindowSensor", logicWindow, FusedSil_m);
838  window->SetLineColor(kBlue-4);
839 
840  // create grease layers between MCP window and the EV back side
841  TGeoBBox* logicMCPgrease = new TGeoBBox("logicMCPgrease", MCPsize/2., MCPsize/2., PDgreaseLayer/2.);
842  TGeoVolume *mcpgrease = new TGeoVolume("DrcMcpGrease", logicMCPgrease, OpticalGrease_m);
843  mcpgrease->SetLineColor(kSpring);
844 
845  //version 1: place MCPs in global cs
846  Double_t hgap = 0.5*(MCPsize - MCPactiveArea + MCPgap);
847  Double_t step = MCPactiveArea + 2.*hgap;
848  Double_t xpos = 0.;
849  Double_t ypos = 0.;
850  Double_t xcurr = 0.;
851  Double_t ycurr = 0.;
852  Double_t xsmall, xlarge, lastxsmall;
853  Int_t i=0;
854  cout<<"hgap = "<<hgap<<", step = "<<step<<endl;
855  xpos = radiusCornerIn + 0.5*MCPactiveArea;
856  ypos = step/2.;
857  int cas1 = 1;
858  int cas2 = 1;
859  Int_t TotalNmcp = 0;
860  for(Int_t nquat =0; nquat<4; nquat++){
861  for(Int_t i=0; i<Int_t((radiusCornerOut)/step); i++){ // loop over y
862  cas1 = findSectorIn(ypos-step/2., dphi_rad, radius, EVdrop, hthick);
863  cout<<"+++++++++++++++ nquat = "<<nquat<<", i = "<<i<<", xpos = "<<xpos<<", ypos = "<<ypos<<endl;
864  cout << "calculated sector1 from y =" << ypos << " is " << cas1 << endl;
865  getXrange(dphi, radiusCornerOut+1.3, radiusCornerIn, MCPactiveArea, cas1, cas2, ypos-(step/2.), step, xsmall, xlarge);
866  xpos = xsmall;
867  if(xsmall < xlarge){
868  for(Int_t j=0; j<TMath::Ceil((xlarge-xsmall)/step); j++){ // loop over x
869  cout<<"j = "<<j<<". quat = "<<nquat<<", xsmall = "<<xsmall<<", xlarge = "<<xlarge<<", cas1 = "<<cas1<<", limit = "<<(Int_t((xlarge-xsmall)/step))+1<<endl;
870 
871  if(nquat == 1 || nquat == 2){
872  xcurr = -xpos;
873  }
874  if(nquat == 0 || nquat == 3){
875  xcurr = xpos;
876  }
877  if(nquat == 0 || nquat == 1){
878  ycurr = ypos;
879  }
880  if(nquat == 2 || nquat == 3){
881  ycurr = -ypos;
882  }
883  // put Photo Cathodes
884  pdbase->AddNode(phcathode, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(0)));
885  //put windows
886  pdbase->AddNode(window, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(0)));
887  //put grease layers
888  pdbase->AddNode(mcpgrease, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(0)));
889  //put pixel holders:
890  pdbase->AddNode(pixelholder, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PDwindowThick-PhCathodeThick-PDsensitiveThick/2., new TGeoRotation(0)));
891  TotalNmcp = TotalNmcp + 1;
892  xpos = xpos + step;
893  //if(xpos >= xlarge){continue;}
894  }
895  cas2 = findSectorOut(ypos+step/2., dphi_rad, radiusCornerOut);
896  cout<<"calculated sector2 from y = "<<ypos+step/2. << " is "<< cas2<<endl;
897  ypos = ypos + step;
898  if(i == 12){ypos = step/2.;cas2 = 1;}
899  } // if
900  }
901  } // for nquat
902  cout<<"version 1: totally "<<TotalNmcp+1<<" mcps on the PD plane"<<endl;
903 
904 /* //version 2: place MCPs in local bar coordinate systems:
905  //create a map of MCPs for one sector (in bar c.s.):
906  Double_t sectorWidth = 0.;
907  Int_t nmcp = 0;
908  Int_t totalnumbering = 1;
909  TVector3 location;
910  Double_t phi_curr1 = 0.;
911  for(Int_t m = 0; m < bbnum; m ++){
912  phi_curr1 = (90. - phi0 - dphi*m)/180.*pi;
913  if(m > bbnum/2-1){ phi_curr1 = (90. - phi0 - dphi*m - 2.*pipehAngle)/180.*pi; }
914  //dx_bbox = radius * cos(phi_curr1);
915  //dy_bbox = radius * sin(phi_curr1);
916  TGeoRotation rot_sector;
917  rot_sector.RotateZ( -phi0 - m*dphi - (TMath::Floor(2.*m/bbnum))*(2.*pipehAngle));
918  cout<<"phi_curr1 = "<<phi_curr1/3.1415*180.<<endl;
919  // placement of MCPs in one sector
920  cout<<"rad corner out = "<<radiusCornerOut<<", step = "<<step<<endl;
921  for(Int_t nrow=0; nrow < Int_t((radiusCornerOut*cos(dphi/2./180.*pi)-radiusMiddleSmall)/step); nrow++){
922  sectorWidth = 2.* (radiusMiddleSmall + step*nrow) * tan(dphi_rad/2.);
923  nmcp = Int_t(sectorWidth/(MCPsize + 0.5*MCPgap));
924  cout<<"sector width = "<<sectorWidth<<", nmcp = "<<nmcp<<", nmcp%2 = "<<nmcp%2<<endl;
925  xpos = (radiusMiddleSmall + 0.5*MCPactiveArea + step*(nrow));
926  for(Int_t ny=0; ny<nmcp; ny++){
927  ypos = ((-Int_t(nmcp/2.) - 0.5*(nmcp%2))*step + (0.5+ny)*step);
928  location.SetXYZ(xpos,ypos,0.);
929  location.RotateZ(phi_curr1);
930  pdbase->AddNode(mcpgrease, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(rot_sector)));
931  pdbase->AddNode(window, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(rot_sector)));
932  pdbase->AddNode(phcathode, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(rot_sector)));
933  pdbase->AddNode(pixelholder, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick-PDwindowThick-PDsensitiveThick/2., new TGeoRotation(rot_sector)));
934  totalnumbering = totalnumbering + 1;
935  }
936  //vLocalMother->AddNode(mcp, m, new TGeoCombiTrans(dx_bbox, dy_bbox, sob_shift-PDthick/2., new TGeoRotation(rot_sector)));
937  }
938  } //for m
939  for(Int_t nrow=0; nrow < Int_t((radiusCornerOut*cos(dphi/2./180.*pi)-radiusMiddleSmall)/step); nrow++){
940  for(Int_t nadd = 0; nadd<2; nadd++){
941  xpos = (radiusMiddleSmall*cosFactor1 + 0.5*MCPactiveArea + step*(nrow));
942  location.SetXYZ(xpos,0.,0.);
943  location.RotateZ(pi/2.+pi*nadd);
944  pdbase->AddNode(mcpgrease, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(0)));
945  pdbase->AddNode(window, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(0)));
946  pdbase->AddNode(phcathode, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(0)));
947  pdbase->AddNode(pixelholder, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick-PDwindowThick-PDsensitiveThick/2., new TGeoRotation(0)));
948  totalnumbering = totalnumbering + 1;
949  }
950  }
951 */
952  gGeoManager->CloseGeometry();
953 
954  cave->CheckOverlaps(0.0001, "");
955  gGeoManager->CheckOverlaps(0.00001,""); // [cm]
956  //gGeoManager->CheckGeometryFull();
957 
958  cave->Write();
959  fi->Close();
960  cave->Draw("ogl");
961  //pdbase->Draw("ogl");
962 
963  TObjArray *listOfOverlaps = gGeoManager->GetListOfOverlaps();
964  cout<<listOfOverlaps->GetEntries()<<endl;
965  //listOfOverlaps->Print();
966  return 0;
967 }
968 
969 int findSectorIn(double y, double dphi_rad, double radius, double EVdrop, double hthick){
970  TVector3 EVcornerIn;
971  EVcornerIn.SetXYZ((radius-EVdrop-hthick)/cos(dphi_rad/2.),0.,0.);
972  Double_t Y0[4];
973  for(Int_t i=0; i<4; i++){
974  EVcornerIn.RotateZ(dphi_rad);
975  Y0[i] = EVcornerIn.Y();
976  }
977  if(y <= Y0[0]){return 1;}
978  if(y > Y0[0] && y <= Y0[1]){return 2;}
979  if(y > Y0[1] && y <= Y0[2]){return 3;}
980  if(y > Y0[2] && y <= Y0[3]){return 4;}
981  if(y > Y0[3]){return 5;}
982 
983 }
984 
985 int findSectorOut(double y, double dphi_rad, double radiusCornerOut){
986  TVector3 EVcorner;
987  EVcorner.SetXYZ(radiusCornerOut,0.,0.);
988  Double_t Y0[4];
989  for(Int_t i=0; i<4; i++){
990  EVcorner.RotateZ(dphi_rad);
991  Y0[i] = EVcorner.Y();
992  cout<<"Y ["<<i<<"] = "<<Y0[i]<<endl;
993  }
994  if(y <= Y0[0]){return 1;}
995  if(y > Y0[0] && y <= Y0[1]){return 2;}
996  if(y > Y0[1] && y <= Y0[2]){return 3;}
997  if(y > Y0[2] && y <= Y0[3]){return 4;}
998  if(y > Y0[3]){return 5;}
999 
1000 }
1001 
1002 bool getXrange(double dphi, double radiusCornerOut, double radiusCornerIn, double MCPactiveArea, int cas1, int cas2, double ylow, double step, double &x1, double &x2){
1003 
1004  Double_t dphi_rad = dphi/180.*pi;
1005 
1006  TVector3 EVcorner;
1007  EVcorner.SetXYZ(radiusCornerOut,0.,0.);
1008  TVector3 EVcorner1 = EVcorner;
1009  TVector3 EVcorner2 = EVcorner;
1010  TVector3 EVcornerIn;
1011  EVcornerIn.SetXYZ(radiusCornerIn,0.,0.);
1012  TVector3 EVcornerIn1 = EVcornerIn;
1013  TVector3 EVcornerIn2 = EVcornerIn;
1014  //cout<<"-I- sob_Rout = "<<sob_Rout<<", dphi = "<<dphi<<", |EVcorner| = "<<EVcorner.Mag()<<", |EVcornerIn| = "<<EVcornerIn.Mag()<<endl;
1015 
1016  Double_t xstart, xfinish;
1017  if(cas1 == 1){ // 0 - 21.6 degrees
1018  EVcornerIn2.RotateZ(dphi_rad);
1019  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
1020  }
1021  if(cas1 == 2){ // 21.6 - 43.2 degrees
1022  EVcornerIn1.RotateZ(1.*dphi_rad);
1023  EVcornerIn2.RotateZ(2.*dphi_rad);
1024  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
1025  }
1026  if(cas1 == 3){ // 43.2 - 63.8 degrees
1027  EVcornerIn1.RotateZ(2.*dphi_rad);
1028  EVcornerIn2.RotateZ(3.*dphi_rad);
1029  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
1030  }
1031  if(cas1 == 4){ // 63.8 - 85.6 degrees
1032  EVcornerIn1.RotateZ(3.*dphi_rad);
1033  EVcornerIn2.RotateZ(4.*dphi_rad);
1034  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
1035  }
1036  if(cas1 == 5){ // 85.6 - 90.
1037  EVcornerIn1.RotateZ(4.*dphi_rad);
1038  EVcornerIn2.RotateZ(pi/2.);
1039  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
1040  }
1041 
1042  if(cas2 == 1){ // 0 - 21.6 degrees
1043  EVcorner2.RotateZ(dphi_rad);
1044  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
1045  }
1046  if(cas2 == 2){ // 21.6 - 43.2 degrees
1047  EVcorner1.RotateZ(1.*dphi_rad);
1048  EVcorner2.RotateZ(2.*dphi_rad);
1049  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
1050  }
1051  if(cas2 == 3){ // 43.2 - 63.8 degrees
1052  EVcorner1.RotateZ(2.*dphi_rad);
1053  EVcorner2.RotateZ(3.*dphi_rad);
1054  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
1055  }
1056  if(cas2 == 4){ // 63.8 - 85.6 degrees
1057  EVcorner1.RotateZ(3.*dphi_rad);
1058  EVcorner2.RotateZ(4.*dphi_rad);
1059  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
1060  }
1061  if(cas2 == 5){ // 85.6 - 90.
1062  EVcorner1.RotateZ(4.*dphi_rad);
1063  EVcorner2.RotateZ(pi/2.);
1064  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
1065  }
1066  if(xstart <= 0.){xstart = step/2.;}
1067  if(xstart < step/2.){xstart = step/2.;}
1068  if(xfinish <= 0.){xfinish = step/2.;}
1069  x1 = xstart;
1070  x2 = xfinish;
1071 
1072  cout<<"-I- xstart = "<<xstart<<", xfinish = "<<xfinish<<endl;
1073  if(x1 != x2)return true;
1074  else return false;
1075 }
1076 
1077 double getX(TVector3 corner1, TVector3 corner2, double yy){
1078  //cout<<"-III- vector1"<<endl;
1079  //corner1.Print();
1080  //cout<<"-III- vector2 "<<endl;
1081  //corner2.Print();
1082  Double_t b = ( corner1.X()*corner2.Y() - corner1.Y()*corner2.X() ) / ( corner1.X() - corner2.X() );
1083  Double_t k = ( corner1.Y() - corner2.Y() ) / ( corner1.X() - corner2.X() );
1084  //cout<<"-II- b = "<<b<<", k = "<<k<<endl;
1085  if(( yy - b )/k >= 0.) return ( yy - b )/k;
1086  if(( yy - b )/k < 0.) return 0.;
1087 }
static T ASin(const T &x)
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
Double_t EVbackAngle()
Definition: PndGeoDrc.h:148
double dy
double r
Definition: RiemannTest.C:14
int createRootGeometry_DIRC_fsEVdroplens_MCPs(Int_t fFocusingSystem=0, Bool_t fprizm=kFALSE)
Double_t BBoxNum()
Definition: PndGeoDrc.h:136
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
TTree * b
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
double getX(TVector3 corner1, TVector3 corner2, double yy)
Double_t EVlen()
Definition: PndGeoDrc.h:127
Double_t par[3]
Double_t BBoxGap()
Definition: PndGeoDrc.h:130
Double_t EVdrop()
Definition: PndGeoDrc.h:142
Double_t boxGap()
Definition: PndGeoDrc.h:116
TGeoManager * gGeoManager
static T Cos(const T &x)
Definition: PndCAMath.h:43
Double_t barBoxZDown()
Definition: PndGeoDrc.h:104
Double_t EVoffset()
Definition: PndGeoDrc.h:145
TGeoRotation * rot1
Double_t PipehAngle()
Definition: PndGeoDrc.h:139
Double_t PrismAngle()
Definition: PndGeoDrc.h:160
Int_t a
Definition: anaLmdDigi.C:126
Double_t barBoxZUp()
Definition: PndGeoDrc.h:108
#define cave
Definition: createSTT.C:62
TGeoShape * shape
TFile * fi
Double_t PrismDrop()
Definition: PndGeoDrc.h:154
Double_t
Double_t phi0
Definition: checkhelixhit.C:60
double eps(TVector3 v1, TVector3 v2)
const Double_t zpos
int findSectorIn(double y, double dphi, double radius, double EVdrop, double hthick)
vLocalMother
Double_t z
Int_t t3
Definition: hist-t7.C:106
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Double_t boxThick()
Definition: PndGeoDrc.h:120
double dx
void getXrange(double dphi, double radiusCornerOut, double radiusCornerIn, double MCPactiveArea, int cas1, int cas2, double ylow, double step, double &x1, double &x2)
Double_t barHalfThick()
Definition: PndGeoDrc.h:96
Int_t t2
Definition: hist-t7.C:106
Double_t PrismOffset()
Definition: PndGeoDrc.h:151
Double_t PrismhLength()
Definition: PndGeoDrc.h:157
TTree * t
Definition: bump_analys.C:13
int findSectorOut(double y, double dphi_rad, double radiusCornerOut)
Double_t y
double alpha
Definition: f_Init.h:9
TString fGeoFile
Double_t barNum()
Definition: PndGeoDrc.h:124
#define air
Definition: createSTT.C:72
double r2
Double_t radius()
Definition: PndGeoDrc.h:92