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