FairRoot/PandaRoot
createRootGeometry_DIRC.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(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  Double_t eps = 0.01; // epsilon
52  Double_t mirr_hthick = 0.02;
53  Double_t PDbaseLayer = 5.; // [cm] thickness of the carbon at the back of the EV
54 
55  Double_t radius = fGeo->radius(); // 50. radius in middle of the barbox (x and y)
56  Double_t hthick = fGeo->barHalfThick(); // 1.7/2. half thickness of the bars
57  Double_t barnum = fGeo->barNum(); // 6 number of bars per barbox
58  Double_t bbnum = fGeo->BBoxNum(); //16. total number of sides = barboxes
59  Double_t bbGap = fGeo->BBoxGap(); //1.5 gap btw the neighboring barboxes (at the middle height)
60  Double_t pipehAngle = fGeo->PipehAngle(); //3.6 [degrees] half of the angular space needed for the target pipe
61  Double_t barWin_hthick = 0.1/2.; // [cm]=15um thickness of the 'glas' window at the readout end of the barbox;
62  Double_t phi0 = (180.-2.*pipehAngle)/bbnum + pipehAngle;
63  Double_t dphi = (180.-2.*pipehAngle)/bbnum*2.;
64  Double_t dphi_rad = dphi/180.*pi;
65 
66  Double_t bbox_zdown = fGeo->barBoxZDown(); // 130. bar box z downstream
67  Double_t bbox_zup = fGeo->barBoxZUp(); // -120. bar box z upstream
68  Double_t bbox_hlen = 0.5*(bbox_zdown - bbox_zup); // 125. bar box half length
69  Double_t bbox_shift = bbox_zup + bbox_hlen; // 5. bar box shift
70  Double_t bargap = fGeo->barGap(); // 0.01 half gap between bars
71  Double_t boxgap = fGeo->boxGap(); // 0.1 gap between bars and the bar box
72  Double_t boxthick = fGeo->boxThick(); // 0.05 thickness of the bar box
73  Double_t gluehthick = 0.0005;//[cm] // glue layer, connects two halfs into one long radiator bar
74  Double_t len = 0.; // length of the lenses block. see further
75  Double_t fSlabEnd = 0.; // [cm] position at which bar in front of EV ends
76  Double_t fdz_mirr1 = 0.;
77  Double_t fdz_mirr2 = 0.;
78  Double_t fdz_lens3 = 0.;
79  Double_t fdz_lens2 = 0.;
80  Double_t fdz_lens1 = 0.;
81 
82  //parameters for MCPs:
83  Double_t MCPsize = 5.78;//5.9; //[cm] width=height of an MCP
84  Double_t MCPactiveArea = 5.3; //[cm] width=height of the active area of an MCP
85  Double_t MCPgap = 0.1; //[cm] gap between MCPs
86  Double_t PixelSize = 0.65; //[cm] size of a pixel
87  Int_t Npix = 8; // number of mcps in a row/column
88  const Int_t NpixTotal = Npix*Npix; // total number of pixels for 1 mcp
89  Double_t PixelGap = (MCPactiveArea - Double_t(Npix)*PixelSize)/(Double_t(Npix)-1.);
90  Double_t PDwindowThick = 0.1; //[cm]
91  Double_t PhCathodeThick = 0.0001; //[cm] = 1 um
92  Double_t PDsensitiveThick = 0.01; //[cm]
93  Double_t PDgreaseLayer = 0.1; //[cm]
94  Double_t hgap = 0.5*(MCPsize - MCPactiveArea + MCPgap); // gap btw MCPs
95  Double_t step = MCPactiveArea + 2.*hgap; // step in which to locate MCPs = 6 for Mcp1, Mcp2; = 5.88 for Mcp2a. This is total size of one MCP with gaps in between
96 
97  // EV parameters
98  Double_t sob_len = fGeo->EVlen(); // 30. in current version
99  Double_t sob_shift = -bbox_hlen + bbox_shift - sob_len; // -150.
100  Double_t sob_angleB = fGeo->EVbackAngle(); //90. [degrees] angle of the EV (usually it is 90)
101  Double_t EVdrop = 0.5;//fGeo->EVdrop(); //0.7 [cm] drop of the EV - inner radius
102  Double_t EVoffset = fGeo->EVoffset(); //0. [cm] offset of the EV - outer radius
103  Double_t EVgreaseLayer = 0.0015; //[cm] grease layer thickness btw the bar window and the EV
104  //Double_t sob_angle = 45.;//fGeo->EVangle(); //80. [degrees] opening angle of the EV
105  Double_t sob_angle = atan((5.*step-2.*hthick-EVoffset-EVdrop)/sob_len)/pi*180.;// [degrees] opening angle of the EV
106  Double_t sob_Rout = (radius + hthick + EVoffset + sob_len*tan(sob_angle/180.*pi));// [cm] radius in the middle of the section
107 
108  Double_t bbAngle = ( 180. - 2.*pipehAngle - bbGap/radius/pi*180.*(bbnum/2.-1.) )/(bbnum/2.); // ~20 degrees
109  Double_t bbX = radius*bbAngle/180.*pi; // 17.45 cm
110 
111  Double_t radiusCornerIn = (radius-hthick-EVdrop)/cos(dphi/2./180.*pi);
112  Double_t radiusCornerOut = sob_Rout/cos(dphi_rad/2.);
113  Double_t radiusMiddleSmall = radius-hthick-EVdrop;
114 
115  cout<<"radius = "<<radius<<"rad corner in = "<< radiusCornerIn<<", rad corner out = "<< radiusCornerOut<<endl; cout<<"pixel gap = "<< PixelGap<<endl;
116  cout<<"sob_angle = "<<sob_angle<<endl;
117  cout<<"sob_Rout = "<<sob_Rout<<endl;
118  cout<<"bbAngle = "<<bbAngle<<", bbX = "<<bbX<<endl;
119  cout<<"dphi = "<<dphi<<", phi0 = "<<phi0<<endl;
120  //----------------------------------------------------------------------
121 
122  // prizm:
123  Double_t phlength = fGeo->PrismhLength(); //4.5 [cm] half length of the prizm
124  Double_t pdrop = fGeo->PrismDrop(); //0.5 [cm] drop of the prizm - inner side
125  Double_t pangle = fGeo->PrismAngle(); //30. [degrees] angle of the edge
126  Double_t poffset = fGeo->PrismOffset(); //1.[cm] prizm offset - outer side
127  Double_t pheight = 2.*phlength * tan(pangle/180.*pi); // [cm] half heigth of the prism
128 
129  Double_t sob_Rprizm = radius + hthick + poffset + pheight + EVoffset + (sob_len-2.*phlength)*tan(60./180.*pi);
130 
131  cout<<"prizm length = "<<phlength*2.<<endl;
132  cout<<"prizm height = "<<2.*phlength * tan(pangle/180.*pi)<<endl;
133  //----------------------------------------------------------
134 
135  // Rotations:
136  TGeoRotation rot1;
137  rot1.RotateZ(90.);
138 
139  TString fGeoFile= Form("../../geometry/dirc_l%d_p%d_mirrorGap_Mcp2a_GG_PC_correctedEVsensor.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  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
151  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
152  cout << "geoface setmediafile" << endl;
153  geoFace->setMediaFile("../../geometry/media_pnd.geo");
154  cout << "geoface readmedia" << endl;
155  geoFace->readMedia();
156  //cout << "geoface print" << endl;
157  //geoFace->print();
158  //cout << "geoface done" << endl;
159 
160  FairGeoMedia *Media = geoFace->getMedia();
161  FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();
162 
163  FairGeoMedium *DrcAir = Media->getMedium("air");
164  FairGeoMedium *DrcAirNoSens = Media->getMedium("DIRCairNoSens");
165  FairGeoMedium *DrcEpotek301_2 = Media->getMedium("Epotek301_2");
166  FairGeoMedium *DrcOpticalGrease = Media->getMedium("OpticalGrease");
167  FairGeoMedium *DrcCarbonFiber = Media->getMedium("DIRCcarbonFiber");
168  FairGeoMedium *DrcFusedSil = Media->getMedium("FusedSil");
169  FairGeoMedium *DrcMirror = Media->getMedium("Mirror");
170  FairGeoMedium *DrcMarcol82_7 = Media->getMedium("Marcol82-7");
171  FairGeoMedium *DrcNLAK33A = Media->getMedium("NLAK33A");
172  FairGeoMedium *DrcPhotocathode= Media->getMedium("Photocathode");
173 
174  Int_t nmed=geobuild->createMedium(DrcAir);
175  nmed=geobuild->createMedium(DrcAirNoSens);
176  nmed=geobuild->createMedium(DrcEpotek301_2);
177  nmed=geobuild->createMedium(DrcOpticalGrease);
178  nmed=geobuild->createMedium(DrcCarbonFiber);
179  nmed=geobuild->createMedium(DrcFusedSil);
180  nmed=geobuild->createMedium(DrcMirror);
181  nmed=geobuild->createMedium(DrcMarcol82_7);
182  nmed=geobuild->createMedium(DrcNLAK33A);
183  nmed=geobuild->createMedium(DrcPhotocathode);
184  //-------------------------------------------------------------------------------------------------------------
185  // focusing systems:
186  // 0 no ficusing:
187  if(fFocusingSystem == 0){
188  Double_t flen = 0.;
189  fSlabEnd = -bbox_hlen + bbox_shift + flen;
190  cout<<"bar ends at = "<<fSlabEnd<<endl;
191 
192  //fAtBarEnd = "DrcBar";
193  }
194 // 1 lens
195  if (fFocusingSystem == 1){ // L E N S E S (no airgaps, thin nlak33)
196  // some notes to lens operations (revision 9649) is in
197  // ~carsten/work/documents/software/PandaRoot/lens_definitions_2.pdf
198 
199  Double_t r = 75.18; // first lens radius (cm)
200  Double_t alpha = TMath::ASin(bbX/12./r); // lside -> bbX, assuming that there are 6 bars in the bar box??
201  Double_t a = r - r*TMath::Cos(alpha);
202  Double_t b = a + 0.5; // box dimension .6 instead of .5 due to strong curvature
203 
204  cout<<" DIRC a,b = "<<a<<" "<<b<<endl;
205 
206 
207  Double_t r2 = 75.18; // radius second lens (cm)
208  Double_t b2 = 0.2; // + a2;
209 
210  Double_t r3 = 17.95; // third lens radius (cm)
211  Double_t alpha3 = TMath::ASin(bbX/12./r3); // lside -> bbX, assuming that there are 6 bars in the bar box??
212  Double_t a3 = r3 - r3*TMath::Cos(alpha3);
213  Double_t b3 = a3 + 0.0; // box dimension .6 instead of .5 due to strong curvature
214 
215  // lens1 (b) + lens2 (0.6) + lens3 (b3) + gap (0.5)
216 
217  len = b + 0.2 + b3 + 0.5;//+ a2; // dimension of the box containing both lenses
218 
219  cout<<"DIRC len= "<<len<<endl;
220 
221  fSlabEnd = -bbox_hlen + bbox_shift + len; // used in processHits
222  cout<<"bar ends at = "<<fSlabEnd<<endl;
223 
224  // Lenses
225 
226  // Lens 1
227  Double_t t = -r +b/2;
228 
229  TGeoSphere* logicSphere= new TGeoSphere("S",0.,r, 0. ,180.,0.,360.);
230  TGeoBBox* lBox = new TGeoBBox("B", (bbX/barnum)/2-bargap, hthick, b/2.);
231  TGeoTranslation *tr1 = new TGeoTranslation("tr1", 0.,0., t);
232  tr1->RegisterYourself();
233  TGeoCompositeShape *cs = new TGeoCompositeShape("cs","S*(B:tr1)");
234  TGeoVolume *lens1 = new TGeoVolume("DrcLENS1Sensor",cs, gGeoManager->GetMedium("FusedSil"));
235  lens1->SetLineColor(kRed-8);
236  lens1->SetTransparency(40);
237 
238  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r + b
239  // with -0.01 one can make a gap visible (.1mm) for orientation
240  //barContainer->(lens1, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r+b) /*-0.01*/, new TGeoRotation (0)));
241 /* barContainer->AddNode(lens1, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r+b), new TGeoRotation (0)));
242 */ fdz_lens1 = -(bbox_hlen)+len -(-r+b);
243 
244 
245  //Lens 2
246  Double_t t2 = -r2;// +b2/2 r2 is the reference point (concave lens)
247  TGeoSphere* logicSphere2 = new TGeoSphere("S2",0 ,r2, 0. ,180.,0.,360.);
248  TGeoBBox* lBox2 = new TGeoBBox("B2", (bbX/barnum)/2.-bargap, hthick, b2/2.); // lside -> bbX
249  TGeoTranslation *tr2 = new TGeoTranslation("tr2", 0.,0., t2);
250  tr2->RegisterYourself();
251  TGeoCompositeShape *cs2 = new TGeoCompositeShape("cs2","(B2:tr2)-S2");
252  TGeoVolume *lens2 = new TGeoVolume("DrcLENS2Sensor",cs2, gGeoManager->GetMedium("NLAK33A"));
253  lens2->SetLineColor(kRed+2);
254  lens2->SetTransparency(40);
255 
256  // place the lens exactly on lens1
257  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r2
258  // the tip of lens1 is at b
259  // with -0.02 one can make a gap visible (.1mm due to lens1) for orientation
260  //barContainer->AddNode(lens2, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r2) -b /*-0.02*/, new TGeoRotation (0)));
261 /* barContainer->AddNode(lens2, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r2) -b , new TGeoRotation (0)));
262 */ fdz_lens2 = -(bbox_hlen)+len -(-r2) -b;
263 
264  //Lens3 (like lens1, same treatment)
265  Double_t t3 = -r3+b3/2;
266  TGeoSphere* logicSphere3= new TGeoSphere("S3",0.,r3, 0. ,180.,0.,360.);
267  TGeoBBox* lBox3 = new TGeoBBox("B3", (bbX/barnum)/2-bargap, hthick, b3/2.); // lside -> bbX
268  TGeoTranslation *tr3 = new TGeoTranslation("tr3", 0.,0., t3);
269  tr3->RegisterYourself();
270  TGeoCompositeShape *cs3 = new TGeoCompositeShape("cs3","S3*(B3:tr3)");
271  TGeoVolume *lens3 = new TGeoVolume("DrcLENS3Sensor",cs3, gGeoManager->GetMedium("NLAK33A"));
272  lens3->SetLineColor(kRed-6);
273  lens3->SetTransparency(40);
274 
275  // place the lens exactly on lens2 plane side
276  // position lens within already shifted bar container at -(bbox_hlen-eps)+len, the lens base is -r3 + b3
277  // with -0.03 one can make a gap visible (.1mm due to lens 1&2) for orientation
278  // b2/2 is the thickness of lens2 in the middle
279  //barContainer->AddNode(lens3, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen-eps)+len -(-r3+b3) -b - b2/2 /*-0.03*/, new TGeoRotation (0)));
280 /* barContainer->AddNode(lens3, 1,new TGeoCombiTrans(0., 0., -(bbox_hlen)+len -(-r3+b3) -b - b2/2 , new TGeoRotation (0)));
281  //cout<<" DIRC r,r2,r3 = "<<r<<" "<<r2<<" "<<r3<<endl;
282 */ fdz_lens3 = -(bbox_hlen)+len -(-r3+b3) -b - b2/2. ;
283 
284  //fAtBarEnd = "DrcLENS3";
285  } // E N D O F N E W L E N S E S
286 
287  // 2 - mirror:
288  if (fFocusingSystem == 2){ // Mirrors at front
289  // Put some mirrors at the downstream end with a focal plane at the PD.
290 
291  double zpos=120.;
292 
293  // The bar is produced with "bbox_hlen-fabs(len)/2.-mirr_hthick"
294  len = -130. + zpos + 1.;
295  //len = -247.0; // negative to make space at downstream end of bar.
296 
297  Double_t len1 = 1.5; // 1st block
298 
299  Double_t mirror_angle = 0.0;// 0.0 is pointing upstream, -90.0 is pointing to the beam axis
300  Double_t focal_length = 2. * bbox_hlen + sob_len - fabs(len) + len1;
301  Double_t mirror_radius = 2. * focal_length;
302 
303  cout<<" mirror radius: "<<mirror_radius<<endl;
304 
305 
306  // no angle for the moment
307  // block
308  TGeoSphere* logicSphere = new TGeoSphere("S",0.,mirror_radius, 0. ,180.,0.,360.);
309  TGeoBBox* lBox = new TGeoBBox("B", (bbX/barnum)/2-bargap, hthick, fabs(len1)/2.); // lside -> bbX
310 
311  Double_t t = mirror_radius - len1/2.;
312 
313  TGeoTranslation *tr1 = new TGeoTranslation("tr1", 0.,0., t);
314  tr1->RegisterYourself();
315  TGeoCompositeShape *cs = new TGeoCompositeShape("cs","S*(B:tr1)");
316 
317  TGeoVolume *block1 = new TGeoVolume("DrcBlock1",cs, gGeoManager->GetMedium("FusedSil"));
318  block1->SetLineColor(kRed);
319  block1->SetTransparency(40);
320 
321 
322  Double_t shift1 = len1-mirror_radius; // Now the start of block1 is at zero
323  shift1 += bbox_hlen-fabs(len)-2.*mirr_hthick;
324 
325 /*
326  barContainer->AddNode(block1, 1, new TGeoCombiTrans(0., 0., shift1, new TGeoRotation (0)));
327 */ fdz_mirr1 = shift1;
328 
329 
330  Double_t gap = 0;
331  Double_t len2 = fabs(len) - len1 - gap;
332 
333  // no angle for the moment
334  // block
335  TGeoSphere* logicSphere2 = new TGeoSphere("S2",0,mirror_radius, 0. ,180.,0.,360.);
336  TGeoBBox* lBox2 = new TGeoBBox("B2", (bbX/barnum)/2-bargap, hthick, fabs(len2)/2.); // lside -> bbX
337 
338  // make radius part of block
339  Double_t t2 = mirror_radius + len2/2 - 1;
340 
341  TGeoTranslation *tr2 = new TGeoTranslation("tr2", 0.,0., t2);
342  tr2->RegisterYourself();
343  TGeoCompositeShape *cs2 = new TGeoCompositeShape("cs2","(B2:tr2)-S2");
344 
345  TGeoVolume *block2 = new TGeoVolume("DrcBlock2",cs2, gGeoManager->GetMedium("Mirror"));
346  block2->SetLineColor(kGreen);
347  block2->SetTransparency(40);
348 
349 
350  Double_t shift2 = -mirror_radius; // Now the start of the block is at zero
351  shift2 += bbox_hlen-fabs(len)-2*mirr_hthick; // at end of bar
352  shift2 += len1 + gap;
353 
354  fdz_mirr2 = shift2;
355 
356 /*
357  barContainer->AddNode(block2,
358  1,
359  new TGeoCombiTrans(0.,
360  0.,
361  shift2,
362  // bbox_hlen-mirror_radius-2*mirr_hthick-len2+len1+1,
363  //bbox_hlen-mirror_radius-len2+0.1,
364  new TGeoRotation (0)
365  )
366  );
367 */
368 
369 
370  Double_t flen = 0.;
371  fSlabEnd = -bbox_hlen + bbox_shift + flen;
372  cout<<"bar ends at = "<<fSlabEnd<<endl;
373 
374  //fAtBarEnd = "DrcBar";
375  } // E N D O F MIRRORS
376 
377  if(fFocusingSystem == 3){ // cylindrical lens w/o airgap
378 
379  // main parameters:
380 
381  Double_t Hcyl2 = 0.1; //[cm] thickness in the middle of the second (FusedSil) lens
382  Double_t Rcyl = 7.35; // [cm]
383  Double_t Lcyl = (bbX/barnum)/2.-bargap; //cylinder length/2.
384  Double_t Acyl = TMath::ASin(hthick/Rcyl); // angle
385  Double_t Tcyl1 = 0.5*Rcyl*(1. + TMath::Cos(Acyl)); // distance between centers of the cylinder and the box
386 
387  len = Hcyl2 + Rcyl*(1. - TMath::Cos(Acyl));
388  cout<<"DIRC: len = "<<len<<endl;
389 
390  cout<<"length = "<< Lcyl<<", angle = "<<Acyl/3.1415*180.<<", translation = "<<Tcyl1<<endl;
391  cout<<"lens 0.5*width = "<<hthick<<", 0.5*thickness = "<< 0.5*Rcyl*(1.-TMath::Cos(Acyl))<<", 0.5*heigth = "<<Lcyl <<endl;
392 
393  //Lens1
394  TGeoEltu* lCylinder = new TGeoEltu("Cyl",Rcyl,Rcyl,Lcyl);
395  TGeoBBox* lCylBox = new TGeoBBox("CylBox", hthick, hthick, Lcyl);
396  //TGeoBBox* lCylBox = new TGeoBBox("CylBox",36.75, 36.75, Lcyl);
397  //TGeoTranslation *trCyl = new TGeoTranslation("trCyl", 0., 100.4, 0.);
398  TGeoTranslation *trCyl1 = new TGeoTranslation("trCyl", 0., Rcyl*TMath::Cos(Acyl)+hthick, 0.);// !!!
399  trCyl1->RegisterYourself();
400  TGeoCompositeShape *llens1 = new TGeoCompositeShape("llens1","Cyl*(CylBox:trCyl)");
401  TGeoVolume *CylLens1 = new TGeoVolume("DrcLENS1Sensor",llens1, gGeoManager->GetMedium("NLAK33A"));
402  CylLens1->SetLineColor(kRed-8);
403  CylLens1->SetTransparency(40);
404 
405  fdz_lens1 = -(bbox_hlen-barWin_hthick) + len - (-Rcyl + Rcyl*(1.-TMath::Cos(Acyl)));
406 
407  //Lens2
408  TGeoTranslation *trCyl2 = new TGeoTranslation("trCyl2", 0., Rcyl + Hcyl2 - hthick, 0.);
409  trCyl2->RegisterYourself();
410  TGeoCompositeShape *llens2 = new TGeoCompositeShape("llens2", "(CylBox:trCyl2) - Cyl");
411  TGeoVolume* CylLens2 = new TGeoVolume("DrcLENS2Sensor", llens2, gGeoManager->GetMedium("FusedSil"));
412  CylLens2->SetLineColor(kRed+2);
413  CylLens2->SetTransparency(40);
414 
415  // TGeoRotation rot_lens3;
416  // rot_lens3.RotateZ(-90.);
417  // rot_lens3.RotateX( 90.);
418 
419 
420 
421  fdz_lens2 = -(bbox_hlen-barWin_hthick) + len + (hthick - Rcyl*(1.-TMath::Cos(Acyl)) - Hcyl2);
422  fSlabEnd = -bbox_hlen + bbox_shift;
423  cout<<"bar ends at = "<<fSlabEnd<<endl;
424 
425  }
426  //---------------------------------------------------------------------------------------------
427 
428 
429  // create top volume:
430  TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
431 
432  TGeoBBox* lTop = new TGeoBBox(500,500,300);
433  TGeoVolume* top = new TGeoVolume("DIRC", lTop, gGeoManager->GetMedium("air"));
434  gGeoManager->SetTopVolume(top);
435 
436  // create pre-top volume:
437  TGeoVolume* vLocalMother;
438  TGeoPcon* shape = new TGeoPcon("BarrelDIRCShape", 0, 360., 4);
439  shape->DefineSection(0, bbox_zdown, 35., 60.);
440  shape->DefineSection(1, bbox_zup, 35., 60.);
441  shape->DefineSection(2, bbox_zup - sob_len, 35., sob_Rout+poffset+pheight+EVoffset+1.);
442  shape->DefineSection(3, bbox_zup - sob_len - PDbaseLayer - EVgreaseLayer, 35., sob_Rout+poffset+pheight+EVoffset+1.);
443  vLocalMother = new TGeoVolume("BarrelDIRC", shape, gGeoManager->GetMedium("DIRCairNoSens")); //DIRCairNoSens_m); //################
444  top->AddNode(vLocalMother, 0,0);
445 
446  cout<<"bbox length = "<<2.*(bbox_hlen-0.5*(boxgap+boxthick))<<", prizm length = "<<2.*(phlength+0.5*(boxgap+boxthick))<<endl;
447  cout<<"prizm shift = "<<-bbox_hlen+0.5*(boxgap+boxthick)-phlength<<endl;
448 
449  // create BarBoxes:
450  TGeoBBox* logicbbL;
451  TGeoVolume *bbox;
452  logicbbL = new TGeoBBox("logicbbL", bbX/2.+boxthick, hthick+boxgap+boxthick, bbox_hlen);
453  bbox = new TGeoVolume("DrcBarBox", logicbbL,gGeoManager->GetMedium("DIRCcarbonFiber"));
454  bbox->SetLineColor(30);
455 
456 
457  TGeoBBox* logicbbS;
458  TGeoVolume *abox;
459  if(fFocusingSystem == 0 || fFocusingSystem == 1 || fFocusingSystem == 2){
460  logicbbS = new TGeoBBox("logicbbS", bbX/2., hthick+boxgap, bbox_hlen-barWin_hthick);
461  }
462  if(fFocusingSystem == 3){
463  logicbbS = new TGeoBBox("logicbbS", bbX/2., hthick+boxgap, bbox_hlen-barWin_hthick/*-len/2.*/);
464  }
465  abox = new TGeoVolume("DrcAirBox", logicbbS, gGeoManager->GetMedium("DIRCairNoSens"));
466  if(fFocusingSystem == 0 || fFocusingSystem == 1 || fFocusingSystem == 2){
467  bbox->AddNode(abox, 1, new TGeoCombiTrans(0., 0., barWin_hthick, new TGeoRotation(0)));
468  }
469  if(fFocusingSystem == 3){
470  bbox->AddNode(abox, 1, new TGeoCombiTrans(0., 0., barWin_hthick/*+len/2.*/, new TGeoRotation(0)));
471  }
472  abox->SetLineColor(19); // gray
473 
474  //create windows at the readout end of the bar boxes:
475  TGeoBBox* logicBarWin = new TGeoBBox("logicBarWin", bbX/2., hthick+boxgap, barWin_hthick);
476  TGeoVolume* barwin = new TGeoVolume("DrcBarboxWindowSensor", logicBarWin, gGeoManager->GetMedium("FusedSil"));
477  barwin->SetLineColor(kBlue-4);
478  bbox->AddNode(barwin, 1, new TGeoCombiTrans(0.,0.,-bbox_hlen+barWin_hthick,new TGeoRotation(0)));
479  /*
480  if(fFocusingSystem == 3){
481  TGeoRotation rot_lens3;
482  rot_lens3.RotateZ(-90.);
483  rot_lens3.RotateY( 90.);
484  bbox->AddNode(CylLens1, 1, new TGeoCombiTrans(0., 0., fdz_lens1+2.*barWin_hthick, new TGeoRotation(rot_lens3)));
485  bbox->AddNode(CylLens2, 1, new TGeoCombiTrans(0., 0., fdz_lens1+2.*barWin_hthick, new TGeoRotation(rot_lens3)));
486  }
487  */
488  //create layers of grease at the readout end of the bar boxes, between the windows and the EV:
489  TGeoBBox* logicEVgrease = new TGeoBBox("logicEVgrease", bbX/2., hthick+boxgap, EVgreaseLayer/2.);
490  TGeoVolume* evgrease = new TGeoVolume("DrcEVgrease", logicEVgrease, gGeoManager->GetMedium("OpticalGrease"));
491  evgrease->SetLineColor(kSpring);
492 
493  // put barboxes into right positions:
494  Double_t dx_bbox, dy_bbox, dz_bbox, phi_curr;
495 
496  for(Int_t m = 0; m < bbnum; m ++){
497  phi_curr = (90. - phi0 - dphi*m)/180.*pi;
498  if(m > bbnum/2-1){ phi_curr = (90. - phi0 - dphi*m - 2.*pipehAngle)/180.*pi; }
499  dx_bbox = radius * cos(phi_curr);
500  dy_bbox = radius * sin(phi_curr);
501  dz_bbox = bbox_shift;
502 
503  TGeoRotation rot_bbox;
504  rot_bbox.RotateZ( -phi0 - m*dphi - (TMath::Floor(2.*m/bbnum))*(2.*pipehAngle));
505  vLocalMother->AddNode(bbox, m+1, new TGeoCombiTrans(dx_bbox, dy_bbox, dz_bbox, new TGeoRotation(rot_bbox)));
506  vLocalMother->AddNode(evgrease, m+1, new TGeoCombiTrans(dx_bbox, dy_bbox, bbox_zup-EVgreaseLayer/2.,new TGeoRotation(rot_bbox)));
507  }
508 
509  cout<<"bar width = "<<2.*(((bbX/barnum)/2.)-bargap) << ", bar with gaps = "<<(bbX/barnum)<<endl;
510  cout<<"barboxL width = "<<bbX/2.+boxgap+boxthick<<", barboxS width = "<<bbX/2.+boxgap<<endl;
511 
512  // create logic bar:
513  TGeoBBox* logicBar;
514  if(fFocusingSystem == 0 || fFocusingSystem == 1 || fFocusingSystem == 2){
515  logicBar = new TGeoBBox("logicBar", ((bbX/barnum)/2.)-bargap, hthick, bbox_hlen-fabs(len)/2.-mirr_hthick-barWin_hthick);
516  }
517  if(fFocusingSystem == 3){
518  logicBar = new TGeoBBox("logicBar", ((bbX/barnum)/2.)-bargap, hthick, bbox_hlen-fabs(len)/2.-mirr_hthick-barWin_hthick);
519  }
520 
521  TGeoVolume* bar;
522  bar = new TGeoVolume("DrcBarSensor",logicBar, gGeoManager->GetMedium("FusedSil"));
523  bar->SetLineColor(kCyan-9);
524  bar->SetTransparency(60);
525 
526  // create logic mirror:
527  TGeoBBox* logicMirror = new TGeoBBox("logicMirror", bbX/barnum/2.-bargap, hthick, mirr_hthick /2.); //$$$$$$$
528  TGeoVolume *mirr = new TGeoVolume("DrcMirr", logicMirror, gGeoManager->GetMedium("Mirror"));
529  mirr->SetLineColor(5);
530 
531  // create glue layer inside the bar (connects two halves):
532  TGeoBBox* logicBarGlue = new TGeoBBox("logicBarGlue", bbX/barnum/2.-bargap, hthick, gluehthick); //$$$$$$$
533  TGeoVolume *barglue = new TGeoVolume("DrcBarGlueSensor", logicBarGlue, gGeoManager->GetMedium("Epotek301_2"));
534  barglue->SetLineColor(kSpring-5);
535  bar->AddNode(barglue, 1, new TGeoCombiTrans(0., 0., 0., new TGeoRotation (0)));
536 
537  Double_t dx, dy, dz_bar, dz_mirr;
538 
539  for(Int_t j=0; j<barnum; j++){
540  dx = - (bbX/2.) + (bbX/barnum)/2. + j * (bbX/barnum);
541  dy = 0.;//len/2. - mirr_hthick;
542  if(fFocusingSystem == 0 || fFocusingSystem == 1 || fFocusingSystem == 2){
543  dz_mirr = bbox_hlen - barWin_hthick - mirr_hthick;
544  dz_bar = len/2.-mirr_hthick;
545  }
546  if(fFocusingSystem == 3){
547  dz_mirr = bbox_hlen - barWin_hthick -fabs(len)/2. - mirr_hthick;
548  dz_bar = -mirr_hthick;
549  }
550  if(fFocusingSystem == 1){ // lens
551  abox->AddNode(lens1, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation (0)));
552  abox->AddNode(lens2, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens2, new TGeoRotation (0)));
553  abox->AddNode(lens3, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens3, new TGeoRotation (0)));
554  }
555  if(fFocusingSystem == 2){ // mirror
556  abox->AddNode(block1, 1+j, new TGeoCombiTrans(dx, dy, fdz_mirr1, new TGeoRotation (0)));
557  abox->AddNode(block2, 1+j, new TGeoCombiTrans(dx, dy, fdz_mirr2, new TGeoRotation (0)));
558  }
559  if(fFocusingSystem == 3){
560  TGeoRotation rot_lens3;
561  rot_lens3.RotateZ(-90.);
562  rot_lens3.RotateY( 90.);
563  abox->AddNode(CylLens1, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation(rot_lens3)));
564  abox->AddNode(CylLens2, 1+j, new TGeoCombiTrans(dx, dy, fdz_lens1, new TGeoRotation(rot_lens3)));
565  }
566  abox->AddNode(bar, 1+j, new TGeoCombiTrans(dx, dy, dz_bar, new TGeoRotation(0)));
567  if(fFocusingSystem != 2){ // not a forward mirror
568  abox->AddNode(mirr, 1+j, new TGeoCombiTrans(dx, dy, dz_mirr, new TGeoRotation(0)));
569  }
570  }
571 
572  // Expansion volume:
573  Double_t dR;
574  Double_t xEV;
575  Double_t cosFactor1;
576 
577  TGeoPgon* logicEV1, * logicEV2, *logicEV3, * logicEV4;
578  cosFactor1 = cos(pipehAngle/180.*pi)/cos(dphi/180.*pi/2.);
579  dR = (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi) - (radius-hthick);
580  xEV = (dR + sob_len*tan(sob_angle/180.*pi))/ (tan(sob_angle/180.*pi) + tan(sob_angleB/180.*pi));
581  if(sob_angleB == 90.){
582  logicEV1 = new TGeoPgon("logicEV1", 93.6, 172.8, bbnum/2, 2);
583  logicEV1->DefineSection(0, 0., radiusMiddleSmall, sob_Rout);
584  logicEV1->DefineSection(1, sob_len, radiusMiddleSmall, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi));
585  logicEV2 = new TGeoPgon("logicEV2", -86.4, 172.8, bbnum/2, 2);
586  logicEV2->DefineSection(0, 0., radiusMiddleSmall, sob_Rout);
587  logicEV2->DefineSection(1, sob_len, radiusMiddleSmall, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi));
588  logicEV3 = new TGeoPgon("logicEV3", 86.4, 7.2, 1, 2);
589  logicEV3->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1, sob_Rout*cosFactor1);
590  logicEV3->DefineSection(1, sob_len, (radiusMiddleSmall)*cosFactor1, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi)*cosFactor1);
591  logicEV4 = new TGeoPgon("logicEV4", -93.6, 7.2, 1, 2);
592  logicEV4->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1, sob_Rout*cosFactor1);
593  logicEV4->DefineSection(1, sob_len, (radiusMiddleSmall)*cosFactor1, (radius+hthick+boxgap+boxthick+EVoffset)/cos(dphi/2./180.*pi)*cosFactor1);
594 
595  }
596  if(sob_angleB != 90.){
597  logicEV1 = new TGeoPgon("logicEV1", 93.6, 172.8, bbnum/2, 3);
598  logicEV1->DefineSection(0, 0., radius-hthick, radius-hthick+eps);
599  logicEV1->DefineSection(1, xEV, radius-hthick, sob_Rout - xEV*tan(sob_angle/180.*pi));
600  logicEV1->DefineSection(2, sob_len, radius-hthick, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi));
601  logicEV2 = new TGeoPgon("logicEV2", -86.4, 172.8, bbnum/2, 3);
602  logicEV2->DefineSection(0, 0., radius-hthick, radius-hthick+eps);
603  logicEV2->DefineSection(1, xEV, radius-hthick, sob_Rout - xEV*tan(sob_angle/180.*pi));
604  logicEV2->DefineSection(2, sob_len, radius-hthick, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi));
605  logicEV3 = new TGeoPgon("logicEV3", 86.4, 7.2, 1, 3);
606  logicEV3->DefineSection(0, 0., (radius-hthick)*cosFactor1, (radius-hthick+eps)*cosFactor1);
607  logicEV3->DefineSection(1, xEV, (radius-hthick)*cosFactor1, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1);
608  logicEV3->DefineSection(2, sob_len, (radius-hthick)*cosFactor1, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi)*cosFactor1);
609  logicEV4 = new TGeoPgon("logicEV4", -93.6, 7.2, 1, 3);
610  logicEV4->DefineSection(0, 0., (radius-hthick)*cosFactor1, (radius-hthick+eps)*cosFactor1);
611  logicEV4->DefineSection(1, xEV, (radius-hthick)*cosFactor1, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1);
612  logicEV4->DefineSection(2, sob_len, (radius-hthick)*cosFactor1, (radius+hthick+boxgap+boxthick)/cos(dphi/2./180.*pi)*cosFactor1);
613  }
614 
615  TGeoCompositeShape *logicEV = new TGeoCompositeShape("logicEV","logicEV1 + logicEV2 + logicEV3 + logicEV4");
616  TGeoVolume* baseEV = new TGeoVolume("DrcEVSensor", logicEV, gGeoManager->GetMedium("FusedSil"));//Marcol82_7_m); //####################################
617  baseEV->SetLineColor(kMagenta+2);
618  baseEV->SetTransparency(50);
619  vLocalMother->AddNode(baseEV, 1, new TGeoCombiTrans(0.,0.,sob_shift - EVgreaseLayer, new TGeoRotation(0)));
620 
621  // 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:
622  TGeoPgon *logicPDbase1 = new TGeoPgon("logicPDbase1", 93.6, 172.8, bbnum/2, 2);
623  TGeoPgon *logicPDbase2 = new TGeoPgon("logicPDbase2", -86.4, 172.8, bbnum/2, 2);
624  TGeoPgon *logicPDbase3 = new TGeoPgon("logicPDbase3", 86.4, 7.2, 1, 2);
625  TGeoPgon *logicPDbase4 = new TGeoPgon("logicPDbase4", -93.6, 7.2, 1, 2);
626  Double_t rad_delta = (MCPsize-MCPactiveArea)/2./cos(45./180.*pi);
627  if(sob_angleB == 90.){
628  logicPDbase1->DefineSection(0, 0., radiusMiddleSmall-rad_delta, sob_Rout);
629  logicPDbase1->DefineSection(1, PDbaseLayer, radiusMiddleSmall-rad_delta, sob_Rout);
630  logicPDbase2->DefineSection(0, 0., radiusMiddleSmall-rad_delta, sob_Rout);
631  logicPDbase2->DefineSection(1, PDbaseLayer, radiusMiddleSmall-rad_delta, sob_Rout);
632  logicPDbase3->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
633  logicPDbase3->DefineSection(1, PDbaseLayer, (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
634  logicPDbase4->DefineSection(0, 0., (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
635  logicPDbase4->DefineSection(1, PDbaseLayer, (radiusMiddleSmall)*cosFactor1-rad_delta, sob_Rout*cosFactor1);
636  }
637 /* if(sob_angleB != 90.){
638  logicPD1->DefineSection(0, 0., radius-hthick+eps, radius-hthick+eps+0.1);
639  logicPD1->DefineSection(1, xEV, sob_Rout - xEV*tan(sob_angle/180.*pi), sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1);
640  logicPD2->DefineSection(0, 0., radius-hthick+eps, radius-hthick+eps+0.1);
641  logicPD2->DefineSection(1, xEV, sob_Rout - xEV*tan(sob_angle/180.*pi), sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1);
642  logicPD3->DefineSection(0, 0., (radius-hthick+eps)*cosFactor1, (radius-hthick+eps+0.1)*cosFactor1);
643  logicPD3->DefineSection(1, xEV, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1, (sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1)*cosFactor1);
644  logicPD4->DefineSection(0, 0., (radius-hthick+eps)*cosFactor1, (radius-hthick+eps+0.1)*cosFactor1);
645  logicPD4->DefineSection(1, xEV, (sob_Rout - xEV*tan(sob_angle/180.*pi))*cosFactor1, (sob_Rout- xEV*tan(sob_angle/180.*pi)+0.1)*cosFactor1);
646  }*/
647 
648  TGeoCompositeShape *logicPDbase = new TGeoCompositeShape("logicPDbase","logicPDbase1 + logicPDbase2 + logicPDbase3 + logicPDbase4");
649  TGeoVolume *pdbase = new TGeoVolume("DrcPDbase", logicPDbase, gGeoManager->GetMedium("DIRCcarbonFiber"));
650  pdbase->SetLineColor(kGreen-6);
651  if(sob_angleB == 90.){
652  vLocalMother->AddNode(pdbase, 1, new TGeoCombiTrans(0., 0., sob_shift-EVgreaseLayer-PDbaseLayer, new TGeoRotation(0)));
653  }
654 /* if(sob_angleB != 90.){
655  vLocalMother->AddNode(pd, 1, new TGeoCombiTrans(0., 0., sob_shift, new TGeoRotation(0)));
656  }
657 */
658 
659 /* //MCPs as PD
660  //++++++++++++++++++++++++++++++
661  // create individual pixels:
662  TGeoBBox* logicPixel = new TGeoBBox("logicPixel", PixelSize/2., PixelSize/2., PDsensitiveThick/2.);
663  TGeoVolume* pix = new TGeoVolume("DrcPDSensor",logicPixel, gGeoManager->GetMedium("FusedSil"));
664  pix->SetLineColor(kGreen+1);
665 
666  // create pixel holders:
667  TGeoBBox* logicPixelHolder = new TGeoBBox("logicPixelHolder", MCPactiveArea/2., MCPactiveArea/2., PDsensitiveThick/2.);
668  TGeoVolume *pixelholder = new TGeoVolume("DrcPixelHolder", logicPixelHolder, gGeoManager->GetMedium("DIRCcarbonFiber"));
669  pixelholder->SetLineColor(kOrange-8);
670 
671  //put pixels into holders:
672  for(Int_t itr=0; itr<Npix; itr++){
673  Double_t dy_row = (-Npix/2 + itr)*(PixelSize + PixelGap) + 0.5*(PixelSize+PixelGap);
674  for(Int_t jtr=0; jtr<Npix; jtr++){
675  Double_t dx_row = (-Npix/2 + jtr)*(PixelSize + PixelGap) + 0.5*(PixelSize+PixelGap);
676  pixelholder->AddNode(pix, itr*Npix+jtr, new TGeoCombiTrans(dx_row, dy_row, 0., new TGeoRotation(0)));
677  }
678  }
679  //++++++++++++++++++++++++++++++
680 */
681  //create pixel plates: one for each MCP:
682  TGeoBBox* logicPD = new TGeoBBox("logicPD", MCPactiveArea/2., MCPactiveArea/2., PDsensitiveThick/2.);
683  TGeoVolume *pixelholder = new TGeoVolume("DrcPDSensor", logicPD, gGeoManager->GetMedium("FusedSil"));
684  pixelholder->SetLineColor(kGreen+1);
685 
686  // create photo cathodes for each MCP
687  TGeoBBox* logicPhCathode = new TGeoBBox("logicPhCathode", MCPactiveArea/2., MCPactiveArea/2., PhCathodeThick/2.);
688  TGeoVolume *phcathode = new TGeoVolume("DrcPhCathodeSensor", logicPhCathode, gGeoManager->GetMedium("Photocathode"));//("FusedSil"));
689  phcathode->SetLineColor(kGray+1);
690 
691  // create Windows for each MCP
692  TGeoBBox* logicWindow = new TGeoBBox("logicWindow", MCPsize/2., MCPsize/2., PDwindowThick/2.);
693  TGeoVolume *window = new TGeoVolume("DrcPDwindowSensor", logicWindow, gGeoManager->GetMedium("FusedSil"));
694  window->SetLineColor(kBlue-4);
695 
696  // create grease layers between MCP window and the EV back side
697  TGeoBBox* logicMCPgrease = new TGeoBBox("logicMCPgrease", MCPsize/2., MCPsize/2., PDgreaseLayer/2.);
698  TGeoVolume *mcpgrease = new TGeoVolume("DrcMcpGreaseSensor", logicMCPgrease, gGeoManager->GetMedium("OpticalGrease"));
699  mcpgrease->SetLineColor(kSpring);
700 
701 /* //version 0: full PD as a single layer:
702  TGeoCompositeShape *logicPD = new TGeoCompositeShape("logicPD","logicPDbase1 + logicPDbase2 + logicPDbase3 + logicPDbase4");
703  TGeoVolume *pd = new TGeoVolume("DrcPDSensor", logicPD, gGeoManager->GetMedium("FusedSil"));
704  pdbase->SetLineColor(kGreen-6);
705  if(sob_angleB == 90.){
706  vLocalMother->AddNode(pd, 1, new TGeoCombiTrans(0., 0., sob_shift-EVgreaseLayer-PDbaseLayer, new TGeoRotation(0)));
707  }
708 */
709 
710  //version 1: place MCPs in global cs
711 
712  cout<<"STEP = "<<step<<endl;
713 /* Double_t xpos = 0.;
714  Double_t ypos = 0.;
715  Double_t xcurr = 0.;
716  Double_t ycurr = 0.;
717  Double_t xsmall, xlarge, lastxsmall;
718  Int_t i=0;
719  cout<<"hgap = "<<hgap<<", step = "<<step<<endl;
720  xpos = radiusCornerIn + 0.5*MCPactiveArea;
721  ypos = step/2.;
722  int cas1 = 1;
723  int cas2 = 1;
724  Int_t TotalNmcp = 0;
725  for(Int_t nquat =0; nquat<4; nquat++){
726  for(Int_t i=0; i<Int_t((radiusCornerOut)/step); i++){ // loop over y
727  cas1 = findSectorIn(ypos-step/2., dphi_rad, radius, EVdrop, hthick);
728  cout<<"+++++++++++++++ nquat = "<<nquat<<", i = "<<i<<", xpos = "<<xpos<<", ypos = "<<ypos<<endl;
729  cout << "calculated sector1 from y =" << ypos << " is " << cas1 << endl;
730  getXrange(dphi, radiusCornerOut+1.3, radiusCornerIn, MCPactiveArea, cas1, cas2, ypos-(step/2.), step, xsmall, xlarge);
731  xpos = xsmall;
732  if(xsmall < xlarge){
733  for(Int_t j=0; j<TMath::Ceil((xlarge-xsmall)/step); j++){ // loop over x
734  cout<<"j = "<<j<<". quat = "<<nquat<<", xsmall = "<<xsmall<<", xlarge = "<<xlarge<<", cas1 = "<<cas1<<", limit = "<<(Int_t((xlarge-xsmall)/step))+1<<endl;
735 
736  if(nquat == 1 || nquat == 2){
737  xcurr = -xpos;
738  }
739  if(nquat == 0 || nquat == 3){
740  xcurr = xpos;
741  }
742  if(nquat == 0 || nquat == 1){
743  ycurr = ypos;
744  }
745  if(nquat == 2 || nquat == 3){
746  ycurr = -ypos;
747  }
748  // put Photo Cathodes
749  pdbase->AddNode(phcathode, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(0)));
750  //put windows
751  pdbase->AddNode(window, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(0)));
752  //put grease layers
753  pdbase->AddNode(mcpgrease, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(0)));
754  //put pixel holders:
755  pdbase->AddNode(pixelholder, TotalNmcp, new TGeoCombiTrans(xcurr, ycurr, PDbaseLayer-PDgreaseLayer-PDwindowThick-PhCathodeThick-PDsensitiveThick/2., new TGeoRotation(0)));
756  TotalNmcp = TotalNmcp + 1;
757  xpos = xpos + step;
758  //if(xpos >= xlarge){continue;}
759  }
760  cas2 = findSectorOut(ypos+step/2., dphi_rad, radiusCornerOut);
761  cout<<"calculated sector2 from y = "<<ypos+step/2. << " is "<< cas2<<endl;
762  ypos = ypos + step;
763  if(i == 12){ypos = step/2.;cas2 = 1;}
764  } // if
765  }
766  } // for nquat
767  cout<<"version 1: totally "<<TotalNmcp+1<<" mcps on the PD plane"<<endl;
768  */
769  //version 2: place MCPs in local bar coordinate systems:
770  //create a map of MCPs for one sector (in bar c.s.):
771  Double_t sectorWidth = 0.;
772  Int_t nmcp = 0;
773  Int_t totalnumbering = 1;
774  TVector3 location;
775  Double_t phi_curr1 = 0.;
776  for(Int_t m = 0; m < bbnum; m ++){
777  phi_curr1 = (90. - phi0 - dphi*m)/180.*pi;
778  if(m > bbnum/2-1){ phi_curr1 = (90. - phi0 - dphi*m - 2.*pipehAngle)/180.*pi; }
779  //dx_bbox = radius * cos(phi_curr1);
780  //dy_bbox = radius * sin(phi_curr1);
781  TGeoRotation rot_sector;
782  rot_sector.RotateZ( -phi0 - m*dphi - (TMath::Floor(2.*m/bbnum))*(2.*pipehAngle));
783  cout<<"phi_curr1 = "<<phi_curr1/3.1415*180.<<endl;
784  // placement of MCPs in one sector
785  cout<<"rad corner out = "<<radiusCornerOut<<", step = "<<step<<", nrowMax = "<<Int_t((sob_Rout-radiusMiddleSmall)/step)<<endl;
786  for(Int_t nrow=0; nrow < Int_t((sob_Rout-radiusMiddleSmall)/step); nrow++){
787  sectorWidth = 2.* (radiusMiddleSmall + step*nrow) * tan(dphi_rad/2.);
788  //nmcp = Int_t(sectorWidth/(MCPsize + 0.5*MCPgap));
789  nmcp = Int_t(sectorWidth/step);
790  cout<<"tg(dphi_rad/2.)= "<<tan(dphi_rad/2.)<<", nrow = "<<nrow<<", 2.* (radiusMiddleSmall + step*nrow) = "<<2.* (radiusMiddleSmall + step*nrow)<<endl;
791  cout<<"sector width = "<<sectorWidth<<", nmcp = "<<nmcp<<", nmcp%2 = "<<nmcp%2<<endl;
792  xpos = (radiusMiddleSmall + 0.5*MCPactiveArea + step*(nrow));
793  for(Int_t ny=0; ny<nmcp; ny++){
794  ypos = ((-Int_t(nmcp/2.) - 0.5*(nmcp%2))*step + (0.5+ny)*step);
795  location.SetXYZ(xpos,ypos,0.);
796  location.RotateZ(phi_curr1);
797  pdbase->AddNode(mcpgrease, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(rot_sector)));
798  pdbase->AddNode(window, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(rot_sector)));
799  pdbase->AddNode(phcathode, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(rot_sector)));
800  pdbase->AddNode(pixelholder, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick-PDwindowThick-PDsensitiveThick/2., new TGeoRotation(rot_sector)));
801  totalnumbering = totalnumbering + 1;
802  }
803  //vLocalMother->AddNode(mcp, m, new TGeoCombiTrans(dx_bbox, dy_bbox, sob_shift-PDthick/2., new TGeoRotation(rot_sector)));
804  }
805  } //for m
806  for(Int_t nrow=0; nrow < Int_t((sob_Rout-radiusMiddleSmall)/step); nrow++){
807  for(Int_t nadd = 0; nadd<2; nadd++){
808  xpos = (radiusMiddleSmall*cosFactor1 + 0.5*MCPactiveArea + step*(nrow));
809  location.SetXYZ(xpos,0.,0.);
810  location.RotateZ(pi/2.+pi*nadd);
811  pdbase->AddNode(mcpgrease, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer/2., new TGeoRotation(0)));
812  pdbase->AddNode(window, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PDwindowThick/2., new TGeoRotation(0)));
813  pdbase->AddNode(phcathode, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick/2.-PDwindowThick, new TGeoRotation(0)));
814  pdbase->AddNode(pixelholder, totalnumbering, new TGeoCombiTrans(location.X(), location.Y(), PDbaseLayer-PDgreaseLayer-PhCathodeThick-PDwindowThick-PDsensitiveThick/2., new TGeoRotation(0)));
815  totalnumbering = totalnumbering + 1;
816  }
817  }
818 
819  gGeoManager->CloseGeometry();
820 
821  top->CheckOverlaps(0.0001, "");
822  gGeoManager->CheckOverlaps(0.00001,""); // [cm]
823  //gGeoManager->CheckGeometryFull();
824  gGeoManager->SetVisLevel(4);
825  top->Write();
826  fi->Close();
827  top->Draw("ogl");
828  //pdbase->Draw("ogl");
829 
830  TObjArray *listOfOverlaps = gGeoManager->GetListOfOverlaps();
831  cout<<listOfOverlaps->GetEntries()<<endl;
832  listOfOverlaps->Print();
833  return 0;
834 }
835 
836 int findSectorIn(double y, double dphi_rad, double radius, double EVdrop, double hthick){
837  TVector3 EVcornerIn;
838  EVcornerIn.SetXYZ((radius-EVdrop-hthick)/cos(dphi_rad/2.),0.,0.);
839  Double_t Y0[4];
840  for(Int_t i=0; i<4; i++){
841  EVcornerIn.RotateZ(dphi_rad);
842  Y0[i] = EVcornerIn.Y();
843  }
844  if(y <= Y0[0]){return 1;}
845  if(y > Y0[0] && y <= Y0[1]){return 2;}
846  if(y > Y0[1] && y <= Y0[2]){return 3;}
847  if(y > Y0[2] && y <= Y0[3]){return 4;}
848  if(y > Y0[3]){return 5;}
849 
850 }
851 
852 int findSectorOut(double y, double dphi_rad, double radiusCornerOut){
853  TVector3 EVcorner;
854  EVcorner.SetXYZ(radiusCornerOut,0.,0.);
855  Double_t Y0[4];
856  for(Int_t i=0; i<4; i++){
857  EVcorner.RotateZ(dphi_rad);
858  Y0[i] = EVcorner.Y();
859  cout<<"Y ["<<i<<"] = "<<Y0[i]<<endl;
860  }
861  if(y <= Y0[0]){return 1;}
862  if(y > Y0[0] && y <= Y0[1]){return 2;}
863  if(y > Y0[1] && y <= Y0[2]){return 3;}
864  if(y > Y0[2] && y <= Y0[3]){return 4;}
865  if(y > Y0[3]){return 5;}
866 
867 }
868 
869 bool getXrange(double dphi, double radiusCornerOut, double radiusCornerIn, double MCPactiveArea, int cas1, int cas2, double ylow, double step, double &x1, double &x2){
870 
871  Double_t dphi_rad = dphi/180.*pi;
872 
873  TVector3 EVcorner;
874  EVcorner.SetXYZ(radiusCornerOut,0.,0.);
875  TVector3 EVcorner1 = EVcorner;
876  TVector3 EVcorner2 = EVcorner;
877  TVector3 EVcornerIn;
878  EVcornerIn.SetXYZ(radiusCornerIn,0.,0.);
879  TVector3 EVcornerIn1 = EVcornerIn;
880  TVector3 EVcornerIn2 = EVcornerIn;
881  //cout<<"-I- sob_Rout = "<<sob_Rout<<", dphi = "<<dphi<<", |EVcorner| = "<<EVcorner.Mag()<<", |EVcornerIn| = "<<EVcornerIn.Mag()<<endl;
882 
883  Double_t xstart, xfinish;
884  if(cas1 == 1){ // 0 - 21.6 degrees
885  EVcornerIn2.RotateZ(dphi_rad);
886  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
887  }
888  if(cas1 == 2){ // 21.6 - 43.2 degrees
889  EVcornerIn1.RotateZ(1.*dphi_rad);
890  EVcornerIn2.RotateZ(2.*dphi_rad);
891  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
892  }
893  if(cas1 == 3){ // 43.2 - 63.8 degrees
894  EVcornerIn1.RotateZ(2.*dphi_rad);
895  EVcornerIn2.RotateZ(3.*dphi_rad);
896  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
897  }
898  if(cas1 == 4){ // 63.8 - 85.6 degrees
899  EVcornerIn1.RotateZ(3.*dphi_rad);
900  EVcornerIn2.RotateZ(4.*dphi_rad);
901  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
902  }
903  if(cas1 == 5){ // 85.6 - 90.
904  EVcornerIn1.RotateZ(4.*dphi_rad);
905  EVcornerIn2.RotateZ(pi/2.);
906  xstart = getX(EVcornerIn1, EVcornerIn2, ylow) + MCPactiveArea/2.;
907  }
908 
909  if(cas2 == 1){ // 0 - 21.6 degrees
910  EVcorner2.RotateZ(dphi_rad);
911  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
912  }
913  if(cas2 == 2){ // 21.6 - 43.2 degrees
914  EVcorner1.RotateZ(1.*dphi_rad);
915  EVcorner2.RotateZ(2.*dphi_rad);
916  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
917  }
918  if(cas2 == 3){ // 43.2 - 63.8 degrees
919  EVcorner1.RotateZ(2.*dphi_rad);
920  EVcorner2.RotateZ(3.*dphi_rad);
921  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
922  }
923  if(cas2 == 4){ // 63.8 - 85.6 degrees
924  EVcorner1.RotateZ(3.*dphi_rad);
925  EVcorner2.RotateZ(4.*dphi_rad);
926  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
927  }
928  if(cas2 == 5){ // 85.6 - 90.
929  EVcorner1.RotateZ(4.*dphi_rad);
930  EVcorner2.RotateZ(pi/2.);
931  xfinish = getX(EVcorner1, EVcorner2, ylow+step) - MCPactiveArea/2.;
932  }
933  if(xstart <= 0.){xstart = step/2.;}
934  if(xstart < step/2.){xstart = step/2.;}
935  if(xfinish <= 0.){xfinish = step/2.;}
936  x1 = xstart;
937  x2 = xfinish;
938 
939  cout<<"-I- xstart = "<<xstart<<", xfinish = "<<xfinish<<endl;
940  if(x1 != x2)return true;
941  else return false;
942 }
943 
944 double getX(TVector3 corner1, TVector3 corner2, double yy){
945  //cout<<"-III- vector1"<<endl;
946  //corner1.Print();
947  //cout<<"-III- vector2 "<<endl;
948  //corner2.Print();
949  Double_t b = ( corner1.X()*corner2.Y() - corner1.Y()*corner2.X() ) / ( corner1.X() - corner2.X() );
950  Double_t k = ( corner1.Y() - corner2.Y() ) / ( corner1.X() - corner2.X() );
951  //cout<<"-II- b = "<<b<<", k = "<<k<<endl;
952  if(( yy - b )/k >= 0.) return ( yy - b )/k;
953  if(( yy - b )/k < 0.) return 0.;
954 }
CAMath::ASin
static T ASin(const T &x)
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndGeoDrc::EVbackAngle
Double_t EVbackAngle()
Definition: PndGeoDrc.h:148
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
r
double r
Definition: RiemannTest.C:14
PndGeoDrc::BBoxNum
Double_t BBoxNum()
Definition: PndGeoDrc.h:136
geoLoad
FairGeoLoader * geoLoad
Definition: createRootGeoFile.C:16
i
Int_t i
Definition: run_full.C:25
Media
FairGeoMedia * Media
Definition: createRootFscGeometryFile.C:68
m
__m128 m
Definition: P4_F32vec4.h:28
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
gGeoMan
TGeoManager * gGeoMan
Definition: createRootFscGeometryFile.C:85
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndGeoDrc::EVlen
Double_t EVlen()
Definition: PndGeoDrc.h:127
par
Double_t par[3]
Definition: anaLmdCluster.C:146
PndGeoDrc::BBoxGap
Double_t BBoxGap()
Definition: PndGeoDrc.h:130
PndGeoDrc::boxGap
Double_t boxGap()
Definition: PndGeoDrc.h:116
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndGeoDrc::barBoxZDown
Double_t barBoxZDown()
Definition: PndGeoDrc.h:104
PndGeoDrc::EVoffset
Double_t EVoffset()
Definition: PndGeoDrc.h:145
top
TGeoVolume * top
Definition: createRootFscGeometryFile.C:89
createRootGeometry_DIRC
int createRootGeometry_DIRC(Int_t fFocusingSystem=0, Bool_t fprizm=kFALSE)
Definition: createRootGeometry_DIRC.C:31
rot1
TGeoRotation * rot1
Definition: createDesySetupAli.C:115
PndGeoDrc::PipehAngle
Double_t PipehAngle()
Definition: PndGeoDrc.h:139
PndGeoDrc::PrismAngle
Double_t PrismAngle()
Definition: PndGeoDrc.h:160
b2
Double_t b2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:22
a
Int_t a
Definition: anaLmdDigi.C:126
PndGeoDrc::barBoxZUp
Double_t barBoxZUp()
Definition: PndGeoDrc.h:108
geobuild
FairGeoBuilder * geobuild
Definition: createRootFscGeometryFile.C:69
shape
TGeoShape * shape
Definition: createDummySensor.C:18
fi
TFile * fi
Definition: createRootGeoFile.C:238
PndGeoDrc::PrismDrop
Double_t PrismDrop()
Definition: PndGeoDrc.h:154
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
eps
double eps(TVector3 v1, TVector3 v2)
Definition: softtrigger_kin5.C:163
zpos
const Double_t zpos
Definition: createRootFscGeometryFile.C:11
pi
const Double_t pi
Definition: createRootGeometry_DIRC.C:24
vLocalMother
vLocalMother
Definition: createRootGeoFile.C:121
findSectorOut
int findSectorOut(double y, double dphi_rad, double radiusCornerOut)
Definition: createRootGeometry_DIRC.C:852
z
Double_t z
Definition: createRootFscGeometryFile.C:146
t3
Int_t t3
Definition: hist-t7.C:106
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndGeoDrc::boxThick
Double_t boxThick()
Definition: PndGeoDrc.h:120
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
PndGeoDrc::barHalfThick
Double_t barHalfThick()
Definition: PndGeoDrc.h:96
t2
Int_t t2
Definition: hist-t7.C:106
PndGeoDrc::PrismOffset
Double_t PrismOffset()
Definition: PndGeoDrc.h:151
PndGeoDrc::PrismhLength
Double_t PrismhLength()
Definition: PndGeoDrc.h:157
t
TTree * t
Definition: bump_analys.C:13
geoFace
FairGeoInterface * geoFace
Definition: createRootGeoFile.C:17
getXrange
void getXrange(double dphi, double radiusCornerOut, double radiusCornerIn, double MCPactiveArea, int cas1, int cas2, double ylow, double step, double &x1, double &x2)
Definition: createRootGeometry_DIRC.C:869
findSectorIn
int findSectorIn(double y, double dphi, double radius, double EVdrop, double hthick)
Definition: createRootGeometry_DIRC.C:836
y
Double_t y
Definition: createRootFscGeometryFile.C:159
alpha
double alpha
Definition: f_Init.h:9
fGeoFile
TString fGeoFile
Definition: createRootGeoFile.C:26
getX
double getX(TVector3 corner1, TVector3 corner2, double yy)
Definition: createRootGeometry_DIRC.C:944
PndGeoDrc::barNum
Double_t barNum()
Definition: PndGeoDrc.h:124
r2
double r2
Definition: RiemannTrackFinderTest.C:19
PndGeoDrc::radius
Double_t radius()
Definition: PndGeoDrc.h:92
Generated on Sun Aug 18 2019 03:01:09 for FairRoot/PandaRoot by   doxygen 1.8.8