FairRoot/PandaRoot
PndLmdDim.cxx
Go to the documentation of this file.
1 /*
2  * Pndcpp
3  *
4  * this class gives you methods to retrieve parameters of the geometry
5  *
6  * Created on: Oct 5, 2012
7  * Author: promme
8  * as of 2015 maintained by Roman Klasen, klasen@kph.uni-mainz.de
9  */
10 
11 #include <PndLmdDim.h>
12 #include<TGeoMatrix.h>
13 #include <stdlib.h>
14 #include <TGeoPhysicalNode.h>
15 // //work with DB
16 // #include<PndLmdContFact.h>
17 // #include<TList.h>
18 // #include<PndLmdAlignPar.h>
19 PndLmdDim* PndLmdDim::pinstance = 0;
20 
21 #include <cmath>
22 
23 // version 3: navigation path names changed
24 int PndLmdDim::geometry_version = 3;
25 
26 #include <TROOT.h>
27 PndLmdDim::PndLmdDim() {
28  sensIDoffset = 0;
29  double test_mult_fact = 1.; //100.; // should be 1 when not debugging code
30  pi = 3.141592654;
31  // number of detector planes
32  n_planes = 4;
33  // number of modules per plane half
34  nmodules = 5;
35  // position of planes where the first plane defines the origin
36  plane_pos_z = new double[4];
37  plane_pos_z[0] = 0.0;
38  plane_pos_z[1] = 20.0;
39  plane_pos_z[2] = 30.0;
40  plane_pos_z[3] = 40.0;
41  half_offset_x = 0.005 * test_mult_fact; // 50 mum
42  half_offset_y = 0.005 * test_mult_fact;
43  half_offset_z = 0.005 * test_mult_fact;
44  half_tilt_phi = 0.; // negligible
45  half_tilt_theta = 0.;
46  half_tilt_psi = 0.;
47  plane_half_offset_x = 0.005 * test_mult_fact;
48  plane_half_offset_y = 0.005 * test_mult_fact;
49  plane_half_offset_z = 0.005 * test_mult_fact;
50  plane_half_tilt_phi = 0.0; // negligible
51  plane_half_tilt_theta = 0.0;
52  plane_half_tilt_psi = 0.0;
53  cvd_offset_x = 0.005 * test_mult_fact;
54  cvd_offset_y = 0.005 * test_mult_fact;
55  cvd_offset_z = 0.0;
56  cvd_tilt_phi = 0.001 * test_mult_fact;
57  cvd_tilt_theta = 0.;
58  cvd_tilt_psi = 0.;
59  side_offset_x = 0.005 * test_mult_fact;
60  side_offset_y = 0.005 * test_mult_fact;
61  side_offset_z = 0.; // do not use! -> clashing volumes
62  side_tilt_phi = 0.; // do not use! -> clashing volumes
63  side_tilt_theta = 0.;
64  side_tilt_psi = 0.; // do not use! -> clashing volumes
65  sensor_offset_x = 0*100e-4 * test_mult_fact; // 100 mum
66  sensor_offset_y = 0*100e-4 * test_mult_fact; // 100 mum
67  sensor_offset_z = 0.; // do not use! -> clashing volumes
68  sensor_tilt_phi = 0*2500e-6 * test_mult_fact; // 150 murad
69  sensor_tilt_theta = 0.; // do not use! -> clashing volumes
70  sensor_tilt_psi = 0.; // do not use! -> clashing volumes
71 
72 
73 
74  // cvd_diamond is cut out of 79.5 mm discs of 200 micron thickness
75  // inner min. radius due to beam pipe + a safety margin
76  inner_rad = 3.7;
77  // not used yet but should be the outer acceptance
78  outer_rad = 12.;
79  // number of CVD diamond discs per plane
80  n_cvd_discs = 10;
81  // radius of a CVD diamond disc
82  cvd_disc_rad = 7.95 / 2.;
83  // the half of the diamond thickness
84  cvd_disc_thick_half = 0.01;
85  // even and odd discs in a plane will be shifted in z in order to prevent
86  // mechanical damage during assembly
87  cvd_disc_even_odd_offset = 0.25;
88  // angle from the division of a circle into n_cvd_discs
89  delta_phi = 2. * pi / ((double) (n_cvd_discs));
90  // a polygon of n_cvd_discs sides fitting a radius of inner_rad
91  // has a side length pol_side_lg of
92  pol_side_lg_half = inner_rad * sin(delta_phi / 2.);
93  // the minimum distance to the center of the polygone is given by
94  pol_side_dist_min = inner_rad * cos(delta_phi / 2.);
95  // the cvd disc has to be placed such that the disc crosses
96  // the inner ring at an angle of 0 and delta_phi
97  // this defines the distance to the center according to pythagoras
98  cvd_disc_dist = pol_side_dist_min
99  + sqrt(cvd_disc_rad * cvd_disc_rad - pol_side_lg_half * pol_side_lg_half);
100 
101  kapton_disc_thick_half = 0.004454 / 2.; // Using aluminum in equivalent thickness to cover flex cable and gluing at once
102 
103  maps_n_col = 3;
104  maps_n_row = 2;
105  // enabled [row][col]
106  enabled = new bool*[maps_n_row];
107  enabled[0] = new bool[maps_n_col];
108  enabled[1] = new bool[maps_n_col];
109  enabled[0][0] = true;
110  enabled[0][1] = true;
111  enabled[0][2] = true;
112  enabled[1][0] = false;
113  enabled[1][1] = true;
114  enabled[1][2] = true;
115  n_sensors = 5;
116  // NOTE: MOST of the following VARIABLES are HALF of it
117  // due to geometry construction in GEANT
118  maps_thickness = 0.005 / 2.;
119  maps_passive_top = 0.01 / 2.; //*12000.;
120  maps_passive_bottom = 0.05 / 2.;
121  maps_passive_left = 0.01 / 2.; //*6000.;
122  maps_passive_right = 0.01 / 2.; //*3000.;
123  maps_active_width = 2.0 / 2. - maps_passive_left - maps_passive_right;
124  maps_active_height = 2.0 / 2. - maps_passive_top - maps_passive_bottom;
125 
126  maps_active_pixel_size = 0.008; // 80µm pixel size
127 
128  maps_width = maps_passive_left + maps_active_width + maps_passive_right;
129  maps_height = maps_passive_top + maps_active_height + maps_passive_bottom;
130 
131  maps_active_offset_x = (maps_passive_left - maps_passive_right) / 2.;
132  maps_active_offset_y = (maps_passive_top - maps_passive_bottom) / 2.;
133 
134  die_gap = 0.05; // (cm)
135 
136  maps_die_width = maps_width * maps_n_col;
137  maps_die_height = maps_height * maps_n_row + die_gap * maps_n_row - 1;
138 
139  die_offset_x = 0.;
140  die_offset_y = 0.;
141  die_offset_z = 0.;
142  die_tilt_phi = 0.;
143  die_tilt_theta = 0.;
144  die_tilt_psi = 0.;
145  //*********************************** lumi box parameters ***********************************
146  // see CAD files for details
147  // https://edms.cern.ch/nav/P:FAIR-000000719:V0/P:FAIR-000000726:V0/TAB3
148  // width
149  box_size_x = 36. / 2.;
150  // height
151  box_size_y = 60.6 / 2.;
152  // length
153  box_size_z = 95.6 / 2.;
154  // thickness of the V2A steel plates
155  box_thickness = 0.5;
156  // position of the inner rib
157  pos_rib = 36.4;
158  // beam pipe radius at entrance
159  rad_entrance = 9. + 1.;
160  // beam pipe radius at exit
161  rad_exit = 5.;
162  // beam pipe separating non interacting paricles
163  rad_pipe = 3.5;
164  // beam pipe thickness;
165  //pipe_thickness = 0.1;
166  // cone height of the transition region
167  //length_transision = 36.4;
168  // length of the inner pipe
169  //length_pipe = 50.;
170  //*********************************** global parameters *************************************
171  // where bending starts with
172  end_seg_upstream = 361;
173  // the bending radius
174  r_bend = 5700.;
175  // and the angle of the circle path
176  phi_bend = 40.0e-3;
177  ;
178  // the point where both tangents of the straight beam pipe tubes meet is
179  pos_rot_z = end_seg_upstream + tan(phi_bend / 2.) * r_bend;
180  // z position of the lmd box
181  pos_z = 1050.; //(cm)
182  // position of the first detector plane
183  pos_plane_0 = 74.1;
184  end_seg_bend = end_seg_upstream + sin(phi_bend) * r_bend;
185  // x position of the lmd
186  pos_x = (pos_z - end_seg_upstream - tan(phi_bend / 2.) * r_bend)
187  * tan(phi_bend);
188  // y position of the lmd
189  pos_y = 0.;
190  rot_x = 0.;
191  rot_y = phi_bend;
192  rot_z = 0.;
193 
194  // later transformations will require a relative
195  // navigation structure based on strings
196  // the path to individual volumes is stored here
197 
198  // luminosity detector top box
199  nav_paths.push_back("lmd_vol_vac");
200  // luminosity detector reference system
201  nav_paths.push_back("lmd_vol_ref_sys");
202  // luminosity detector halfs with the
203  // number 0 for top and 1 for bottom
204  nav_paths.push_back("lmd_vol_half");
205  // luminosity detector plane 0 to 3
206  nav_paths.push_back("lmd_vol_plane");
207  // luminosity detector module 0 to 5
208  // clockwise around direction upstream (z)
209  nav_paths.push_back("lmd_vol_module");
210  // luminosity detector front side (upstream = 0)
211  // and backside (downstream = 1)
212  nav_paths.push_back("lmd_vol_side");
213  // luminosity detector die
214  // ( 0 for the 1x3 sensors and 1 for the 1x2 sensors )
215  nav_paths.push_back("lmd_vol_die");
216  // luminosity detector sensor
217  //nav_paths.push_back("lmd_vol_sensor_"); misalignment of individual sensors is not supportet yet
218  // luminosity detector active sensor
219  // 0 is the most inner sensor
220  nav_paths.push_back("LumActivePixelRect"); // the _ comes automatically with the copy number
221 
222  fgGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
223  if (!fgGeoMan)
224  cout << "Error: could not find the geometry manager!" << endl;
225 }
226 
227 PndLmdDim::PndLmdDim(const PndLmdDim & instance) {
228  if (!instance.box_thickness)
229  cout << " fix for pedantic compiler will never appear on cout " << endl;
230 }
231 
232 PndLmdDim::~PndLmdDim() {
233  cout << endl;
234  cout << " Cleaning up PndLmdDim " << endl;
235  cout
236  << " If you see that message several times please check the performance of your code! "
237  << endl;
238  Cleanup();
239  delete pinstance;
240 }
241 
242 PndLmdDim & PndLmdDim::Get_instance() {
243  if (!pinstance) {
244  pinstance = new PndLmdDim();
245  }
246  return (*pinstance);
247 }
248 
249 PndLmdDim* PndLmdDim::Instance() {
250  if (!pinstance) {
251  pinstance = new PndLmdDim();
252  }
253  return (pinstance);
254 }
255 
256 /*
257  void PndLmdDim::transform_sensor_local_to_lmd_local(const int sensor_id, double & x, double & y, double & z, bool misaligned)
258  {
259  }
260 
261 
262 
263  void PndLmdDim::transform_local_lmd_to_global(const int sensor_id, double & x, double & y, double & z, bool misaligned)
264  {
265  }
266 
267 
268 
269  void PndLmdDim::transform_local_sensor()
270  {
271  }*/
272 
273 #include<FairGeoLoader.h>
274 #include<FairGeoInterface.h>
275 #include<FairGeoMedia.h>
276 #include<FairGeoBuilder.h>
277 #include<TGeoTorus.h>
278 #include<TGeoSphere.h>
279 #include<TGeoTube.h>
280 #include<TGeoCompositeShape.h>
281 #include<TGeoCone.h>
282 #include<TGeoBBox.h>
283 #include <TGeoPcon.h>
284 
285 #include<TPad.h>
286 bool PndLmdDim::Retrieve_version_number() {
287  bool result = false;
288  //FairRootManager* ioman = FairRootManager::Instance();
289  TGeoManager* gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
290  if (!gGeoMan) {
291  cout << " Info: no FAIRGeom found, using gGeoManager " << endl;
292  gGeoMan = gGeoManager;
293  }
294  if (gGeoMan) {
295  //gGeoMan->Draw("ogl");
296  //gPad->Update();
297  //gPad->Print("test.pdf");
298  TGeoVolume* vol = gGeoMan->FindVolumeFast(nav_paths[0].c_str());
299  //cout << gGeoMan->GetCurrentNode()->GetName() << endl;
300  if (vol) {
301  string sversion = vol->GetTitle();
302  if (sversion.compare(0, 8, "version ") != 0) {
303  cout
304  << " Warning from PndLmdDim::Retrieve_version_number: no version number encoded in the node title. Setting it to 0. "
305  << endl;
306  geometry_version = 0;
307  } else {
308  // take the number behind "version " string
309  geometry_version = atoi(&(vol->GetTitle()[8]));
310  }
311  cout
312  << " Info from PndLmdDim::Retrieve_version_number: The geometry version was set to "
313  << geometry_version << endl;
314  //cout << vol->GetName() << endl;
315  result = true;
316  } else {
317  cout << " *** Error in PndLmdDim::Retrieve_version_number:" << endl;
318  cout << " Could not find the top volume " << nav_paths[0].c_str()
319  << " to retrieve the version number of the luminosity detector! Is the geometry already loaded? "
320  << endl;
321  }
322  } else {
323  cout << " *** Error in PndLmdDim::Retrieve_version_number:" << endl;
324  cout
325  << " Could not find a GeoManager to retrieve the version number of the luminosity detector! "
326  << endl;
327  }
328  return result;
329 }
330 
331 void PndLmdDim::Generate_rootgeom(TGeoVolume& mothervol, bool misaligned) {
332  Cleanup();
333  FairGeoLoader* geoLoad = FairGeoLoader::Instance();
334  if (!geoLoad) {
335  geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
336  cout << " creating FairGeoLoader instance " << endl;
337  }
338  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
339  string dir = getenv("VMCWORKDIR");
340  geoFace->setMediaFile((dir + "/geometry/media_pnd.geo").c_str()); //("${VMCWORKDIR}/geometry/media_pnd.geo");
341  geoFace->readMedia();
342  //geoFace->print();
343 
344  FairGeoMedia *Media = geoFace->getMedia();
345  FairGeoBuilder *geobuild = geoLoad->getGeoBuilder();
346 
347  // retrieve available media
348  FairGeoMedium* FairMediumAir = Media->getMedium("air");
349  FairGeoMedium* FairMediumSteel = Media->getMedium("steel");
350  FairGeoMedium* FairMediumAl = Media->getMedium("Aluminum");
351  FairGeoMedium *FairMediumSilicon = Media->getMedium("silicon");
352  FairGeoMedium *FairMediumDiamond = Media->getMedium("HYPdiamond");
353  FairGeoMedium *FairMediumVacuum = Media->getMedium("vacuum7");
354  FairGeoMedium *FairMediumKapton = Media->getMedium("kapton");
355  FairGeoMedium *FairMediumCopper = Media->getMedium("copper");
356 
357  if (!FairMediumAir || !FairMediumSteel || !FairMediumAl || !FairMediumKapton
358  || !FairMediumSilicon || !FairMediumVacuum || !FairMediumCopper) {
359  std::cout << " Error: not all media found " << std::endl;
360  return;
361  }
362 
363  int nmed;
364  nmed = geobuild->createMedium(FairMediumAir);
365  nmed = geobuild->createMedium(FairMediumSteel);
366  nmed = geobuild->createMedium(FairMediumAl);
367  nmed = geobuild->createMedium(FairMediumKapton);
368  nmed = geobuild->createMedium(FairMediumSilicon);
369  nmed = geobuild->createMedium(FairMediumDiamond);
370  nmed = geobuild->createMedium(FairMediumVacuum);
371  if (!nmed)
372  cout << " fix for pedantic compiler line " << endl;
373 
374  // no translation nor rotation
375  TGeoRotation* rot_no = new TGeoRotation("rot_no", 0., 0., 0.); // no rotation
376  TGeoTranslation* trans_no = new TGeoTranslation("trans_no", 0., 0., 0.); // no translation
377  TGeoCombiTrans* rottrans_no = new TGeoCombiTrans("rottrans_no", 0., 0., 0.,
378  rot_no);
379  rottrans_no->RegisterYourself();
380 
381  // ************ create the luminosity monitor box ***********
382  // create the bounding box
383  double tube_upstream_length = 33.3 / 2.;
384  double tube_upstream_rad_in = 10.;
385  double tube_upstream_rad_out = tube_upstream_rad_in + 0.2;
386  double tube_dostream_length = 12. / 2.; // length with flange
387 
388  // create a vacuum around the luminosity detector
389  double lmd_total_length = box_size_z + tube_upstream_length
390  + tube_dostream_length;
391 
392  double origin[3] = { 0., 0., lmd_total_length };
393  TGeoBBox* lmd_box_vac = new TGeoBBox("lmd_box_vac", box_size_x, box_size_y,
394  lmd_total_length, origin);
395  TGeoVolume *lmd_vol_vac = new TGeoVolume(nav_paths[0].c_str(), lmd_box_vac,
396  fgGeoMan->GetMedium("vacuum7"));
397  //lmd_vol_vac->SetTransparency(20);
398  lmd_vol_vac->SetLineColor(3);
399  double x, y, z, rottheta, rotphi, rotpsi;
400  Get_pos_lmd_global(x, y, z, rottheta, rotphi, rotpsi);
401  TGeoRotation* lmd_rot = new TGeoRotation("lmd_rot");
402  lmd_rot->RotateX(rottheta / pi * 180.);
403  lmd_rot->RotateY(rotphi / pi * 180.);
404  lmd_rot->RotateZ(rotpsi / pi * 180.);
405  TGeoCombiTrans* lmd_transrot = new TGeoCombiTrans(x, y, z, lmd_rot);
406  lmd_transrot->SetName("lmd_transrot");
407  lmd_transrot->RegisterYourself();
408  // pumping station upstream
409  TGeoTube* pipe_upstream = new TGeoTube("pipe_upstream", tube_upstream_rad_in,
410  tube_upstream_rad_out, tube_upstream_length);
411  TGeoCombiTrans* comb_trans_pipe_upstream = new TGeoCombiTrans(
412  "comb_trans_pipe_upstream", 0., 0., tube_upstream_length, rot_no);
413  //comb_trans_pipe_upstream->RegisterYourself();
414  TGeoVolume* lmd_vol_pipe_up = new TGeoVolume("lmd_vol_pipe_up", pipe_upstream,
415  fgGeoMan->GetMedium("steel"));
416  lmd_vol_pipe_up->SetLineColor(11);
417  lmd_vol_vac->AddNode(lmd_vol_pipe_up, 0, comb_trans_pipe_upstream);
418  // the lmd box
419  TGeoBBox* lmd_box_outer = new TGeoBBox("lmd_box_outer", box_size_x,
420  box_size_y, box_size_z);
421  TGeoBBox* lmd_box_inner = new TGeoBBox("lmd_box_inner",
422  box_size_x - box_thickness, box_size_y - box_thickness,
423  box_size_z - box_thickness);
424  TGeoBBox* lmd_box_rib = new TGeoBBox("lmd_box_rib",
425  box_size_x - box_thickness, box_size_y - box_thickness,
426  box_thickness / 2.);
427  TGeoTube* box_hole_upstream = new TGeoTube("box_hole_upstream", 0.0,
428  rad_entrance, box_thickness);
429  // move the cut pipe upstream
430  TGeoCombiTrans* comb_trans_cut_pipe_upstream = new TGeoCombiTrans(
431  "comb_trans_cut_pipe_upstream", 0., 0., -box_size_z + box_thickness / 2.,
432  rot_no);
433  comb_trans_cut_pipe_upstream->RegisterYourself();
434  TGeoTube* box_hole_downstream = new TGeoTube("box_hole_downstream", 0.0,
435  rad_exit, box_thickness);
436  // move the cut pipe downstream
437  TGeoCombiTrans* comb_trans_cut_pipe_downstream = new TGeoCombiTrans(
438  "comb_trans_cut_pipe_downstream", 0., 0.,
439  +box_size_z - box_thickness / 2., rot_no);
440  comb_trans_cut_pipe_downstream->RegisterYourself();
441  // move rib upstream
442  TGeoCombiTrans* comb_trans_rib = new TGeoCombiTrans("comb_trans_rib", 0., 0.,
443  -box_size_z + pos_rib, rot_no);
444  comb_trans_rib->RegisterYourself();
445  // inner clash protection rods
446  double clash_rod_x = 3.;
447  double clash_rod_y = 0.5;
448  double clash_rod_z = 29.75;
449  double origin_left[3] =
450  { -box_size_x + clash_rod_x + box_thickness, 0., 18.15 };
451  new TGeoBBox("lmd_box_clash_rod_left",
452  clash_rod_x, clash_rod_y, clash_rod_z, origin_left);// TGeoBBox* lmd_box_clash_rod_left = //[R.K.03/2017] unused variable
453  double origin_right[3] = { +box_size_x - clash_rod_x - box_thickness, 0.,
454  18.15 };
455  new TGeoBBox("lmd_box_clash_rod_right",
456  clash_rod_x, clash_rod_y, clash_rod_z, origin_right);//TGeoBBox* lmd_box_clash_rod_right = //[R.K.03/2017] unused variable
457  // compose all the parts into one luminosity vacuum box
458  TGeoCompositeShape *shape_lmd_box =
459  new TGeoCompositeShape("shape_lmd_box",
460  "(lmd_box_outer-lmd_box_inner + ((lmd_box_rib-box_hole_upstream):comb_trans_rib))"
461  "-box_hole_upstream:comb_trans_cut_pipe_upstream"
462  "-box_hole_downstream:comb_trans_cut_pipe_downstream"
463  "+lmd_box_clash_rod_left+lmd_box_clash_rod_right");
464  TGeoVolume *lmd_vol_box = new TGeoVolume("lmd_vol_box", shape_lmd_box,
465  fgGeoMan->GetMedium("steel"));
466  lmd_vol_box->SetLineColor(11);
467  //lmd_vol_box->SetVisibility(false);//TEST
468  //lmd_vol_box->SetTransparency(20);
469  TGeoCombiTrans* comb_trans_lmd_box = new TGeoCombiTrans("comb_trans_lmd_box",
470  0., 0., 2 * tube_upstream_length + box_size_z, rot_no);
471  //comb_trans_pipe_upstream->RegisterYourself();
472  lmd_vol_vac->AddNode(lmd_vol_box, 0, comb_trans_lmd_box);
473 
474  // pipe downstream of the box
475  double lmd_pipe_box_do[18];
476  lmd_pipe_box_do[0] = 0.;
477  lmd_pipe_box_do[1] = 360.;
478  lmd_pipe_box_do[2] = 4.;
479 
480  lmd_pipe_box_do[3] = box_size_z * 2. + tube_upstream_length * 2.;
481  lmd_pipe_box_do[4] = rad_exit - 0.2;
482  lmd_pipe_box_do[5] = rad_exit;
483 
484  lmd_pipe_box_do[6] = lmd_pipe_box_do[3] + tube_dostream_length * 2. - 4.;
485  lmd_pipe_box_do[7] = lmd_pipe_box_do[4];
486  lmd_pipe_box_do[8] = lmd_pipe_box_do[5];
487 
488  lmd_pipe_box_do[9] = lmd_pipe_box_do[6];
489  lmd_pipe_box_do[10] = lmd_pipe_box_do[4];
490  lmd_pipe_box_do[11] = 7.;
491 
492  lmd_pipe_box_do[12] = lmd_pipe_box_do[3] + tube_dostream_length * 2.;
493  lmd_pipe_box_do[13] = lmd_pipe_box_do[10];
494  lmd_pipe_box_do[14] = lmd_pipe_box_do[11];
495 
496  //lmd_pipe_box_do[15] = lmd_pipe_box_do[12];
497  //lmd_pipe_box_do[16] = lmd_pipe_box_do[4];
498  //lmd_pipe_box_do[17] = lmd_pipe_box_do[14];
499 
500  TGeoPcon* shape_pipe_box_do = new TGeoPcon(lmd_pipe_box_do);
501  TGeoVolume* vlum_pipe_box_do = new TGeoVolume("vlum_pipe_box_do",
502  shape_pipe_box_do, fgGeoMan->GetMedium("steel"));
503  vlum_pipe_box_do->SetLineColor(kGray); //39);
504  TGeoCombiTrans* lmd_trans_pipe_box_do = new TGeoCombiTrans(
505  "lmd_trans_pipe_box_do", 0., 0., 0., rot_no);
506  lmd_trans_pipe_box_do->RegisterYourself();
507  lmd_vol_vac->AddNode(vlum_pipe_box_do, 0, lmd_trans_pipe_box_do);
508  // end of pipe downstream of the box
509 
510  // TGeoTube* lmd_flange_upstr = new TGeoTube(
511  // "lmd_flange_upstr", 9.2, 25.3/2., 1.2);
512  // TGeoCombiTrans* lmd_trans_fl_up = new TGeoCombiTrans("lmd_trans_fl_up", 0., 0., 1.2+delta, r1);
513  // lmd_trans_fl_up->RegisterYourself();
514  // upstream flange holding the kapton cone
515  // TGeoTube* lmd_cone_flange_upstr = new TGeoTube(
516  // "lmd_cone_flange_upstr", 9.2, 25.3/2., 1.5);
517  // TGeoCombiTrans* lmd_trans_co_fl_up = new TGeoCombiTrans("lmd_trans_co_fl_up", 0., 0., -1.5+50.-delta, r1);
518  // lmd_trans_co_fl_up->RegisterYourself();
519 
520  // ********************* enhanced beam pipe construct in the box ********************
521 
522  double debug_multiplier = 1.;
523  // define where the inner region of the box starts
524  double box_inner_up_z = 2 * tube_upstream_length + box_thickness;
525 
526  // flange holding the kapton cone upstream
527  double lmd_cone_clamp_length = 13.;
528  double lmd_cone_clamp_params[15];
529  lmd_cone_clamp_params[0] = 0.;
530  lmd_cone_clamp_params[1] = 360.;
531  lmd_cone_clamp_params[2] = 4.;
532 
533  // values taken from drawings starting upstream at the box
534  // pipe extension as well as female part of the clamp are included as well
535  lmd_cone_clamp_params[3] = 0.;
536  lmd_cone_clamp_params[4] = 20. / 2.;
537  lmd_cone_clamp_params[5] = 20.4 / 2.;
538 
539  lmd_cone_clamp_params[6] = 10.;
540  lmd_cone_clamp_params[7] = 20. / 2.;
541  lmd_cone_clamp_params[8] = 20.4 / 2.;
542 
543  lmd_cone_clamp_params[9] = 10.;
544  lmd_cone_clamp_params[10] = 20. / 2.;
545  lmd_cone_clamp_params[11] = 26. / 2.;
546 
547  lmd_cone_clamp_params[12] = lmd_cone_clamp_length;
548  lmd_cone_clamp_params[13] = 20. / 2.;
549  lmd_cone_clamp_params[14] = 26. / 2.;
550 
551  TGeoPcon* lmd_cone_clamp = new TGeoPcon(lmd_cone_clamp_params);
552  TGeoVolume* vlum_cone_clamp = new TGeoVolume("vlum_cone_clamp",
553  lmd_cone_clamp, fgGeoMan->GetMedium("steel"));
554  vlum_cone_clamp->SetLineColor(kGray); //39);
555  TGeoCombiTrans* lmd_trans_lmd_cone_clamp = new TGeoCombiTrans(
556  "lmd_trans_lmd_cone_clamp", 0., 0., box_inner_up_z, rot_no);
557  lmd_trans_lmd_cone_clamp->RegisterYourself();
558  lmd_vol_vac->AddNode(vlum_cone_clamp, 0, lmd_trans_lmd_cone_clamp);
559 
560  /*
561  * polymide |\
562  * aluminum |\\
563  * \\\|\
564  * \\\ \
565  * \\\ \
566  * \\| \ ________
567  * |\|\__________| V2A tube
568  * \____| brass
569  *
570  * y [cm]
571  * ^
572  * |
573  * |-------> z [cm]
574  *
575  *
576  * nomenclature:
577  * transistion region consists of an aluminum (AL) / polymide (pol) laminate cone (c)
578  * with an inner (i) and outer (o) diameter upstream (u) and downstream (d)
579  * cpoluo = cone polymide upstream outer
580  * The tube (t) is made of V2A (A2) and starts with a cone to glue the transistion cone to it
581  * To conduct the EM fields of the non interacting beam a brass (br) ring is inserted
582  * The opening angle of the cone is 15 degree as determined by HESR
583  * the tube must have a diameter of 70 mm due to beam acceptance restrictions
584  */
585 
586  // 10µm aluminum
587  double cALthick = 0.001 * debug_multiplier;
588  // 20µm polimide
589  double cpolthick = 0.002 * debug_multiplier;
590  // 500µm brass ring
591  double tbrthick = 0.05;
592  // 500µm inner beam pipe
593  double tA2thick = 0.05;
594  double cone_angle = 15. / 180. * 3.1416;
595 
596  double cALdiy = 7. / 2.;
597  double cALdiz = 40.;
598  double cALdoy = cALdiy + cALthick / cos(cone_angle);
599  double cALdoz = cALdiz;
600  double cpoldiy = cALdoy;
601  double cpoldiz = cALdoz;
602  double cpoldoy = cpoldiy + cpolthick / cos(cone_angle);
603  double cpoldoz = cpoldiz;
604  double cbrdoy = cALdiy;
605  double cbrdoz = cALdiz;
606  double cbrdiy = cbrdoy - tbrthick; // /cos(cone_angle) does not apply here for simplicity;
607  double cbrdiz = cbrdoz;
608 
609  double cALuiz = lmd_cone_clamp_length;
610  double cALuiy = cALdiy + (cALdiz - cALuiz) * tan(cone_angle);
611  double cALuoy = cALuiy + cALthick / cos(cone_angle);
612  double cALuoz = cALuiz;
613  double cpoluiy = cALuoy;
614  double cpoluiz = cALuoz;
615  double cpoluoy = cpoluiy + cpolthick / cos(cone_angle);
616  double cpoluoz = cpoluiz;
617 
618  double cbruoz = cbrdiz - 0.5;
619  double cbruoy = cbrdoy + 0.5 * tan(cone_angle);
620  double cbruiz = cbruoz;
621  double cbruiy = cbruoy - tbrthick; // /cos(cone_angle) does not apply here for simplicity;
622 
623  double cA2diy = cALdiy;
624  double cA2diz = cpoldoz + (cpoldoy - cALdiy) / tan(cone_angle);
625  double cA2uiz = cA2diz - 0.1 / tan(cone_angle);
626  double cA2uiy = cA2diy + 0.1;
627  double cA2uoy = cA2uiy + tA2thick;
628  //double cA2uoz = cA2uiz; //[R.K. 01/2017] unused variable?
629  double cA2doy = cA2diy + tA2thick;
630  //double cA2doz = cA2diz; //[R.K.02/2017] Unused variable?
631 
632  double tbruiz = cbrdiz;
633  double tbruiy = cbrdiy;
634  //double tbruoz = cbrdoz; //[R.K. 01/2017] unused variable?
635  double tbruoy = cbrdoy;
636 
637  double tbrdiz = tbruiz + 2.;
638  double tbrdiy = tbruiy;
639  //double tbrdoz = tbrdiz; //[R.K. 01/2017] unused variable?
640  double tbrdoy = tbruoy;
641 
642  //double tA2uoz = cA2doz; //[R.K. 01/2017] unused variable?
643  double tA2uoy = cA2doy;
644  double tA2uiz = cA2diz;
645  double tA2uiy = cA2diy;
646 
647  double tA2diz = tA2uiz + 44.;
648  double tA2diy = tA2uiy;
649  //double tA2doz = tA2uiz; //[R.K. 01/2017] unused variable?
650  double tA2doy = tA2uoy;
651 
652  TGeoCone* lmd_capton_cone = new TGeoCone("lmd_capton_cone",
653  (cpoldoz - cpoluoz) / 2., cpoluiy, cpoluoy, cpoldiy, cpoldoy);
654  TGeoCombiTrans* lmd_trans_cap_co = new TGeoCombiTrans("lmd_trans_cap_co", 0.,
655  0., box_inner_up_z + (cpoldoz - cpoluoz) / 2. + cpoluoz, rot_no);
656  lmd_trans_cap_co->RegisterYourself();
657  TGeoVolume *vlum_CaptonCone = new TGeoVolume("vlum_CaptonCone",
658  lmd_capton_cone, fgGeoMan->GetMedium("kapton"));
659  vlum_CaptonCone->SetLineColor(kRed); //39);
660  lmd_vol_vac->AddNode(vlum_CaptonCone, 0, lmd_trans_cap_co); //TEST with/without cone!!!
661 
662  // 10 mu thick kapton foil aluminum coating
663  TGeoCone* lmd_al_cone = new TGeoCone("lmd_al_cone", (cALdiz - cALuiz) / 2.,
664  cALuiy, cALuoy, cALdiy, cALdoy);
665  TGeoVolume *vlum_AlCone = new TGeoVolume("vlum_AlCone", lmd_al_cone,
666  fgGeoMan->GetMedium("Aluminum"));
667  vlum_AlCone->SetLineColor(kGray); //39);
668  lmd_vol_vac->AddNode(vlum_AlCone, 0, lmd_trans_cap_co); //TEST with/without cone!!!
669 
670  double lmd_pipe_params[36];
671  lmd_pipe_params[0] = 0.;
672  lmd_pipe_params[1] = 360.;
673  lmd_pipe_params[2] = 11.;
674 
675  lmd_pipe_params[3] = cA2uiz;
676  lmd_pipe_params[4] = cA2uiy;
677  lmd_pipe_params[5] = cA2uoy;
678 
679  lmd_pipe_params[6] = cA2diz;
680  lmd_pipe_params[7] = cA2diy;
681  lmd_pipe_params[8] = cA2doy;
682 
683  lmd_pipe_params[9] = tA2diz;
684  lmd_pipe_params[10] = tA2diy;
685  lmd_pipe_params[11] = tA2doy;
686 
687  lmd_pipe_params[12] = 0 + lmd_pipe_params[9];
688  lmd_pipe_params[13] = 3.5;
689  lmd_pipe_params[14] = 4.05;
690 
691  lmd_pipe_params[15] = 1.5 + lmd_pipe_params[9];
692  lmd_pipe_params[16] = 3.5;
693  lmd_pipe_params[17] = 4.05;
694 
695  lmd_pipe_params[18] = 3.5 + lmd_pipe_params[9];
696  lmd_pipe_params[19] = 4.2;
697  lmd_pipe_params[20] = 6.;
698 
699  lmd_pipe_params[21] = 10. + lmd_pipe_params[9];
700  lmd_pipe_params[22] = 4.2;
701  lmd_pipe_params[23] = 6.;
702 
703  lmd_pipe_params[24] = 10. + lmd_pipe_params[9];
704  lmd_pipe_params[25] = 4.2;
705  lmd_pipe_params[26] = 4.55;
706 
707  lmd_pipe_params[27] = 20. + lmd_pipe_params[9];
708  lmd_pipe_params[28] = 4.2;
709  lmd_pipe_params[29] = 4.55;
710 
711  lmd_pipe_params[30] = 20. + lmd_pipe_params[9];
712  lmd_pipe_params[31] = 4.2;
713  lmd_pipe_params[32] = 4.55; //7.;
714 
715  lmd_pipe_params[33] = 22 + lmd_pipe_params[9];
716  lmd_pipe_params[34] = 4.2;
717  lmd_pipe_params[35] = 4.55; //7.;
718 
719  TGeoPcon* lmd_V2_pipe = new TGeoPcon(lmd_pipe_params);
720  TGeoVolume* vlum_V2_pipe = new TGeoVolume("vlum_V2_pipe", lmd_V2_pipe,
721  fgGeoMan->GetMedium("steel"));
722  vlum_V2_pipe->SetLineColor(kGray); //39);
723  TGeoCombiTrans* lmd_trans_lmd_V2_pipe = new TGeoCombiTrans(
724  "lmd_trans_lmd_V2_pipe", 0., 0., box_inner_up_z, rot_no);
725  lmd_trans_lmd_V2_pipe->RegisterYourself();
726  lmd_vol_vac->AddNode(vlum_V2_pipe, 0, lmd_trans_lmd_V2_pipe);
727 
728  double lmd_conductor_params[12];
729  lmd_conductor_params[0] = 0.;
730  lmd_conductor_params[1] = 360.;
731  lmd_conductor_params[2] = 3.;
732 
733  lmd_conductor_params[3] = cbruiz;
734  lmd_conductor_params[4] = cbruiy;
735  lmd_conductor_params[5] = cbruoy;
736 
737  lmd_conductor_params[6] = cbrdiz;
738  lmd_conductor_params[7] = cbrdiy;
739  lmd_conductor_params[8] = cbrdoy;
740 
741  lmd_conductor_params[9] = tbrdiz;
742  lmd_conductor_params[10] = tbrdiy;
743  lmd_conductor_params[11] = tbrdoy;
744 
745  TGeoPcon* lmd_cond_ring = new TGeoPcon(lmd_conductor_params);
746  TGeoVolume* vlum_cond_ring = new TGeoVolume("vlum_cond_ring", lmd_cond_ring,
747  fgGeoMan->GetMedium("steel")); // should be brass
748  vlum_cond_ring->SetLineColor(kYellow); //39);
749  TGeoCombiTrans* lmd_trans_lmd_cond_ring = new TGeoCombiTrans(
750  "lmd_trans_lmd_cond_ring", 0., 0., box_inner_up_z, rot_no);
751  lmd_trans_lmd_cond_ring->RegisterYourself();
752  lmd_vol_vac->AddNode(vlum_cond_ring, 0, lmd_trans_lmd_cond_ring);
753  // beam pipe to shield the sensors
754  // TGeoTube* lmd_beam_pipe = new TGeoTube("lmd_beam_pipe", 3.5, 3.6, 50./2.);
755  // TGeoCombiTrans* lmd_trans_p = new TGeoCombiTrans("lmd_trans_p", 0., 0., 50.+23.386+50./2., r1);
756  // lmd_trans_p->RegisterYourself();
757  // beam pipe cone downstream
758  // TGeoCone* lmd_cone_downstr = new TGeoCone("lmd_cone_downstr", 20./2., 3.5, 3.7, 9./2., 9.2/2.);
759  // TGeoCombiTrans* lmd_trans_co_do = new TGeoCombiTrans("lmd_trans_co_do", 0., 0., 50.+23.386+50.+20./2., r1);
760  // lmd_trans_co_do->RegisterYourself();
761  // beam pipe downstream
762  // TGeoTube* lmd_beam_pipe_downstream = new TGeoTube("lmd_beam_pipe_downstream", 9./2., 9.2/2., 56./2.);
763  // TGeoCombiTrans* lmd_trans_p_down = new TGeoCombiTrans("lmd_trans_p_down", 0., 0., 50.+23.386+50.+20.+56./2., r1);
764  // lmd_trans_p_down->RegisterYourself();*/
765 
766  // ****************************** luminosity detector reference system ************************
767  // definition of the luminosity detector local system
768  TGeoCombiTrans* rottrans_lmd_in_box = new TGeoCombiTrans(
769  "rottrans_lmd_in_box", 0., 0., pos_plane_0, rot_no);
770  TGeoVolumeAssembly* lmd_vol_ref_sys = new TGeoVolumeAssembly(
771  nav_paths[1].c_str());
772  // definition of the retractable luminosity detector halves
773 
774  // ****************************** cvd cooling support structure ********************
775 
776  // According to Heinrich's drawings from 07.03.14
777  double lmd_cool_sup_inner_rad = 10.6;
778  double lmd_cool_sup_outer_rad = 16.;
779  double lmd_cool_sup_outer_cut = 14.;
780  double lmd_cool_sup_thick = 1.; // half of it
781 
782  // construct first a tube
783  new TGeoTube("shape_cool_sup_tube",
784  lmd_cool_sup_inner_rad, lmd_cool_sup_outer_rad, lmd_cool_sup_thick);//TGeoTube* shape_cool_sup_tube = //[R.K.03/2017] unused variable
785  // to cut off a half + a little bit
786  new TGeoBBox("shape_cool_sup_cut",
787  lmd_cool_sup_outer_rad, lmd_cool_sup_outer_rad, lmd_cool_sup_thick + 0.1);//TGeoBBox* shape_cool_sup_cut = //[R.K.03/2017] unused variable
788  // set the position for the cut off
789  TGeoTranslation* trans_shape_cool_sup_cut_low = new TGeoTranslation(
790  "trans_shape_cool_sup_cut_low", 0., -lmd_cool_sup_outer_rad + 0.5, 0.); // 0.5 should be 0.3 but for the sake of simplicity not to clash to simple rod description
791  TGeoCombiTrans* combtrans_shape_cool_sup_cut_low = new TGeoCombiTrans(
792  *trans_shape_cool_sup_cut_low, *rot_no);
793  combtrans_shape_cool_sup_cut_low->SetName("combtrans_shape_cool_sup_cut_low");
794  combtrans_shape_cool_sup_cut_low->RegisterYourself();
795 
796  TGeoTranslation* trans_shape_cool_sup_cut_high = new TGeoTranslation(
797  "trans_shape_cool_sup_cut_high", 0., +lmd_cool_sup_outer_rad - 0.5, 0.); // 0.5 should be 0.3 but for the sake of simplicity not to clash to simple rod description
798  TGeoCombiTrans* combtrans_shape_cool_sup_cut_high = new TGeoCombiTrans(
799  *trans_shape_cool_sup_cut_high, *rot_no);
800  combtrans_shape_cool_sup_cut_high->SetName(
801  "combtrans_shape_cool_sup_cut_high");
802  combtrans_shape_cool_sup_cut_high->RegisterYourself();
803 
804  // We need some cut outs for the modules and the outer structure
805  // we give them a little bit more space for misalignment studies without clashing volumes
806  new TGeoTube("shape_module_cutout", 0.,
807  cvd_disc_rad + 0.05,
808  cvd_disc_thick_half + 2 * kapton_disc_thick_half + 0.01);//TGeoTube* shape_module_cutout = //[R.K.03/2017] unused variable
809  for (size_t imodule = 0; imodule < nmodules * 2 /*upper and lower half*/;
810  imodule++) {
811  double angle = delta_phi / 2. + imodule * delta_phi;
812  double add_z = cvd_disc_even_odd_offset;
813  // the offset of the modules in the upper and lower halfs
814  // are opposite
815  if (((imodule) % 2) == 0)
816  add_z = -add_z;
817  double _x = cos(angle) * cvd_disc_dist;
818  double _y = sin(angle) * cvd_disc_dist;
819  TGeoTranslation* trans_shape_module_cutout = new TGeoTranslation(_x, _y,
820  add_z);
821  TGeoCombiTrans* combtrans_shape_module_cutout = new TGeoCombiTrans(
822  *trans_shape_module_cutout, *rot_no);
823  stringstream _cutout_name;
824  _cutout_name << "rottrans_cutout_" << imodule;
825  combtrans_shape_module_cutout->SetName(_cutout_name.str().c_str());
826  combtrans_shape_module_cutout->RegisterYourself();
827 
828  _x = cos(angle) * (lmd_cool_sup_outer_cut + lmd_cool_sup_outer_rad);
829  _y = sin(angle) * (lmd_cool_sup_outer_cut + lmd_cool_sup_outer_rad);
830  TGeoTranslation* trans_cool_sup_cut_outer = new TGeoTranslation(_x, _y, 0.);
831  stringstream _cutshape_rot_name;
832  _cutshape_rot_name << "cutshaperot_" << imodule;
833  TGeoRotation* rot_cool_sup_cut_outer = new TGeoRotation(
834  _cutshape_rot_name.str().c_str(), angle / pi * 180., 0., 0.);
835  TGeoCombiTrans* combtrans_cool_sup_cut_outer = new TGeoCombiTrans(
836  *trans_cool_sup_cut_outer, *rot_cool_sup_cut_outer);
837  stringstream _cutshape_name;
838  _cutshape_name << "cutshape_" << imodule;
839  combtrans_cool_sup_cut_outer->SetName(_cutshape_name.str().c_str());
840  combtrans_cool_sup_cut_outer->RegisterYourself();
841  }
842 
843  // construct the support from basic shape and it's cut outs
844  TGeoCompositeShape *shape_cool_sup_up = new TGeoCompositeShape(
845  "shape_cool_sup_up",
846  "shape_cool_sup_tube-shape_cool_sup_cut:combtrans_shape_cool_sup_cut_low"
847  "-shape_module_cutout:rottrans_cutout_0"
848  "-shape_module_cutout:rottrans_cutout_1"
849  "-shape_module_cutout:rottrans_cutout_2"
850  "-shape_module_cutout:rottrans_cutout_3"
851  "-shape_module_cutout:rottrans_cutout_4"
852  "-shape_cool_sup_cut:cutshape_0"
853  "-shape_cool_sup_cut:cutshape_1"
854  "-shape_cool_sup_cut:cutshape_2"
855  "-shape_cool_sup_cut:cutshape_3"
856  "-shape_cool_sup_cut:cutshape_4");
857 
858  TGeoVolume* lmd_vol_cool_sup_up = new TGeoVolume("lmd_vol_cool_sup_up",
859  shape_cool_sup_up, fgGeoMan->GetMedium("Aluminum"));
860  lmd_vol_cool_sup_up->SetLineColor(17);
861 
862  // the lower cooling support has to be rotated, we create simply a new volume for it
863  TGeoRotation* rot_shape_cool_sup_down = new TGeoRotation(
864  "rot_shape_cool_sup_down", 180., 180., 0.);
865  TGeoCombiTrans* combtrans_shape_cool_sup_down = new TGeoCombiTrans(*trans_no,
866  *rot_shape_cool_sup_down);
867  combtrans_shape_cool_sup_down->SetName("combtrans_shape_cool_sup_down");
868  combtrans_shape_cool_sup_down->RegisterYourself();
869  // construct the support from basic shape and it's cut outs
870  TGeoCompositeShape *shape_cool_sup_down = new TGeoCompositeShape(
871  "shape_cool_sup_down",
872  "shape_cool_sup_tube-shape_cool_sup_cut:combtrans_shape_cool_sup_cut_high"
873  "-shape_module_cutout:rottrans_cutout_5"
874  "-shape_module_cutout:rottrans_cutout_6"
875  "-shape_module_cutout:rottrans_cutout_7"
876  "-shape_module_cutout:rottrans_cutout_8"
877  "-shape_module_cutout:rottrans_cutout_9"
878  "-shape_cool_sup_cut:cutshape_5"
879  "-shape_cool_sup_cut:cutshape_6"
880  "-shape_cool_sup_cut:cutshape_7"
881  "-shape_cool_sup_cut:cutshape_8"
882  "-shape_cool_sup_cut:cutshape_9");
883 
884  TGeoVolume* lmd_vol_cool_sup_down = new TGeoVolume("lmd_vol_cool_sup_down",
885  shape_cool_sup_down, fgGeoMan->GetMedium("Aluminum"));
886  lmd_vol_cool_sup_down->SetLineColor(17);
887 
888  // ****************************** cvd cooling support discs ************************
889  // the cvd disc shape
890  double gap_between_disc_and_support_structure(0.025); // 250 mu gap
891  TGeoTubeSeg* shape_cvd_disc = new TGeoTubeSeg("shape_cvd_disc", inner_rad,
892  lmd_cool_sup_inner_rad - gap_between_disc_and_support_structure,
893  cvd_disc_thick_half, -delta_phi / 2. / pi * 180.,
894  +delta_phi / 2. / pi * 180.);
895 
896  TGeoRotation* cvd_rotation = new TGeoRotation("cvd_rotation", 0, 0, 0);
897  TGeoTranslation* cvd_translation = new TGeoTranslation("cvd_translation",
898  -cvd_disc_dist, 0, 0);
899  TGeoCombiTrans* cvd_combtrans = new TGeoCombiTrans(*cvd_translation,
900  *cvd_rotation);
901  cvd_combtrans->SetName("cvd_combtrans");
902  cvd_combtrans->RegisterYourself();
903 
904  //this next line is pretty stupid but it made the work for the better geometry minimal
905  //otherwise I would have to do some deeper digging and reworking...
906  TGeoCompositeShape *shape_cvd_support = new TGeoCompositeShape(
907  "shape_cvd_support",
908  "(shape_cvd_disc:cvd_combtrans+shape_cvd_disc:cvd_combtrans)");
909 
910  TGeoVolume* lmd_vol_cvd_disc = new TGeoVolume("lmd_vol_cvd_disc",
911  shape_cvd_support, fgGeoMan->GetMedium("HYPdiamond"));
912  lmd_vol_cvd_disc->SetLineColor(9);
913  // ****************************** kapton flexible circuits to the sensors ************************
914  // the cvd disc shape
915  TGeoTubeSeg* shape_kapton_disc = new TGeoTubeSeg("shape_kapton_disc",
916  inner_rad,
917  lmd_cool_sup_inner_rad - gap_between_disc_and_support_structure,
918  kapton_disc_thick_half, -delta_phi / 2. / pi * 180.,
919  +delta_phi / 2. / pi * 180.);
920 
921  TGeoRotation* kapton_rotation = new TGeoRotation("kapton_rotation", 0, 0, 0);
922  TGeoTranslation* kapton_translation = new TGeoTranslation(
923  "kapton_translation", -cvd_disc_dist, 0, 0);
924  TGeoCombiTrans* kapton_combtrans = new TGeoCombiTrans(*kapton_translation,
925  *kapton_rotation);
926  kapton_combtrans->SetName("kapton_combtrans");
927  kapton_combtrans->RegisterYourself();
928 
929  //this next line is pretty stupid but it made the work for the better geometry minimal
930  //otherwise I would have to do some deeper digging and reworking...
931  TGeoCompositeShape *shape_kapton_support =
932  new TGeoCompositeShape("shape_kapton_support",
933  "(shape_kapton_disc:kapton_combtrans+shape_kapton_disc:kapton_combtrans)");
934 
935  TGeoVolume* lmd_vol_kapton_disc = new TGeoVolume("lmd_vol_kapton_disc",
936  shape_kapton_support, fgGeoMan->GetMedium("Aluminum")); //kapton")); // changed to equivalent for glue/flex cable etc.
937  lmd_vol_kapton_disc->SetLineColor(kRed);
938  //lmd_vol_kapton_disc->SetTransparency(50);
939  lmd_vol_kapton_disc->SetVisibility(false);
940  // *********************************** HV-MAPS *************************************
941 
942  // create basic shapes and their positions
943  TGeoBBox *shape_maps_active_centered = new TGeoBBox(
944  "shape_maps_active_centered", maps_active_width, maps_active_height,
945  maps_thickness);
946  TGeoCombiTrans* combtrans_maps_active = new TGeoCombiTrans(
947  "combtrans_maps_active",
948  -maps_width + maps_passive_left * 2. + maps_active_width,
949  -maps_height + maps_passive_bottom * 2. + maps_active_height, 0., rot_no);
950  combtrans_maps_active->RegisterYourself();
951  TGeoBBox *shape_maps_passive_left = new TGeoBBox("shape_maps_passive_left",
952  maps_passive_left, maps_height, maps_thickness);
953  TGeoCombiTrans* combtrans_maps_passive_left = new TGeoCombiTrans(
954  "combtrans_maps_passive_left", -maps_width + maps_passive_left, 0., 0.,
955  rot_no);
956  combtrans_maps_passive_left->RegisterYourself();
957  TGeoBBox *shape_maps_passive_right = new TGeoBBox("shape_maps_passive_right",
958  maps_passive_right, maps_height, maps_thickness);
959  TGeoCombiTrans* combtrans_maps_passive_right = new TGeoCombiTrans(
960  "combtrans_maps_passive_right", maps_width - maps_passive_right, 0., 0.,
961  rot_no);
962  combtrans_maps_passive_right->RegisterYourself();
963  TGeoBBox *shape_maps_passive_top = new TGeoBBox("shape_maps_passive_top",
964  maps_width, maps_passive_top, maps_thickness);
965  TGeoCombiTrans* combtrans_maps_passive_top = new TGeoCombiTrans(
966  "combtrans_maps_passive_top", 0., maps_height - maps_passive_top, 0.,
967  rot_no);
968  combtrans_maps_passive_top->RegisterYourself();
969  TGeoBBox *shape_maps_passive_bottom = new TGeoBBox(
970  "shape_maps_passive_bottom", maps_width, maps_passive_bottom,
971  maps_thickness);
972  TGeoCombiTrans* combtrans_maps_passive_bottom = new TGeoCombiTrans(
973  "combtrans_maps_passive_bottom", 0., -maps_height + maps_passive_bottom,
974  0., rot_no);
975  combtrans_maps_passive_bottom->RegisterYourself();
976 
977  if (!lmd_box_outer || !lmd_box_inner || !lmd_box_rib || !box_hole_upstream
978  || !box_hole_downstream || !shape_cvd_disc || !shape_kapton_disc
979  || !shape_maps_active_centered || !shape_maps_passive_left
980  || !shape_maps_passive_right || !shape_maps_passive_top
981  || !shape_maps_passive_bottom)
982  cout << " pedantic compiler together with root geometries sucks " << endl;
983 
984  TGeoCompositeShape *shape_maps_passive =
985  new TGeoCompositeShape("shape_maps_passive",
986  "(shape_maps_passive_top:combtrans_maps_passive_top+shape_maps_passive_right:combtrans_maps_passive_right+shape_maps_passive_bottom:combtrans_maps_passive_bottom+shape_maps_passive_left:combtrans_maps_passive_left)");
987 
988  TGeoCompositeShape *shape_maps_active = new TGeoCompositeShape(
989  "shape_maps_active",
990  "(shape_maps_active_centered:combtrans_maps_active-shape_maps_passive)");
991 
992  TGeoVolume* _vol_passive = new TGeoVolume("LumPassiveRect",
993  shape_maps_passive, fgGeoMan->GetMedium("silicon"));
994  _vol_passive->SetLineColor(30);
995 
996  TGeoVolume* _vol_active = new TGeoVolume(nav_paths[7].c_str(),
997  shape_maps_active, fgGeoMan->GetMedium("silicon"));
998  _vol_active->SetLineColor(36);
999  // **************************************************************
1000 
1001  // ****************************** loops in the luminosity detector ************************
1002  stringstream name;
1003  stringstream uniqueid; // seems pandaroot has problems when volumes are not uniquely named
1004  double _x(0), _y(0), _z(0), _rotphi(0), _rottheta(0), _rotpsi(0);
1005  double _offset_x(0), _offset_y(0), _offset_z(0), _offset_phi(0),
1006  _offset_theta(0), _offset_psi(0);
1007  unsigned int sensor_id(0);
1008  unsigned int module_id(0);
1009  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
1010  // in order do be able to displace the detector halves those are introduced as
1011  // separate volume assemblies
1012  name.str("");
1013  name << nav_paths[2]; // << ihalf;
1014  uniqueid.str("");
1015  //uniqueid << "_" << ihalf;
1016  TGeoVolumeAssembly* lmd_vol_half_ = new TGeoVolumeAssembly(
1017  name.str().c_str());
1018  //mothervol.AddNode(lmd_vol_half_, 0, rottrans_no);
1019  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
1020  name.str("");
1021  name << nav_paths[3]; // << iplane;
1022  uniqueid.str("");
1023  //uniqueid << "_" << ihalf << iplane;
1024  TGeoVolumeAssembly* lmd_vol_plane_ = new TGeoVolumeAssembly(
1025  (name.str() + uniqueid.str()).c_str());
1026  // move to the position of the corresponding plane
1027  TGeoMatrix* rottrans_plane = new TGeoCombiTrans(0., 0.,
1028  plane_pos_z[iplane], rot_no);
1029  if (misaligned) {
1030  Get_offset(ihalf, iplane, -1, -1, -1, -1, _offset_x, _offset_y,
1031  _offset_z, _offset_phi, _offset_theta, _offset_psi, true);
1032  TGeoRotation* rot_plane_offset = new TGeoRotation("rot_plane_offset",
1033  _offset_phi / pi * 180., _offset_theta / pi * 180.,
1034  _offset_psi / pi * 180.);
1035  TGeoCombiTrans* rottrans_plane_offset = new TGeoCombiTrans(_offset_x,
1036  _offset_y, _offset_z, rot_plane_offset);
1037  // rottrans_plane = new TGeoHMatrix(*rottrans_plane * *rottrans_plane_offset);
1038  rottrans_plane = new TGeoHMatrix(
1039  *rottrans_plane_offset * *rottrans_plane);
1040  }
1041  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
1042  name.str("");
1043  name << nav_paths[4]; // << imodule;
1044  uniqueid.str("");
1045  //uniqueid << "_" << ihalf << iplane << imodule;
1046  TGeoVolumeAssembly* lmd_vol_module_ = new TGeoVolumeAssembly(
1047  (name.str() + uniqueid.str()).c_str());
1048  double angle = delta_phi / 2. + ihalf * pi + imodule * delta_phi;
1049  double add_z = cvd_disc_even_odd_offset;
1050  // the offset of the modules in the upper and lower halfs
1051  // are opposite
1052  if (((imodule + ihalf) % 2) == 0)
1053  add_z = -add_z;
1054  _x = cos(angle) * cvd_disc_dist;
1055  _y = sin(angle) * cvd_disc_dist;
1056  _z = add_z;
1057  _rotphi = 0.;
1058  _rottheta = 0.;
1059  _rotpsi = angle / pi * 180.;
1060  TGeoRotation* rot_module = new TGeoRotation("rot_module", _rotphi,
1061  _rottheta, _rotpsi);
1062  TGeoMatrix* rottrans_module = new TGeoCombiTrans(_x, _y, _z,
1063  rot_module);
1064  if (misaligned) {
1065  Get_offset(ihalf, iplane, imodule, -1, -1, -1, _offset_x, _offset_y,
1066  _offset_z, _offset_phi, _offset_theta, _offset_psi, true);
1067  TGeoRotation* rot_module_offset = new TGeoRotation(
1068  "rot_module_offset", _offset_phi / pi * 180.,
1069  _offset_theta / pi * 180., _offset_psi / pi * 180.);
1070  TGeoCombiTrans* rottrans_module_offset = new TGeoCombiTrans(_offset_x,
1071  _offset_y, _offset_z, rot_module_offset);
1072  rottrans_module = new TGeoHMatrix(
1073  *rottrans_module_offset * *rottrans_module);
1074  }
1075  /*
1076  if (ihalf == 0 && iplane == 2 && imodule == 3){
1077  _x += 2.;
1078  _y += 1.;
1079  _z += -5.;
1080  _rotphi += 20.;
1081  _rottheta += 70.;
1082  }*/
1083  // add the cvd disc into that assembly
1084  //TGeoVolume* lmd_vol_cvd_disc = new TGeoVolume("lmd_vol_cvd_disc",
1085  // shape_cvd_support, fgGeoMan->GetMedium("HYPdiamond"));
1086  lmd_vol_module_->AddNode(lmd_vol_cvd_disc, module_id, rottrans_no);
1087  module_id++;
1088  for (unsigned int iside = 0; iside < 2; iside++) { // loop over the two sides of the modules
1089  name.str("");
1090  name << nav_paths[5]; // << iside;
1091  uniqueid.str("");
1092  //uniqueid << "_" << ihalf << iplane << imodule << iside;
1093  TGeoVolumeAssembly* lmd_vol_side_ = new TGeoVolumeAssembly(
1094  (name.str() + uniqueid.str()).c_str());
1095  // rotation around the y axis for the upstream side
1096  // side 0 is downstream! what may be not so obvious
1097  _x = 0.;
1098  _y = 0.;
1099  _z = 0.;
1100  if (iside == 0) {
1101  _rotphi = 0.;
1102  _rottheta = 0.;
1103  _rotpsi = 0.;
1104  } else {
1105  _rotphi = 0.;
1106  _rottheta = 180.;
1107  _rotpsi = 0.;
1108  }
1109  TGeoRotation* rot_side = new TGeoRotation("rot_side", _rotphi,
1110  _rottheta, _rotpsi);
1111  TGeoMatrix* rottrans_side = new TGeoCombiTrans(_x, _y, _z, rot_side);
1112  if (misaligned) {
1113  Get_offset(ihalf, iplane, imodule, iside, -1, -1, _offset_x,
1114  _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi,
1115  true);
1116  TGeoRotation* rot_side_offset = new TGeoRotation("rot_side_offset",
1117  _offset_phi / pi * 180., _offset_theta / pi * 180.,
1118  _offset_psi / pi * 180.);
1119  TGeoCombiTrans* rottrans_side_offset = new TGeoCombiTrans(_offset_x,
1120  _offset_y, _offset_z, rot_side_offset);
1121  rottrans_side = new TGeoHMatrix(
1122  *rottrans_side_offset * *rottrans_side);
1123  }
1124  // glue the HV-MAPS to the cvd surface
1125  for (unsigned int idie = 0; idie < 2; idie++) { // loop over dies
1126  // it seems we will have only 2 sensors on a die and only a few of them
1127  // -> term die will stay as it is, but may be different from the reality
1128  name.str("");
1129  name << nav_paths[6]; // << idie;
1130  uniqueid.str("");
1131  //uniqueid << "_" << ihalf << iplane << imodule << iside << idie;
1132  TGeoVolumeAssembly* lmd_vol_die_ = new TGeoVolumeAssembly(
1133  (name.str() + uniqueid.str()).c_str());
1134  // rotation to the cut side of the cvd_disc
1135  // the origin is the inner edge
1136  const double _sinhalf = sin(delta_phi / 2.);
1137  //const double _coshalf = cos (delta_phi/2.);
1138  const double _edge_y = +_sinhalf * inner_rad;
1139  // angle between the edge and the half of the cvd disc
1140  const double _edgeangle = asin(-_edge_y / cvd_disc_rad);
1141  const double _edge_x = -cvd_disc_rad * cos(_edgeangle);
1142  _x = _edge_x;
1143  _y = _edge_y;
1144  _z = +cvd_disc_thick_half + maps_thickness; // move to the surface
1145  _rotphi = +delta_phi / 2. / pi * 180.;
1146  _rottheta = 0.;
1147  _rotpsi = 0.; // rotate to the cut edge
1148  TGeoRotation* rot_die = new TGeoRotation("rot_die", _rotphi,
1149  _rottheta, _rotpsi);
1150  TGeoCombiTrans* rottrans_die = new TGeoCombiTrans(_x, _y, _z,
1151  rot_die);
1152  // construct now the sensors on the two dies
1153  for (unsigned int isensor = 0; isensor < 3; isensor++) { // loop over sensors
1154  // the 0th die is oriented at the inner edge
1155  _x = +maps_width + maps_width * 2. * isensor;
1156  _y = -maps_height;
1157  _z = 0.;
1158  // only the inner two sensors will possibly on one physical die
1159  // so take an offset into account
1160  if (isensor == 2)
1161  _x = _x + die_gap;
1162  // next row
1163  if (idie == 1) {
1164  if (isensor == 0)
1165  continue;
1166  _y -= die_gap + 2. * maps_height;
1167  _rotphi = 0.;
1168  } else {
1169  _rotphi = 0.;
1170  }
1171  _rottheta = 0.;
1172  _rotpsi = 0.;
1173  //"LumActiveRect" is the keyword for digitization of hits
1174  /*
1175  name.str("");
1176  name << nav_paths[7] << isensor;
1177  uniqueid.str("");
1178  uniqueid << "_" << ihalf << iplane << imodule << iside << idie << isensor;
1179  TGeoVolume* _vol_active =
1180  new TGeoVolume(
1181  (name.str()+uniqueid.str()).c_str(),
1182  shape_maps_active,
1183  fgGeoMan->GetMedium("silicon"));
1184  _vol_active->SetLineColor(36);
1185  */
1186  /*
1187  name.str("");
1188  name << "LumPassiveRect_" << isensor;
1189  TGeoVolume* _vol_passive =
1190  new TGeoVolume(
1191  (name.str()+uniqueid.str()).c_str(),
1192  shape_maps_passive,
1193  fgGeoMan->GetMedium("silicon"));
1194  _vol_passive->SetLineColor(30);
1195  */
1196  TGeoRotation* rot_sensor = new TGeoRotation("rot_sensor", _rotphi,
1197  _rottheta, _rotpsi);
1198  TGeoMatrix* rottrans_sensor = new TGeoCombiTrans(_x, _y, _z,
1199  rot_sensor);
1200  if (misaligned) {
1201  Get_offset(ihalf, iplane, imodule, iside, idie, isensor,
1202  _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta,
1203  _offset_psi, true);
1204  TGeoRotation* rot_sensor_offset = new TGeoRotation(
1205  "rot_sensor_offset", _offset_phi / pi * 180.,
1206  _offset_theta / pi * 180., _offset_psi / pi * 180.);
1207  TGeoCombiTrans* rottrans_sensor_offset = new TGeoCombiTrans(
1208  _offset_x, _offset_y, _offset_z, rot_sensor_offset);
1209  rottrans_sensor = new TGeoHMatrix(
1210  *rottrans_sensor_offset * *rottrans_sensor);
1211  }
1212  lmd_vol_die_->AddNode(_vol_active, sensor_id, rottrans_sensor);
1213  lmd_vol_die_->AddNode(_vol_passive, sensor_id, rottrans_sensor);
1214  //cout << sensor_id << " " << _vol_active->GetName() << endl << endl;
1215  // save the transformation from the cvd_side reference frame
1216  // into the local frame of the sensors
1217  transformation_matrices[Tkey(ihalf, iplane, imodule, iside, idie,
1218  isensor)] = new TGeoHMatrix(
1219  (*rottrans_die) * (*rottrans_sensor));
1220 
1221  if (0) { // some tests for debugging
1222  unsigned int _sensor_id = Get_sensor_id(ihalf, iplane, imodule,
1223  iside, idie, isensor);
1224  if (sensor_id != _sensor_id) {
1225  cout << " wrong sensor id " << _sensor_id << " != "
1226  << sensor_id << endl;
1227  }
1228  int _ihalf, _iplane, _imodule, _iside, _idie, _isensor;
1229  Get_sensor_by_id(sensor_id, _ihalf, _iplane, _imodule, _iside,
1230  _idie, _isensor);
1231  if ((signed) ihalf != _ihalf)
1232  cout << " wrong half " << _ihalf << endl;
1233  if ((signed) iplane != _iplane)
1234  cout << " wrong plane " << _iplane << endl;
1235  if ((signed) imodule != _imodule)
1236  cout << " wrong module " << _imodule << endl;
1237  if ((signed) iside != _iside)
1238  cout << " wrong side " << _iside << endl;
1239  if ((signed) idie != _idie)
1240  cout << " wrong die " << _idie << endl;
1241  if ((signed) isensor != _isensor)
1242  cout << " wrong sensor " << _isensor << endl;
1243  }
1244  sensor_id++;
1245  } // loop over sensors
1246  lmd_vol_side_->AddNode(lmd_vol_die_, idie, rottrans_die);
1247  } // loop over dies
1248  _x = 0;
1249  _y = 0;
1250  _z = +cvd_disc_thick_half + maps_thickness * 2.
1251  + kapton_disc_thick_half; // move to the surface
1252  _rotphi = 0.;
1253  _rottheta = 0.;
1254  _rotpsi = 0.; // rotate to the cut edge
1255  TGeoRotation* rot_kapton = new TGeoRotation("rot_kapton", _rotphi,
1256  _rottheta, _rotpsi);
1257  TGeoCombiTrans* rottrans_kapton = new TGeoCombiTrans(_x, _y, _z,
1258  rot_kapton);
1259  lmd_vol_side_->AddNode(lmd_vol_kapton_disc, 0, rottrans_kapton); // GGenerate_Tkeyumber(ihalf,iplane, imodule, iside)
1260  lmd_vol_module_->AddNode(lmd_vol_side_, iside, rottrans_side);
1261  // save the transformation from the lumi reference frame
1262  // into the local cvd side reference frame
1263  transformation_matrices[Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1264  new TGeoHMatrix(
1265  (*rottrans_plane) * (*rottrans_module) * (*rottrans_side));
1266  } // loop over the two sides of the modules
1267  lmd_vol_plane_->AddNode(lmd_vol_module_, imodule, rottrans_module);
1268  // transformation_matrices[Generate_Tkey(ihalf, iplane, imodule, -1, -1, -1)] =
1269  // new TGeoHMatrix((*rottrans_plane) * (*rottrans_module));
1270  } // loop over modules
1271  if (ihalf == 0)
1272  lmd_vol_plane_->AddNode(lmd_vol_cool_sup_up, iplane, rottrans_no);
1273  else
1274  lmd_vol_plane_->AddNode(lmd_vol_cool_sup_down, iplane, rottrans_no);
1275  lmd_vol_half_->AddNode(lmd_vol_plane_, iplane, rottrans_plane);
1276  // // save the transformation from the lumi reference frame
1277  // // into the plane reference frame
1278  // transformation_matrices[Generate_Tkey(ihalf, iplane, -1, -1, -1, -1)] =
1279  // new TGeoHMatrix((*rottrans_plane));
1280  } // loop over planes
1281  lmd_vol_ref_sys->AddNode(lmd_vol_half_, ihalf, rottrans_no);
1282  // save the transformation into the lumi reference frame
1283  transformation_matrices[Tkey(-1, -1, -1, -1, -1, -1)] = new TGeoHMatrix(
1284  (*lmd_transrot) * (*rottrans_lmd_in_box));
1285  } // loop over detector halves
1286  // code geometry version in the node title
1287  stringstream nodetitle;
1288  nodetitle << "version " << geometry_version << endl;
1289  lmd_vol_vac->SetTitle(nodetitle.str().c_str());
1290  lmd_vol_vac->AddNode(lmd_vol_ref_sys, 0, rottrans_lmd_in_box);
1291  mothervol.AddNode(lmd_vol_vac, geometry_version, lmd_transrot);
1292 
1293  /*
1294  gGeoMan->CloseGeometry();
1295  //gGeoMan->Get
1296  //cout << gGeoMan->InitTrack(0,30,1050,0,0,1)->GetName() << endl;
1297  const double* current_point = gGeoMan->GetCurrentPoint();
1298  cout << current_point[0] << " " << current_point[1] << " " << current_point[2] << endl;
1299  double new_point[3] = {8.7, 1., 29.7375};
1300  gGeoMan->SetCurrentPoint(new_point);
1301  const double* new_current_point = gGeoMan->GetCurrentPoint();
1302  cout << new_current_point[0] << " " << new_current_point[1] << " " << new_current_point[2] << endl;
1303  cout << gGeoMan->GetCurrentNode()->GetName() << endl;
1304  gGeoMan->FindNode();
1305  cout << gGeoMan->GetCurrentNode()->GetName() << endl;
1306 
1307  string path = "/lmd_HV_MAPS_1/vol_lmd_vac_0/Lumi_HV-MAPS_0/vol_LumPassive_cvd_disc_plane_1_disc_1_0";
1308  gGeoMan->cd(path.c_str());///LumActivePixelRect_plane_0_disc_0_side_0_col_0_row_0_0");
1309  cout << gGeoMan->GetPath()<< endl;
1310  TGeoPhysicalNode* node =
1311  gGeoMan->MakePhysicalNode((path + "/LumActivePixelRect_plane_0_disc_0_side_0_col_0_row_0_0").c_str());///Lumi_HV-MAPS/LumPassive_cvd_disc_plane_2_disc_3");
1312  if (node){
1313  //TGeoCombiTrans* align_trans = new TGeoCombiTrans(1.,0.,0.,georot_no);
1314  //node->Align(align_trans);
1315  //gGeoMan->CloseGeometry();
1316  TGeoHMatrix* matrix = node -> GetMatrix();
1317  matrix -> Print();
1318  }
1319  //gGeoMan->Set
1320  */
1321 }
1322 
1323 //void PndLmdDim::Create_transformation_matrices_aligned(string filename_in, string filename_out){
1324 void PndLmdDim::reCreate_transformation_matrices() {
1325  // no translation nor rotation
1326  TGeoRotation* rot_no = new TGeoRotation("rot_no", 0., 0., 0.); // no rotation
1327  //TGeoCombiTrans* rottrans_no = new TGeoCombiTrans("rottrans_no", 0., 0.,
1328  // 0., rot_no);
1329 
1330  double _x(0), _y(0), _z(0), _rotphi(0), _rottheta(0), _rotpsi(0);
1331  double _offset_x(0), _offset_y(0), _offset_z(0), _offset_phi(0),
1332  _offset_theta(0), _offset_psi(0);
1333 
1334  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
1335  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
1336  // move to the position of the corresponding plane
1337  TGeoMatrix* rottrans_plane = new TGeoCombiTrans(0., 0.,
1338  plane_pos_z[iplane], rot_no);
1339 
1340  // Get_offset(ihalf, iplane, -1, -1, -1, -1,
1341  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1342  // TGeoRotation* rot_plane_offset = new TGeoRotation("rot_plane_offset",
1343  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1344  // TGeoCombiTrans* rottrans_plane_offset =
1345  // new TGeoCombiTrans(_offset_x, _offset_y, _offset_z, rot_plane_offset);
1346  // // rottrans_plane = new TGeoHMatrix(*rottrans_plane * *rottrans_plane_offset);
1347  // rottrans_plane = new TGeoHMatrix(*rottrans_plane_offset * *rottrans_plane);
1348 
1349  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
1350  double angle = delta_phi / 2. + ihalf * pi + imodule * delta_phi;
1351  double add_z = cvd_disc_even_odd_offset;
1352  // the offset of the modules in the upper and lower halfs
1353  // are opposite
1354  if (((imodule + ihalf) % 2) == 0)
1355  add_z = -add_z;
1356  _x = cos(angle) * cvd_disc_dist;
1357  _y = sin(angle) * cvd_disc_dist;
1358  _z = add_z;
1359  _rotphi = 0.;
1360  _rottheta = 0.;
1361  _rotpsi = angle / pi * 180.;
1362  TGeoRotation* rot_module = new TGeoRotation("rot_module", _rotphi,
1363  _rottheta, _rotpsi);
1364  TGeoMatrix* rottrans_module = new TGeoCombiTrans(_x, _y, _z,
1365  rot_module);
1366  // Set_offset(ihalf, iplane, imodule, -1, -1, -1,
1367  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1368  Get_offset(ihalf, iplane, imodule, -1, -1, -1, _offset_x, _offset_y,
1369  _offset_z, _offset_phi, _offset_theta, _offset_psi);
1370  // TGeoRotation* rot_module_offset = new TGeoRotation("rot_module_offset",
1371  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1372  // cout<<"_offset_x: "<<_offset_x<<" cm"<<endl;
1373  TGeoRotation* rot_module_offset = new TGeoRotation("rot_module_offset",
1374  0, 0, 0);
1375  // // ///TEST
1376  // rot_module_offset->RotateZ(-_offset_psi/pi*180.);
1377  // rot_module_offset->RotateY(_offset_theta/pi*180.);
1378  // rot_module_offset->RotateX(-_offset_phi/pi*180.);
1379 
1380  rot_module_offset->RotateZ(_offset_psi / pi * 180.);
1381  rot_module_offset->RotateY(_offset_theta / pi * 180.);
1382  rot_module_offset->RotateX(_offset_phi / pi * 180.);
1383  // rot_module_offset->Print();
1384 
1385  // cout<<"_rot_x_theta: "<<_offset_phi/pi*180.<<" degree"<<endl;
1386  TGeoCombiTrans* rottrans_module_offset = new TGeoCombiTrans(_offset_x,
1387  _offset_y, _offset_z, rot_module_offset);
1388  // rottrans_module = new TGeoHMatrix(*rottrans_module_offset * *rottrans_module);
1389  rottrans_module = new TGeoHMatrix(
1390  *rottrans_module * *rottrans_module_offset);
1391  for (int iside = 0; iside < 2; iside++) { // loop over the two sides of the modules
1392  // rotation around the y axis for the downstream side
1393  _x = 0.;
1394  _y = 0.;
1395  _z = 0.;
1396  if (iside == 0) {
1397  _rotphi = 0.;
1398  _rottheta = 0.;
1399  _rotpsi = 0.;
1400  } else {
1401  _rotphi = 0.;
1402  _rottheta = 180.;
1403  _rotpsi = 0.;
1404  }
1405  TGeoRotation* rot_side = new TGeoRotation("rot_side", _rotphi,
1406  _rottheta, _rotpsi);
1407  TGeoMatrix* rottrans_side = new TGeoCombiTrans(_x, _y, _z, rot_side);
1408 
1409  // Get_offset(ihalf, iplane, imodule, iside, -1, -1,
1410  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1411  // TGeoRotation* rot_side_offset = new TGeoRotation("rot_side_offset",
1412  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1413  // TGeoCombiTrans* rottrans_side_offset =
1414  // new TGeoCombiTrans(_offset_x, _offset_y, _offset_z, rot_side_offset);
1415  // rottrans_side = new TGeoHMatrix(*rottrans_side_offset * *rottrans_side);
1416  // save the transformation from the lumi reference frame
1417  // into the local cvd side reference frame
1418  transformation_matrices[Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1419  new TGeoHMatrix(
1420  (*rottrans_plane) * (*rottrans_module) * (*rottrans_side));
1421  // transformation_matrices_aligned[Generate_Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1422  // new TGeoHMatrix((*rottrans_plane) * (*rottrans_module) * (*rottrans_side)); //TEST
1423 
1424  } // loop over the two sides of the modules
1425  // transformation_matrices[Generate_Tkey(ihalf, iplane, imodule, -1, -1, -1)] =
1426  // new TGeoHMatrix((*rottrans_plane) * (*rottrans_module));
1427  } // loop over modules
1428  } // loop over planes
1429  // // save the transformation into the lumi reference frame
1430  // transformation_matrices[Generate_Tkey(-1, -1, -1, -1, -1, -1)] = new TGeoHMatrix((*lmd_transrot) * (*rottrans_lmd_in_box));
1431  } // loop over detector halves
1432 }
1433 
1434 void PndLmdDim::Correct_transformation_matrices() {
1435  // Read_transformation_matrices(filename_in,true);//read initial matrices
1436 
1437  // ****************************** loops in the luminosity detector ************************
1438  double _offset_x(0), _offset_y(0), _offset_z(0), _offset_phi(0),
1439  _offset_theta(0), _offset_psi(0);
1440  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
1441  //TODO: "Transformation matrix not existent!" =
1442  // TGeoMatrix* rottrans_half = Get_matrix(ihalf, -1, -1 , -1, -1, -1, false);
1443  // Get_offset(ihalf, -1, -1, -1, -1, -1,
1444  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1445  // TGeoRotation* rot_half_offset = new TGeoRotation("rot_half_offset",
1446  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1447  // TGeoCombiTrans* rottrans_half_offset =
1448  // new TGeoCombiTrans(_offset_x, _offset_y, _offset_z, rot_half_offset);
1449  // rottrans_half = new TGeoHMatrix(*rottrans_half_offset * *rottrans_half);
1450 
1451  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
1452  // TODO: "Transformation matrix not existent!" =
1453  // cout<<"Correct matrix plane #"<<iplane<<endl;
1454  // TGeoMatrix* rottrans_plane = Get_matrix(ihalf, iplane, -1 , -1, -1, -1, false);
1455  // Get_offset(ihalf, iplane, -1, -1, -1, -1,
1456  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1457  // TGeoRotation* rot_plane_offset = new TGeoRotation("rot_plane_offset",
1458  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1459  // TGeoCombiTrans* rottrans_plane_offset =
1460  // new TGeoCombiTrans(_offset_x, _offset_y, _offset_z, rot_plane_offset);
1461  // rottrans_plane = new TGeoHMatrix(*rottrans_plane_offset * *rottrans_plane);
1462 
1463  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
1464  //TODO: "Transformation matrix not existent!" =
1465  // cout<<"Correct matrix module #"<<imodule<<" "<<ihalf<<iplane<<imodule<<"-1-1-1"<<endl;
1466  //TGeoMatrix* rottrans_module = Get_matrix(ihalf, iplane, imodule, -1, -1, -1, false);
1467  // Get_offset(ihalf, iplane, imodule, -1, -1, -1,
1468  // _offset_x, _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1469  // cout<<"offsets:"<<_offset_x<<" "<<_offset_y<<" "<<_offset_z<<" "<<_offset_phi<<" "<<_offset_theta<<" "<<_offset_psi<<endl;
1470  // TGeoRotation* rot_module_offset = new TGeoRotation("rot_module_offset",
1471  // _offset_phi/pi*180., _offset_theta/pi*180., _offset_psi/pi*180.);
1472  // TGeoCombiTrans* rottrans_module_offset =
1473  // new TGeoCombiTrans(_offset_x, _offset_y, _offset_z, rot_module_offset);
1474  // rottrans_module = new TGeoHMatrix(*rottrans_module_offset * *rottrans_module);
1475  // // rottrans_module = new TGeoHMatrix(*rottrans_module * *rottrans_module_offset);
1476  for (unsigned int iside = 0; iside < 2; iside++) { // loop over the two sides of the modules
1477  TGeoMatrix* rottrans_side = Get_matrix(ihalf, iplane, imodule, iside,
1478  -1, -1, false);
1479  Get_offset(ihalf, iplane, imodule, iside, -1, -1, _offset_x,
1480  _offset_y, _offset_z, _offset_phi, _offset_theta, _offset_psi);
1481  //cout<<"offsets:"<<_offset_x<<" "<<_offset_y<<" "<<_offset_z<<" "<<_offset_phi<<" "<<_offset_theta<<" "<<_offset_psi<<endl;
1482  TGeoRotation* rot_side_offset = new TGeoRotation("rot_side_offset",
1483  _offset_phi / pi * 180., _offset_theta / pi * 180.,
1484  _offset_psi / pi * 180.);
1485  TGeoCombiTrans* rottrans_side_offset = new TGeoCombiTrans(_offset_x,
1486  _offset_y, _offset_z, rot_side_offset);
1487  rottrans_side = new TGeoHMatrix(
1488  *rottrans_side_offset * *rottrans_side);
1489  // rottrans_side = new TGeoHMatrix(*rottrans_side * *rottrans_side_offset);
1490 
1491  // // save the transformation from the lumi reference frame
1492  // // into the local cvd side reference frame
1493  // transformation_matrices[Generate_Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1494  // new TGeoHMatrix((*rottrans_plane) * (*rottrans_module) * (*rottrans_side));
1495 
1496  // save the transformation from the lumi reference frame
1497  // into the local cvd side reference frame
1498  // transformation_matrices[Generate_Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1499  // new TGeoHMatrix((*rottrans_side)); //TODO: ????
1500  transformation_matrices_aligned[Tkey(ihalf, iplane, imodule, iside,
1501  -1, -1)] = new TGeoHMatrix((*rottrans_side)); //TEST
1502 
1503  } // loop over the two sides of the modules
1504  // transformation_matrices[Generate_Tkey(ihalf, iplane, imodule, -1, -1, -1)] =
1505  // new TGeoHMatrix((*rottrans_plane) * (*rottrans_module)); //TODO: ????
1506  } // loop over modules
1507  } // loop over planes
1508  // // save the transformation into the lumi reference frame
1509  // transformation_matrices[Generate_Tkey(-1, -1, -1, -1, -1, -1)] = new TGeoHMatrix((*lmd_transrot) * (*rottrans_lmd_in_box));
1510  } // loop over detector halves
1511 }
1512 
1513 #include <fstream>
1514 
1515 void PndLmdDim::Read_transformation_matrices(string filename, bool aligned,
1516  int version_number) {
1517  Retrieve_version_number();
1518  map<Tkey, TGeoMatrix*>* matrices = NULL;
1519  if (aligned) {
1520  matrices = &transformation_matrices_aligned;
1521  } else {
1522  matrices = &transformation_matrices;
1523  }
1524  // kill existing matrices
1525  for (it_transformation_matrices = matrices->begin();
1526  it_transformation_matrices != matrices->end();
1527  it_transformation_matrices++) {
1528  delete (it_transformation_matrices->second);
1529  }
1530  matrices->clear();
1531 
1532  string dir = getenv("VMCWORKDIR");
1533  if (filename == ""){
1534  filename = dir + "/geometry/trafo_matrices_lmd.dat";
1535  }
1536  else{
1537  filename = dir + filename;
1538  }
1539 
1540  ifstream file(filename.c_str());
1541  int matrices_counter(0);
1542  if (file.is_open()) {
1543  // read the first line to get additional information
1544  // in case of no (backward compatibility) seek to the first line
1545  string line0;
1546  getline(file, line0);
1547  // first the version info is expected by "version " followed by an integer
1548  string keyword = "version ";
1549  if (line0.compare(0, 8, keyword) != 0) {
1550  // jump to beginning
1551  cout << " **** Warning in PndLmdDim: the input file " << filename
1552  << " has no version! **** " << endl;
1553  file.seekg(0, file.beg);
1554  } else {
1555  int _version = 0;
1556  istringstream inputstream0(&line0[8]);
1557  inputstream0 >> _version;
1558  if (_version != version_number) {
1559  cout << " **** Warning in PndLmdDim: the version " << _version
1560  << " in " << filename << " does not match the geometry version "
1561  << version_number << "! **** " << endl;
1562  }
1563  }
1564  while (1) {
1565  // read the key
1566  string key;
1567  getline(file, key);
1568  //cout << "next key " << key << endl;
1569  //char newline;
1570  if (file.good()) {
1571  // generate rotation and translation matrices with it
1572  // error treatment is not foreseen yet
1573  string line1;
1574  getline(file, line1);
1575  istringstream inputstream1(line1);
1576  double translation[3];
1577  for (int i = 0; i < 3; i++) {
1578  inputstream1 >> translation[i];
1579  }
1580  string line2;
1581  getline(file, line2);
1582  istringstream inputstream2(line2);
1583  double rotation[9];
1584  for (int i = 0; i < 9; i++) {
1585  inputstream2 >> rotation[i];
1586  }
1587  TGeoRotation* rot = new TGeoRotation("rotation");
1588  rot->SetMatrix(rotation);
1589  TGeoCombiTrans* rottrans = new TGeoCombiTrans(translation[0],
1590  translation[1], translation[2], rot);
1591  (*matrices)[Tkey(key)] = rottrans;
1592  matrices_counter++;
1593  } else {
1594  file.close();
1595  break;
1596  }
1597  }
1598  file.close();
1599  cout << " Read " << matrices_counter << " matrices from " << filename
1600  << endl;
1601  } else {
1602  cout
1603  << " Error in PndLmdDim::Read_transformation_matrices: could not read from "
1604  << filename << endl;
1605  }
1606 }
1607 
1608 TGeoHMatrix* PndLmdDim::Get_matrix(string path, bool aligned, int ihalf,
1609  int iplane, int imodule, int iside, int idie, int isensor) {
1610  TGeoHMatrix* result = NULL;
1611  if (!gGeoManager->CheckPath(path.c_str())) {
1612  cout << " Error in PndLmdDim::Get_matrix: path " << path << " is not valid "
1613  << endl;
1614  return result;
1615  }
1616  if (!aligned) { // try to obtain the original matrix, if available
1617  TGeoPNEntry *pne;
1618  //TGeoHMatrix *ide; // ideal
1619  //TGeoHMatrix *mis; // misaligned
1620  string uname = Generate_key(ihalf, iplane, imodule, iside, idie, isensor);
1621  pne = gGeoManager->GetAlignableEntry(uname.c_str());
1622  if (pne) {
1623  TGeoPhysicalNode *node = pne->GetPhysicalNode();// new TGeoPhysicalNode(path.c_str());//
1624  if (node) {
1625  result = new TGeoHMatrix(*node->GetOriginalMatrix());
1626  //if (aligned) result = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
1627  return result;
1628  } else {
1629  cout << " no node found in pn entry at " << path << endl;
1630  cout << " obtaining default matrix " << endl;
1631  }
1632  } else {
1633  cout << " no pn entry (alignable node) found in " << path << endl;
1634  cout << " obtaining default matrix " << endl;
1635  }
1636  }
1637  // I don't care if it was aligned or not, use the last one
1638  gGeoManager->cd(path.c_str());
1639  TGeoMatrix* matrix = gGeoManager->GetCurrentNode()->GetMatrix();
1640  if (!matrix)
1641  return NULL;
1642  result = new TGeoHMatrix(*(matrix));
1643  //result->Print();
1644  //if (result){
1645  // cout << " found a transformation matrix for " << path << endl;
1646  // result->Print();
1647  // return result;
1648  //}
1649 
1650  //
1651  //if (!pne){
1652  // cout << " no pn entry found in " << path << endl;
1653  // return 0;
1654  //}
1655  //TGeoPhysicalNode *node = new TGeoPhysicalNode(path.c_str());//pne->GetPhysicalNode();
1656  //if (!node){
1657  // cout << " no node found in pn entry at " << path << endl;
1658  // return 0;
1659  //}
1660  //result = new TGeoHMatrix(*node->GetOriginalMatrix());
1661  //if (aligned) result = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
1662  //else result = new TGeoHMatrix(*node->GetOriginalMatrix());
1663  //if (!result){
1664  // cout << " warning: no matrix found for the path " << path << endl;
1665  //}
1666  //delete node;
1667  return result;
1668 }
1669 
1670 bool PndLmdDim::Set_matrix(string path, TGeoHMatrix* matrix, int ihalf,
1671  int iplane, int imodule, int iside, int idie, int isensor) {
1672  if (!matrix)
1673  return false;
1674  if (!gGeoManager->CheckPath(path.c_str())) {
1675  cout << " Error in PndLmdDim::Set_matrix: path " << path << " is not valid "
1676  << endl;
1677  return false;
1678  }
1679  TGeoPNEntry *pne;
1680  TGeoPhysicalNode *node;
1681  string uname = Generate_key(ihalf, iplane, imodule, iside, idie, isensor);
1682  pne = gGeoManager->GetAlignableEntry(uname.c_str());
1683  if (!pne) {
1684  cout << " creating alignable entry for " << path << endl;
1685  if (isensor > 0) {
1686  int sensID = Get_sensor_id(ihalf, iplane, imodule, iside, idie, isensor);
1687  pne = gGeoManager->SetAlignableEntry(uname.c_str(), path.c_str(), sensID);
1688  } else {
1689  pne = gGeoManager->SetAlignableEntry(uname.c_str(), path.c_str());
1690  }
1691  if (pne) {
1692  node = new TGeoPhysicalNode(path.c_str());
1693  pne->SetPhysicalNode(node);
1694  }
1695  }
1696  if (!pne) {
1697  cout << " Error: no pn entry (alignable node) at " << path << " created "
1698  << endl;
1699  return false;
1700  }
1701  node = pne->GetPhysicalNode();
1702  if (!node) {
1703  cout << " no node found for pn entry (alignable node) at " << path << endl;
1704  return false;
1705  }
1706  node->Align(matrix); //GetNode(node->GetLevel())->Align();
1707  //cout << " matrix at " << path << " successfully aligned " << endl;
1708  /*
1709  if (aligned) return new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
1710  else return new TGeoHMatrix(*node->GetOriginalMatrix()); */
1711  return true;
1712 }
1713 
1714 bool PndLmdDim::Read_transformation_matrices_from_geometry(bool aligned) {
1715  if (!Retrieve_version_number()) {
1716  return false;
1717  }
1718  if (geometry_version < 3) {
1719  cout
1720  << " *** Error in PndLmdDim::Read_transformation_matrices_from_geometry:"
1721  << endl;
1722  cout << " geometry version " << geometry_version
1723  << " is not compatible with this method! " << endl;
1724  return false;
1725  }
1726  bool result = false;
1727  //FairRootManager* ioman = FairRootManager::Instance();
1728  TGeoManager* gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
1729  if (!gGeoMan) {
1730  cout << " Info: no FAIRGeom found, using gGeoManager " << endl;
1731  gGeoMan = gGeoManager;
1732  }
1733  if (gGeoMan) {
1734  vector<string> list_of_sensors;
1735  string path_to_top = Get_List_of_Sensors(list_of_sensors);
1736  int offset;
1737  if (!Test_List_of_Sensors(list_of_sensors, offset)) {
1738  cout
1739  << " *** Error in PndLmdDim::Read_transformation_matrices_from_geometry:"
1740  << endl;
1741  cout << " could not retrieve list of sensors from geometry " << endl;
1742  return false;
1743  }
1744 
1745  map<Tkey, TGeoMatrix*>* matrices = NULL;
1746  if (aligned) {
1747  matrices = &transformation_matrices_aligned;
1748  } else {
1749  matrices = &transformation_matrices;
1750  }
1751  // kill existing matrices
1752  for (it_transformation_matrices = matrices->begin();
1753  it_transformation_matrices != matrices->end();
1754  it_transformation_matrices++) {
1755  delete (it_transformation_matrices->second);
1756  }
1757  matrices->clear();
1758 
1759  //int matrices_counter(0); //[R.K. 01/2017] unused variable?
1760  //TGeoHMatrix* matrix; //[R.K. 01/2017] unused variable?
1761 
1762  TGeoVolume* vol = gGeoMan->FindVolumeFast(nav_paths[0].c_str());
1763  if (vol) { // check only the existence of the top node, everything else will be not checked, so earlier simulations will fail here
1764  /*
1765  string sversion = vol->GetTitle();
1766  if (sversion.compare(0,8,"version ")!=0){
1767  cout << " Warning from PndLmdDim::Retrieve_version_number: no version number encoded in the node title. Setting it to 0. " << endl;
1768  geometry_version = 0;
1769  } else {
1770  // take the number behind "version " string
1771  geometry_version = atoi(&(vol->GetTitle()[8]));
1772  }
1773  cout << " Info from PndLmdDim::Retrieve_version_number: The geometry version was set to " << geometry_version << endl;
1774  //cout << vol->GetName() << endl;*/
1775  int sensor_id = offset;
1776  stringstream path_to_ref;
1777  path_to_ref << path_to_top << "/" << nav_paths[1] << "_0";
1778  gGeoMan->cd(path_to_ref.str().c_str());
1779  TGeoHMatrix* top_matrix = new TGeoHMatrix(*gGeoMan->GetCurrentMatrix());
1780  //top_matrix->Print();
1781  if (top_matrix)
1782  (*matrices)[Tkey(-1, -1, -1, -1, -1, -1)] = new TGeoHMatrix(
1783  *top_matrix);
1784  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
1785  stringstream path_to_half;
1786  path_to_half << path_to_ref.str() << "/" << nav_paths[2] << "_"
1787  << ihalf;
1788  //gGeoMan->cd(path_to_half.str().c_str());
1789  TGeoHMatrix* half_matrix = Get_matrix(path_to_half.str(), aligned,
1790  ihalf, -1, -1, -1, -1, -1); //gGeoMan->GetCurrentNode()->GetMatrix();
1791  // here is some mess, sorry for that!
1792  //if (half_matrix && top_matrix) (*matrices)[Tkey(-1, -1, -1, -1, -1, -1)] =
1793  // new TGeoHMatrix(*half_matrix * *top_matrix);
1794  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
1795  stringstream path_to_plane;
1796  path_to_plane << path_to_half.str() << "/" << nav_paths[3] << "_"
1797  << iplane;
1798  //gGeoMan->cd(path_to_plane.str().c_str());
1799  TGeoHMatrix* plane_matrix = Get_matrix(path_to_plane.str(), aligned,
1800  ihalf, iplane, -1, -1, -1, -1);
1801  //if (matrix) (*matrices)[Tkey(ihalf, iplane, imodule, iside, idie, isensor)] = new TGeoHMatrix(*matrix);
1802  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
1803  stringstream path_to_module;
1804  path_to_module << path_to_plane.str() << "/" << nav_paths[4] << "_"
1805  << imodule;
1806  //gGeoMan->cd(path_to_module.str().c_str());
1807  TGeoHMatrix* module_matrix = Get_matrix(path_to_module.str(),
1808  aligned, ihalf, iplane, imodule, -1, -1, -1);
1809  //if (matrix) (*matrices)[Tkey(ihalf, iplane, imodule, iside, idie, isensor)] = new TGeoHMatrix(*matrix);
1810  for (unsigned int iside = 0; iside < 2; iside++) { // loop over the two sides of the modules
1811  stringstream path_to_side;
1812  path_to_side << path_to_module.str() << "/" << nav_paths[5] << "_"
1813  << iside;
1814  //gGeoMan->cd(path_to_side.str().c_str());
1815  TGeoHMatrix* side_matrix = Get_matrix(path_to_side.str(), aligned,
1816  ihalf, iplane, imodule, iside, -1, -1); //gGeoMan->GetCurrentNode()->GetMatrix();
1817  if (plane_matrix && module_matrix && side_matrix)
1818  (*matrices)[Tkey(ihalf, iplane, imodule, iside, -1, -1)] =
1819  new TGeoHMatrix(
1820  *plane_matrix * *module_matrix * *side_matrix);
1821  for (unsigned int idie = 0; idie < 2; idie++) { // loop over dies
1822  stringstream path_to_die;
1823  path_to_die << path_to_side.str() << "/" << nav_paths[6] << "_"
1824  << idie;
1825  TGeoHMatrix* die_matrix = Get_matrix(path_to_die.str(), aligned,
1826  ihalf, iplane, imodule, iside, idie, -1);
1827  for (unsigned int isensor = 0; isensor < 3; isensor++) { // loop over sensors
1828  if (idie == 1) {
1829  if (isensor == 0)
1830  continue;
1831  }
1832  stringstream path_to_sensor;
1833  path_to_sensor << path_to_die.str() << "/" << nav_paths[7]
1834  << "_" << sensor_id;
1835  TGeoHMatrix* sensor_matrix = Get_matrix(path_to_sensor.str(),
1836  aligned, ihalf, iplane, imodule, iside, idie, isensor); //gGeoMan->GetCurrentNode()->GetMatrix();
1837  if (sensor_matrix && die_matrix)
1838  (*matrices)[Tkey(ihalf, iplane, imodule, iside, idie,
1839  isensor)] = new TGeoHMatrix(
1840  *die_matrix * *sensor_matrix); // do I have to delete it, or should I copy it?
1841  //new TGeoHMatrix((*rottrans_die) * (*rottrans_sensor));
1842  sensor_id++;
1843  delete sensor_matrix;
1844  } // loop over sensors
1845  delete die_matrix;
1846  } // loop over dies
1847  delete side_matrix;
1848  } // loop over the two sides of the modules
1849  delete module_matrix;
1850  } // loop over modules
1851  delete plane_matrix;
1852  } // loop over planes
1853  delete half_matrix;
1854  } // loop over detector halves
1855  delete top_matrix;
1856  gGeoManager->RefreshPhysicalNodes(false);
1857  cout << " read " << (*matrices).size() << " matrices from root geometry "
1858  << endl;
1859  result = true;
1860  } else {
1861  cout
1862  << " *** Error in PndLmdDim::Read_transformation_matrices_from_geometry:"
1863  << endl;
1864  cout << " Could not find the top volume " << nav_paths[0].c_str()
1865  << " to retrieve the transformation matrix for the luminosity detector! Is the geometry already loaded? "
1866  << endl;
1867  }
1868  } else {
1869  cout
1870  << " *** Error in PndLmdDim::Read_transformation_matrices_from_geometry:"
1871  << endl;
1872  cout
1873  << " Could not find a GeoManager to load the luminosity detector matrices from it! "
1874  << endl;
1875  }
1876  return result;
1877 }
1878 
1879 bool PndLmdDim::Write_transformation_matrices_to_geometry(bool aligned) {
1880  if (!Retrieve_version_number()) {
1881  return false;
1882  }
1883  if (geometry_version < 3) {
1884  cout
1885  << " *** Error in PndLmdDim::Write_transformation_matrices_to_geometry:"
1886  << endl;
1887  cout << " geometry version " << geometry_version
1888  << " is not compatible with this method! " << endl;
1889  return false;
1890  }
1891  bool result = false;
1892  //FairRootManager* ioman = FairRootManager::Instance();
1893  TGeoManager* gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
1894  if (!gGeoMan) {
1895  cout << " Info: no FAIRGeom found, using gGeoManager " << endl;
1896  gGeoMan = gGeoManager;
1897  }
1898  if (gGeoMan) {
1899  vector<string> list_of_sensors;
1900  string path_to_top = Get_List_of_Sensors(list_of_sensors);
1901  int offset;
1902  if (!Test_List_of_Sensors(list_of_sensors, offset)) {
1903  cout
1904  << " *** Error in PndLmdDim::Write_transformation_matrices_to_geometry:"
1905  << endl;
1906  cout << " could not retrieve list of sensors from geometry " << endl;
1907  return false;
1908  }
1909 
1910  map<Tkey, TGeoMatrix*>* matrices = NULL;
1911  if (aligned) {
1912  matrices = &transformation_matrices_aligned;
1913  } else {
1914  matrices = &transformation_matrices;
1915  }
1916  // check matrix existence
1917  if (matrices->size() == 0) {
1918  cout
1919  << " *** Error in PndLmdDim::Write_transformation_matrices_to_geometry:"
1920  << endl;
1921  cout << " no matrices to apply to the geometry! " << endl;
1922  return false;
1923  }
1924 
1925  //int matrices_counter(0); //[R.K. 01/2017] unused variable?
1926  TGeoHMatrix* matrix;
1927  TGeoHMatrix* matrix_mother;
1928 
1929  TGeoVolume* vol = gGeoMan->FindVolumeFast(nav_paths[0].c_str());
1930  if (vol) { // check only the existence of the top node, everything else will be not checked, so earlier simulations will fail here
1931  /*
1932  string sversion = vol->GetTitle();
1933  if (sversion.compare(0,8,"version ")!=0){
1934  cout << " Warning from PndLmdDim::Retrieve_version_number: no version number encoded in the node title. Setting it to 0. " << endl;
1935  geometry_version = 0;
1936  } else {
1937  // take the number behind "version " string
1938  geometry_version = atoi(&(vol->GetTitle()[8]));
1939  }
1940  cout << " Info from PndLmdDim::Retrieve_version_number: The geometry version was set to " << geometry_version << endl;
1941  //cout << vol->GetName() << endl;*/
1942  int sensor_id = offset;
1943  stringstream path_to_ref;
1944  path_to_ref << path_to_top << "/" << nav_paths[1] << "_0";
1945  matrix = (TGeoHMatrix*) Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
1946  // obtain matrix up to the mother volume to subtract it
1947  // from the total one in the search tree
1948  gGeoMan->cd(path_to_top.c_str());
1949  matrix_mother = (TGeoHMatrix*) gGeoMan->GetCurrentMatrix();
1950  if (matrix && matrix_mother) {
1951  TGeoHMatrix* _matrix = new TGeoHMatrix(
1952  matrix_mother->Inverse() * *matrix);
1953  Set_matrix(path_to_ref.str(), _matrix, -1, -1, -1, -1, -1, -1);
1954  }
1955  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
1956  stringstream path_to_half;
1957  path_to_half << path_to_ref.str() << "/" << nav_paths[2] << "_"
1958  << ihalf;
1959  gGeoMan->cd(path_to_half.str().c_str());
1960  // get the original matrix to subtract it
1961  // I'm not sure if it is always correct
1962  //TGeoHMatrix* half_matrix = Get_matrix(path_to_half.str(), false, ihalf,
1963  //-1, -1, -1, -1, -1); //[R.K. 01/2017] unused variable?
1964  //matrix = (TGeoHMatrix*) Get_matrix(ihalf,-1,-1,-1,-1,-1, aligned);
1965  //if (matrix) Set_matrix(path_to_half.str(), matrix,
1966  // ihalf, -1, -1, -1, -1, -1);
1967  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
1968  stringstream path_to_plane;
1969  path_to_plane << path_to_half.str() << "/" << nav_paths[3] << "_"
1970  << iplane;
1971  //gGeoMan->cd(path_to_plane.str().c_str());
1972  TGeoHMatrix* plane_matrix = Get_matrix(path_to_plane.str(), false,
1973  ihalf, iplane, -1, -1, -1, -1);
1974  //matrix = (TGeoHMatrix*) Get_matrix(ihalf,iplane,-1,-1,-1,-1, aligned);
1975  //if (matrix) Set_matrix(path_to_plane.str(), matrix,
1976  // ihalf, iplane, -1, -1, -1, -1);
1977  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
1978  stringstream path_to_module;
1979  path_to_module << path_to_plane.str() << "/" << nav_paths[4] << "_"
1980  << imodule;
1981  TGeoHMatrix* module_matrix = Get_matrix(path_to_module.str(), false,
1982  ihalf, iplane, imodule, -1, -1, -1);
1983  //gGeoMan->cd(path_to_module.str().c_str());
1984  //matrix = (TGeoHMatrix*) Get_matrix(ihalf,iplane,imodule,-1,-1,-1, aligned);
1985  //if (matrix) Set_matrix(path_to_module.str(), matrix,
1986  // ihalf, iplane, imodule, -1, -1, -1);
1987  for (unsigned int iside = 0; iside < 2; iside++) { // loop over the two sides of the modules
1988  stringstream path_to_side;
1989  path_to_side << path_to_module.str() << "/" << nav_paths[5] << "_"
1990  << iside;
1991  //gGeoMan->cd(path_to_side.str().c_str());
1992  TGeoHMatrix* side_matrix = Get_matrix(path_to_side.str(), false,
1993  ihalf, iplane, imodule, iside, -1, -1);
1994  matrix = (TGeoHMatrix*) Get_matrix(ihalf, iplane, imodule, iside,
1995  -1, -1, aligned);
1996  if (matrix && plane_matrix && module_matrix && side_matrix) {
1997  TGeoHMatrix* _matrix = new TGeoHMatrix(
1998  (*plane_matrix * *module_matrix).Inverse() * *matrix);
1999  Set_matrix(path_to_side.str(), _matrix, ihalf, iplane, imodule,
2000  iside, -1, -1);
2001  }
2002  for (unsigned int idie = 0; idie < 2; idie++) { // loop over dies
2003  stringstream path_to_die;
2004  path_to_die << path_to_side.str() << "/" << nav_paths[6] << "_"
2005  << idie;
2006  //gGeoMan->cd(path_to_die.str().c_str());
2007  //matrix = (TGeoHMatrix*) Get_matrix(ihalf,iplane,imodule,iside,idie,-1, aligned);
2008  TGeoHMatrix* die_matrix = Get_matrix(path_to_die.str(), aligned,
2009  ihalf, iplane, imodule, iside, idie, -1);
2010  //if (matrix) Set_matrix(path_to_die.str(), matrix,
2011  // ihalf, iplane, imodule, iside, idie, -1);
2012  for (unsigned int isensor = 0; isensor < 3; isensor++) { // loop over sensors
2013  if (idie == 1) {
2014  if (isensor == 0)
2015  continue;
2016  }
2017  stringstream path_to_sensor;
2018  stringstream path_to_sensor_passive;
2019  path_to_sensor << path_to_die.str() << "/" << nav_paths[7]
2020  << "_" << sensor_id;
2021  // here the passive part of the sensor must be also shifted
2022  // this was not reflected in the original idea to store the path in the nav_paths
2023  path_to_sensor_passive << path_to_die.str() << "/"
2024  << "LumPassiveRect" << "_" << sensor_id - sensIDoffset;
2025  //gGeoMan->cd(path_to_sensor.str().c_str());
2026  TGeoHMatrix* sensor_matrix = Get_matrix(path_to_sensor.str(),
2027  aligned, ihalf, iplane, imodule, iside, idie, isensor);
2028  matrix = (TGeoHMatrix*) Get_matrix(ihalf, iplane, imodule,
2029  iside, idie, isensor, aligned);
2030  if (sensor_matrix && die_matrix && matrix) {
2031  TGeoHMatrix* _matrix = new TGeoHMatrix(
2032  die_matrix->Inverse() * *matrix);
2033  //TGeoHMatrix* _matrix_passive = new TGeoHMatrix(die_matrix->Inverse() * *matrix);
2034  Set_matrix(path_to_sensor.str(), _matrix, ihalf, iplane,
2035  imodule, iside, idie, isensor);
2036  Set_matrix(path_to_sensor_passive.str(), _matrix, ihalf,
2037  iplane, imodule, iside, idie, isensor + 1000);
2038  }
2039  sensor_id++;
2040  } // loop over sensors
2041  //lmd_vol_side_->AddNode(lmd_vol_die_, 0, rottrans_die);
2042  } // loop over dies
2043  } // loop over the two sides of the modules
2044  //lmd_vol_plane_->AddNode(lmd_vol_module_, 0, rottrans_module);
2045  } // loop over modules
2046  } // loop over planes
2047  //lmd_vol_ref_sys->AddNode(lmd_vol_half_, 0, rottrans_no);
2048  //transformation_matrices[Tkey(-1, -1, -1, -1, -1, -1)] = new TGeoHMatrix((*lmd_transrot) * (*rottrans_lmd_in_box));
2049  } // loop over detector halves
2050  cout << " set " << (*matrices).size() << " matrices to root geometry "
2051  << endl;
2052  result = true;
2053  } else {
2054  cout
2055  << " *** Error in PndLmdDim::Write_transformation_matrices_to_geometry:"
2056  << endl;
2057  cout << " Could not find the top volume " << nav_paths[0].c_str()
2058  << " to retrieve the transformation matrix for the luminosity detector! Is the geometry already loaded? "
2059  << endl;
2060  }
2061  gGeoMan->RefreshPhysicalNodes();
2062  } else {
2063  cout
2064  << " *** Error in PndLmdDim::Write_transformation_matrices_to_geometry:"
2065  << endl;
2066  cout
2067  << " Could not find a GeoManager to load the luminosity detector matrices from it! "
2068  << endl;
2069  }
2070 
2071  return result;
2072 }
2073 
2074 string PndLmdDim::Get_List_of_Sensors(vector<string>& list_of_sensors,
2075  bool found_lmd, bool first_call) {
2076  string result("");
2077  if (first_call) {
2078  gGeoManager->CdTop();
2079  first_call = false;
2080  }
2081  //cout << gGeoManager->GetPath() << " ";
2082  string nodename(
2083  gGeoManager->GetCurrentNavigator()->GetCurrentNode()->GetName());
2084  //cout << "current node is " << nodename << endl;
2085  if (nodename.compare(0, nav_paths[7].size(), nav_paths[7]) == 0) {
2086  //cout << " found a sensor " << nav_paths[7] << endl;
2087  list_of_sensors.push_back(gGeoManager->GetPath());
2088  }
2089  int nnodes = gGeoManager->GetCurrentNode()->GetNdaughters();
2090  for (int i = 0; i < nnodes; i++) {
2091  //cout << " navigating into node " << i << endl;
2092  gGeoManager->CdDown(i);
2093  result = Get_List_of_Sensors(list_of_sensors, found_lmd, first_call);
2094  if (nodename == nav_paths[1] || found_lmd) {
2095  cout << " found the lmd node! Aborting recursive search. " << endl;
2096  gGeoManager->CdUp();
2097  break;
2098  }
2099  gGeoManager->CdUp();
2100  }
2101  if (nodename.compare(0, nav_paths[0].size(), nav_paths[0]) == 0) {
2102  result = gGeoManager->GetPath();
2103  cout << " top volume is " << result << endl;
2104  }
2105  return result;
2106 }
2107 
2108 bool PndLmdDim::Test_List_of_Sensors(vector<string> list_of_sensors,
2109  int& offset) {
2110  bool result = true;
2111  offset = 0;
2112  int nsensorstotal = n_sensors * 2 * nmodules * 2 * n_planes;
2113  int max_sensID = nsensorstotal - 1;
2114  int checksum = 0;
2115  int sumtocheck = 0;
2116  if ((int)list_of_sensors.size() != nsensorstotal) {
2117  cout << " PndLmdDim::Test_List_of_Sensors: number of sensors is wrong "
2118  << list_of_sensors.size() << " != " << nsensorstotal << endl;
2119  return false;
2120  }
2121  for (int isensor = 0; isensor < nsensorstotal; isensor++) {
2122  gGeoManager->cd(list_of_sensors[isensor].c_str());
2123  string path(gGeoManager->GetCurrentNavigator()->GetPath());
2124  if (path != list_of_sensors[isensor]) {
2125  cout << " PndLmdDim::Test_List_of_Sensors: Could not navigate to "
2126  << list_of_sensors[isensor] << endl;
2127  return false;
2128  }
2129  // strip the copy number in the path
2130  char sensID_char[5] = { '0', '0', '0', '\0' };
2131  sensID_char[4] = '\0';
2132  for (unsigned ichar = 3; ichar > 0; ichar--) { //path.size()-1; ichar > path.size()-5; ichar--){
2133  char digit = path[path.size() - 4 + ichar];
2134  if (digit == '_')
2135  break;
2136  sensID_char[ichar] = digit;
2137  }
2138  int sensID = atoi(sensID_char);
2139  if (sensID > max_sensID) {
2140  max_sensID = sensID;
2141  }
2142  checksum += isensor;
2143  sumtocheck += sensID;
2144  }
2145  offset = max_sensID - (nsensorstotal - 1);
2146  if (offset != 0) {
2147  cout
2148  << " PndLmdDim::Test_List_of_Sensors: Found an offset in the sensorIDs of "
2149  << offset << endl;
2150  checksum += (nsensorstotal * offset);
2151  }
2152  if (checksum != sumtocheck) {
2153  cout
2154  << " PndLmdDim::Test_List_of_Sensors: check sum of sensor id's does not match "
2155  << checksum << " != " << sumtocheck << endl;
2156  return false;
2157  }
2158  return result;
2159 }
2160 
2161 bool PndLmdDim::Set_sensIDoffset(int offset) {
2162  if (offset >= 0) {
2163  sensIDoffset = (int) offset;
2164  return true;
2165  }
2166  vector<string> sensors;
2167  Get_List_of_Sensors(sensors);
2168  if (Test_List_of_Sensors(sensors, offset)) {
2169  sensIDoffset = offset;
2170  return true;
2171  } else {
2172  cout << " Could not set sensIDoffset from root geometry " << endl;
2173  return false;
2174  }
2175 }
2176 
2177 #include <iomanip>
2178 
2179 void PndLmdDim::Write_transformation_matrices(string filename, bool aligned,
2180  int version_number) {
2181  map<Tkey, TGeoMatrix*>* matrices = NULL;
2182  if (aligned) {
2183  matrices = &transformation_matrices_aligned;
2184  } else {
2185  matrices = &transformation_matrices;
2186  }
2187  int matrices_counter(0);
2188  string dir = getenv("VMCWORKDIR");
2189  if (filename == "")
2190  filename = dir + "/geometry/trafo_matrices_lmd.dat";
2191  ofstream file(filename.c_str());
2192  if (file.is_open()) {
2193  // write the first line to add additional information
2194  file << "version " << version_number << "\n";
2195  for (it_transformation_matrices = matrices->begin();
2196  it_transformation_matrices != matrices->end();
2197  it_transformation_matrices++) {
2198  // write the key
2199  Tkey key = it_transformation_matrices->first;
2200  file
2201  << Generate_key(key.half, key.plane, key.module, key.side, key.die,
2202  key.sensor) << '\n';
2203  const double * translation =
2204  it_transformation_matrices->second->GetTranslation();
2205  const double * rotation =
2206  it_transformation_matrices->second->GetRotationMatrix();
2207  // write the numbers for the translation
2208  for (int i = 0; i < 3; i++) {
2209  file << setw(14) << setprecision(14) << scientific << translation[i];
2210  if (i == 2)
2211  file << '\n';
2212  else
2213  file << ' ';
2214  }
2215  // write the numbers for the rotation
2216  for (int i = 0; i < 9; i++) {
2217  file << setw(14) << setprecision(14) << scientific << rotation[i];
2218  if (i == 8)
2219  file << '\n';
2220  else
2221  file << ' ';
2222  }
2223  matrices_counter++;
2224  }
2225  file.close();
2226  cout << matrices_counter << " Transformation matrices were written to "
2227  << filename << endl;
2228  } else
2229  cout
2230  << " Error in PndLmdDim::Write_transformation_matrices: could not write to "
2231  << filename << endl;
2232 }
2233 
2234 void PndLmdDim::Cleanup() {
2235  for (it_transformation_matrices = transformation_matrices.begin();
2236  it_transformation_matrices != transformation_matrices.end();
2237  it_transformation_matrices++) {
2238  delete (it_transformation_matrices->second);
2239  }
2240  transformation_matrices.clear();
2241  for (it_transformation_matrices = transformation_matrices_aligned.begin();
2242  it_transformation_matrices != transformation_matrices_aligned.end();
2243  it_transformation_matrices++) {
2244  delete (it_transformation_matrices->second);
2245  }
2246  transformation_matrices_aligned.clear();
2247 }
2248 
2249 void PndLmdDim::Read_DB_offsets(PndLmdAlignPar *lmdalignpar) {
2250  // ana = FairRun::Instance();
2251  // rtdb=ana->GetRuntimeDb();
2252  // TList* fAlignParamList = new TList();// alignment params fom DB
2253  double fShiftX[40], fShiftY[40], fShiftZ[40];
2254  double fRotateX[40], fRotateY[40], fRotateZ[40];
2255  // //read params for lumi alignment
2256  // PndLmdContFact* thelmdcontfact = (PndLmdContFact*)rtdb->getContFactory("PndLmdContFact");
2257  // TList* theAlignLMDContNames = thelmdcontfact->GetAlignParNames();
2258  // Info("SetParContainers()","AlignLMD The container names list contains %i entries",theAlignLMDContNames->GetEntries());
2259  // TIter cfAlIter(theAlignLMDContNames);
2260  // while (TObjString* contname = (TObjString*)cfAlIter()) {
2261  // TString parsetname = contname->String();
2262  // Info("SetParContainers()",parsetname.Data());
2263  // PndLmdAlignPar *lmdalignpar = (PndLmdAlignPar*)(rtdb->getContainer(parsetname.Data()));
2264  // if(!lmdalignpar) Fatal("SetParContainers","No ALIGN parameter found: %s",parsetname.Data());
2265  // fAlignParamList->Add(lmdalignpar);
2266  // }
2267 
2268  // TIter alignparams(fAlignParamList);
2269  // PndLmdAlignPar* lmdalignpar=(PndLmdAlignPar*)alignparams();
2270  // // lmdalignpar->getParams(alignparams);
2271  // lmdalignpar->Print();
2272  if (0 == lmdalignpar) {
2273  Error("PndLmdStripClusterTask::SetCalculators()",
2274  "A ALIGN Parameter Set does not exist properly.");
2275  } else {
2276  // cout<<"&&&&& it's from LmdDim &&&&&"<<endl;
2277  // lmdalignpar->Print();
2278  for (int ik = 0; ik < 40; ik++) {
2279  // cout<<"ik = "<<ik<<endl;
2280  fShiftX[ik] = lmdalignpar->GetShiftX(ik);
2281  // cout<<"fShiftX[ik] = "<<fShiftX[ik]<<endl;
2282  fShiftY[ik] = lmdalignpar->GetShiftY(ik);
2283  fShiftZ[ik] = lmdalignpar->GetShiftZ(ik);
2284  fRotateX[ik] = lmdalignpar->GetRotateX(ik);
2285  fRotateY[ik] = lmdalignpar->GetRotateY(ik);
2286  fRotateZ[ik] = lmdalignpar->GetRotateZ(ik);
2287  // if (fVerbose > 2) cout<<"fShiftX["<<ik<<"]="<<fShiftX[ik]<<" fRotateX["<<ik<<"]="<<fRotateX[ik]
2288  // <<" fRotateY["<<ik<<"]="<<fRotateY[ik]<<" fRotateZ["<<ik<<"]="<<fRotateZ[ik]<<endl;
2289  }
2290  }
2291  // lmdalignpar->Print();
2292 
2293  //so far for modules alignment only
2294  //TODO: individual sensor and\or die alignment???
2295  for (unsigned int ihalf = 0; ihalf < 2; ihalf++) { // loop over detector halves
2296  for (unsigned int iplane = 0; iplane < n_planes; iplane++) { // loop over planes
2297  for (unsigned int imodule = 0; imodule < nmodules; imodule++) { // loop over modules
2298  Tkey key(ihalf, iplane, imodule, -1, -1, -1);
2299  int ikey = (ihalf * n_planes * nmodules) + (iplane * nmodules)
2300  + imodule;
2301  // cout<<"for: "<<ihalf<<iplane<<imodule<<": ikey="<<ikey<<endl;
2302  offsets[key].push_back(fShiftX[ikey]);
2303  offsets[key].push_back(fShiftY[ikey]);
2304  offsets[key].push_back(fShiftZ[ikey]);
2305  offsets[key].push_back(fRotateX[ikey]);
2306  offsets[key].push_back(fRotateY[ikey]);
2307  offsets[key].push_back(fRotateZ[ikey]);
2308 
2309  Tkey key2(ihalf, iplane, imodule, 0, -1, -1); //top side
2310  // cout<<"for: "<<ihalf<<iplane<<imodule<<": ikey="<<ikey<<endl;
2311  offsets[key2].push_back(fShiftX[ikey]);
2312  offsets[key2].push_back(fShiftY[ikey]);
2313  offsets[key2].push_back(fShiftZ[ikey]);
2314  offsets[key2].push_back(fRotateX[ikey]);
2315  offsets[key2].push_back(fRotateY[ikey]);
2316  offsets[key2].push_back(fRotateZ[ikey]);
2317 
2318  // //TODO: is it correct???
2319  Tkey key3(ihalf, iplane, imodule, 1, -1, -1); // bottom side
2320  offsets[key3].push_back(fShiftX[ikey]);
2321  offsets[key3].push_back(fShiftY[ikey]);
2322  offsets[key3].push_back(fShiftZ[ikey]);
2323  offsets[key3].push_back(fRotateX[ikey]);
2324  offsets[key3].push_back(fRotateY[ikey]);
2325  offsets[key3].push_back(fRotateZ[ikey]);
2326  }
2327  }
2328  }
2329 }
2330 
2331 void PndLmdDim::Set_offset(int ihalf, int iplane, int imodule, int iside,
2332  int idie, int isensor, double x, double y, double z, double rotphi,
2333  double rottheta, double rotpsi) {
2334  Tkey key = Tkey(ihalf, iplane, imodule, iside, idie, isensor);
2335  itoffset = offsets.find(key);
2336  if (itoffset != offsets.end()) {
2337  cout
2338  << " **** Warning in PndLmdDim::Set_offset: offset exists already! Replacing it! *** "
2339  << endl;
2340  } else {
2341  offsets[key].clear();
2342  offsets[key].push_back(x);
2343  offsets[key].push_back(y);
2344  offsets[key].push_back(z);
2345  offsets[key].push_back(rotphi);
2346  offsets[key].push_back(rottheta);
2347  offsets[key].push_back(rotpsi);
2348  }
2349 
2350 }
2351 
2352 void PndLmdDim::Get_offset(int ihalf, int iplane, int imodule, int iside,
2353  int idie, int isensor, double& x, double& y, double& z, double& rotphi,
2354  double& rottheta, double& rotpsi, bool random) {
2355  Tkey key = Tkey(ihalf, iplane, imodule, iside, idie, isensor);
2356  //cout<<"setting Offset for: "<<ihalf<<", "<<iplane<<", "<<imodule<<", "<<iside<<", "<<idie<<", "<<isensor<<endl;
2357  itoffset = offsets.find(key);
2358  if (itoffset != offsets.end()) {
2359  x = itoffset->second[0];
2360  y = itoffset->second[1];
2361  z = itoffset->second[2];
2362  rotphi = itoffset->second[3];
2363  rottheta = itoffset->second[4];
2364  rotpsi = itoffset->second[5];
2365  } else {
2366  if (random) {
2367  cout << " generating offset for";
2368  if (ihalf >= 0 && iplane >= 0 && imodule >= 0 && iside >= 0 && idie >= 0
2369  && isensor >= 0) {
2370  cout << " one sensor ";
2371  x = gRandom->Gaus(0, sensor_offset_x);
2372  y = gRandom->Gaus(0, sensor_offset_y);
2373  z = gRandom->Gaus(0, sensor_offset_z);
2374  rotphi = gRandom->Gaus(0, sensor_tilt_phi);
2375  rottheta = gRandom->Gaus(0, sensor_tilt_theta);
2376  rotpsi = gRandom->Gaus(0, sensor_tilt_psi);
2377  // check for clashing sensors starting from the inner
2378  // sensor and correct for it
2379  // check is only performed in the glueing plane
2380  // check along x
2381  // calculate the own additional space required by an offset in x and the rotation
2382  /*
2383  double dx = maps_height * (sqrt(2.)*sin(fabs(rotphi)+3.1416/4.)-1);
2384  if ((isensor > 0 && idie == 0) || (isensor > 1 && idie == 1)){
2385  // check the space required by the preceding sensor
2386  double _x, _y, _z, _rotphi, _rottheta, _rotpsi;
2387  Get_offset(ihalf, iplane, imodule, iside, idie, isensor-1, _x, _y, _z, _rotphi, _rottheta, _rotpsi);
2388  double _dx = maps_height * (sqrt(2.)*sin(fabs(_rotphi)+3.1416/4.)-1);
2389  // check if there is space left
2390  if ( (x+_x) < (_dx+dx) ){
2391  x = (_dx+dx)-(x+_x);
2392  }
2393  }
2394  // same for y neighbors
2395  if (isensor > 0 && idie > 0){
2396  // check the space required by the preceding sensor
2397  double _x, _y, _z, _rotphi, _rottheta, _rotpsi;
2398  Get_offset(ihalf, iplane, imodule, iside, idie-1, isensor, _x, _y, _z, _rotphi, _rottheta, _rotpsi);
2399  double _dy = maps_height * (sqrt(2.)*sin(fabs(_rotphi)+3.1416/4.)-1);
2400  // check if there is space left
2401  if ( (y+_y) < (_dy+dx) ){
2402  y = (_dy+dx)-(y+_y);
2403  }
2404  }*/
2405  //x = (cvd_offset_x);
2406  //y = (cvd_offset_y);
2407  //z = (cvd_offset_z);
2408  //rotphi = (cvd_tilt_phi);
2409  //rottheta = (cvd_tilt_theta);
2410  //rotpsi = (cvd_tilt_psi);
2411  }
2412  // the precision of dies containing sensors is requested
2413  // when requested offset applies to all sensors
2414  // on one side
2415  if (ihalf >= 0 && iplane >= 0 && imodule >= 0 && iside >= 0 && idie < 0
2416  && isensor < 0) {
2417  cout << " one module side ";
2418  x = gRandom->Gaus(0, side_offset_x);
2419  y = gRandom->Gaus(0, side_offset_y);
2420  z = gRandom->Gaus(0, side_offset_z);
2421  rottheta = gRandom->Gaus(0, side_tilt_theta);
2422  rotphi = gRandom->Gaus(0, side_tilt_phi);
2423  rotpsi = gRandom->Gaus(0, side_tilt_psi);
2424  }
2425  // the precision of cvd discs is requested
2426  // when requested offset applies also to the sensors
2427  // on both sides
2428  if (ihalf >= 0 && iplane >= 0 && imodule >= 0 && iside < 0 && idie < 0
2429  && isensor < 0) {
2430  cout << " one module ";
2431  x = gRandom->Gaus(0, cvd_offset_x);
2432  y = gRandom->Gaus(0, cvd_offset_y);
2433  z = gRandom->Gaus(0, cvd_offset_z);
2434  rotphi = gRandom->Gaus(0, cvd_tilt_phi);
2435  rottheta = gRandom->Gaus(0, cvd_tilt_theta);
2436  rotpsi = gRandom->Gaus(0, cvd_tilt_psi);
2437  //x = (cvd_offset_x);
2438  //y = (cvd_offset_y);
2439  //z = (cvd_offset_z);
2440  //rotphi = (cvd_tilt_phi);
2441  //rottheta = (cvd_tilt_theta);
2442  //rotpsi = (cvd_tilt_psi);
2443  }
2444  // the precision of plane halves containing sensors is requested
2445  // when requested offset applies to all modules
2446  if (ihalf >= 0 && iplane >= 0 && imodule < 0 && iside < 0 && idie < 0
2447  && isensor < 0) {
2448  cout << " one plane half ";
2449  x = gRandom->Gaus(0, plane_half_offset_x);
2450  y = gRandom->Gaus(0, plane_half_offset_y);
2451  z = gRandom->Gaus(0, plane_half_offset_z);
2452  rottheta = gRandom->Gaus(0, plane_half_tilt_theta);
2453  rotphi = gRandom->Gaus(0, plane_half_tilt_phi);
2454  rotpsi = gRandom->Gaus(0, plane_half_tilt_psi);
2455  }
2456  // the precision of halves of planes is requested
2457  // when requested offset applies to module supports
2458  if (ihalf >= 0 && iplane < 0 && imodule < 0 && iside < 0 && idie < 0
2459  && isensor < 0) {
2460  cout << " one half ";
2461  x = gRandom->Gaus(0, half_offset_x);
2462  y = gRandom->Gaus(0, half_offset_y);
2463  z = gRandom->Gaus(0, half_offset_z);
2464  rottheta = gRandom->Gaus(0, half_tilt_theta);
2465  cout << " wtf ? " << half_tilt_phi << endl;
2466  rotphi = gRandom->Gaus(0, half_tilt_phi);
2467  rotpsi = gRandom->Gaus(0, half_tilt_psi);
2468  }
2469  // // the precision of the luminosity detector is requested
2470  // // when no degrees of freedom are available to the rest
2471  if (ihalf < 0 && imodule < 0 && iplane < 0 && iside < 0 && idie < 0
2472  && isensor < 0) {
2473  cout << " the luminosity detector ";
2474  // not implemented yet
2475  x = gRandom->Gaus(0, 0);
2476  y = gRandom->Gaus(0, 0);
2477  z = gRandom->Gaus(0, 0);
2478  rottheta = gRandom->Gaus(0, 0);
2479  rotphi = gRandom->Gaus(0, 0);
2480  rotpsi = gRandom->Gaus(0, 0);
2481  }
2482  cout << " of x = " << x << " cm y = " << y << " cm z = " << z;
2483  cout << " cm theta = " << rottheta << " rad phi = " << rotphi;
2484  cout << " rad psi = " << rotpsi << " rad " << endl;
2485  } else {
2486  x = 0;
2487  y = 0;
2488  z = 0;
2489  rotphi = 0;
2490  rottheta = 0;
2491  rotpsi = 0;
2492  }
2493  offsets[key].push_back(x);
2494  offsets[key].push_back(y);
2495  offsets[key].push_back(z);
2496  offsets[key].push_back(rotphi);
2497  offsets[key].push_back(rottheta);
2498  offsets[key].push_back(rotpsi);
2499  }
2500  // cout<<"x,y,z,rotphi,rottheta,rotpsi: "<<x<<", "<<y<<", "<<z<<", "<<rotphi<<", "<<rottheta<<", "<<rotpsi<<endl;
2501 }
2502 
2503 TVector3 PndLmdDim::Decode_hit(const int sensorID, const double column,
2504  const double row, const bool aligned, bool newVersion) {
2505 
2506  if (newVersion) {
2507 
2508  double x, y;
2509  x = column;
2510  y = row;
2511 
2512  //correct for pixel corner to center
2513  x += 0.5;
2514  y += 0.5;
2515 
2516  //shift from corner to center (this considers inactive area!)
2517  x -= (247.5 / 2.0);
2518  y -= (242.5 / 2.0);
2519 
2520  //scale for pixel size
2521  x *= 80e-4;
2522  y *= 80e-4;
2523 
2524  TVector3 hit(x, y, 0.);
2525  int ihalf, iplane, imodule, iside, idie, isensor;
2526  Get_sensor_by_id(sensorID, ihalf, iplane, imodule, iside, idie, isensor);
2527  return Transform_sensor_to_global(hit, ihalf, iplane, imodule, iside, idie,
2528  isensor, false, aligned);
2529  }
2530 
2531  // an additional offset of about 0.5 pixels in x and 2 pixels in y
2532  // seems to appear during digitization via SDS class
2533  // this is taken here into account
2534  double x = (column - 250. / 2. + .5) * maps_active_pixel_size
2535  + maps_passive_left * 2.;
2536  double y = (row - 250. / 2. + 0.5 - 2.5) * maps_active_pixel_size
2537  + maps_passive_bottom * 2.;
2538  TVector3 hit(x, y, 0.);
2539  int ihalf, iplane, imodule, iside, idie, isensor;
2540  Get_sensor_by_id(sensorID, ihalf, iplane, imodule, iside, idie, isensor);
2541  //hit = Transform_sensor_to_lmd_local(hit, ihalf, iplane, imodule, iside, idie, isensor, aligned);
2542  //hit =
2543  return Transform_sensor_to_global(hit, ihalf, iplane, imodule, iside, idie,
2544  isensor, false, aligned);
2545 }
2546 
2547 void PndLmdDim::Propagate_fast_ip_to_lmd(TVector3& pos, TVector3& mom,
2548  double pbeam) {
2549  static bool first_call(true);
2550  if (geometry_version != 2 && !first_call) {
2551  first_call = false;
2552  cout
2553  << " Warning in PndLmdDim::Propagate_fast_ip_to_lmd: Wrong geometry version "
2554  << geometry_version << endl;
2555  cout << " Transformation matrix was fit to data with geometry version 2 ! "
2556  << endl;
2557  }
2558  if (pbeam != 1.5 && pbeam != 15) {
2559  cout << " Error in PndLmdDim::Propagate_fast_ip_to_lmd: no matrix fit for "
2560  << pbeam << " GeV/c beam momentum";
2561  cout << " nor interpolation is implemented yet " << endl;
2562  return;
2563  }
2564  // check if to create and initialize a matrix first
2565  if (propagation_matrices.find(pbeam) == propagation_matrices.end()) {
2566  PndLmdDimPropMat propmat(pbeam);
2567  propagation_matrices[pbeam] = propmat;
2568  }
2569  propagation_matrices[pbeam].Propagate(pos, mom);
2570 }
2571 
2572 bool PndLmdDim::Get_overlapping_sensor(const TVector3& point, int& ihalf,
2573  int& iplane, int& imodule, int& iside, int& idie, int& isensor,
2574  bool aligned) {
2575  bool result = false;
2576  int _ihalf = ihalf;
2577  ihalf = -1;
2578  int _iplane = iplane;
2579  iplane = -1;
2580  int _imodule = imodule;
2581  imodule = -1;
2582  int _iside = iside;
2583  iside = -1;
2584  int _idie = idie;
2585  idie = -1;
2586  int _isensor = isensor;
2587  isensor = -1;
2588  // check if point lies on the proposed sensor
2589  if (!Is_on_Sensor(point, _ihalf, _iplane, _imodule, _iside, _idie, _isensor,
2590  aligned))
2591  return false;
2592  // check all other sensors on the opposite side
2593  if (_iside == 0)
2594  _iside = 1;
2595  else
2596  _iside = 0;
2597  // make it faster by the knowledge, which sensors can only overlap
2598  vector<int> jdie;
2599  vector<int> jsensor;
2600  int nchecks = Get_overlapping_sensor(_idie, _isensor, jdie, jsensor);
2601  //cout << " nchecks " << nchecks << endl;
2602  if (nchecks <= 0)
2603  return false;
2604  for (int icheck = 0; icheck < nchecks; icheck++) {
2605  if (Is_on_Sensor(point, _ihalf, _iplane, _imodule, _iside, jdie[icheck],
2606  jsensor[icheck], aligned)) {
2607  ihalf = _ihalf;
2608  iplane = _iplane;
2609  imodule = _imodule;
2610  iside = _iside;
2611  idie = jdie[icheck];
2612  isensor = jsensor[icheck];
2613  //cout << " true " << endl;
2614  return true;
2615  }
2616  }
2617  return result;
2618 }
2619 
2620 int PndLmdDim::Get_overlapping_sensor(int idie, int isensor, vector<int> &jdie,
2621  vector<int> &jsensor) {
2622  jdie.empty();
2623  jsensor.empty();
2624  if (idie == 0 && isensor == 0) {
2625  jdie.push_back(0);
2626  jsensor.push_back(0);
2627  jdie.push_back(0);
2628  jsensor.push_back(1);
2629  return 2;
2630  }
2631  if (idie == 0 && isensor == 1) {
2632  jdie.push_back(1);
2633  jsensor.push_back(1);
2634  jdie.push_back(0);
2635  jsensor.push_back(1);
2636  jdie.push_back(0);
2637  jsensor.push_back(0);
2638  return 3;
2639  }
2640  if (idie == 0 && isensor == 2) {
2641  jdie.push_back(1);
2642  jsensor.push_back(2);
2643  jdie.push_back(1);
2644  jsensor.push_back(1);
2645  return 2;
2646  }
2647  if (idie == 1 && isensor == 1) {
2648  jdie.push_back(0);
2649  jsensor.push_back(1);
2650  jdie.push_back(0);
2651  jsensor.push_back(2);
2652  jdie.push_back(1);
2653  jsensor.push_back(1);
2654  jdie.push_back(1);
2655  jsensor.push_back(2);
2656  return 4;
2657  }
2658  if (idie == 1 && isensor == 2) {
2659  jdie.push_back(1);
2660  jsensor.push_back(2);
2661  jdie.push_back(0);
2662  jsensor.push_back(2);
2663  jdie.push_back(1);
2664  jsensor.push_back(1);
2665  return 3;
2666  }
2667  return 0;
2668 }
2669 
2670 bool PndLmdDim::Is_on_Sensor(const TVector3& point, int ihalf, int iplane,
2671  int imodule, int iside, int idie, int isensor, bool aligned) {
2672  bool result_x = false;
2673  bool result_y = false;
2674  const TVector3& point_on_sensor = Transform_global_to_sensor(point, ihalf,
2675  iplane, imodule, iside, idie, isensor, false, aligned);
2676  // take only the projection
2677  double x = point_on_sensor.X();
2678  double y = point_on_sensor.Y();
2679  //cout << " checking x " << x << " and y " << y << " on sensor " << ihalf << " " << iplane << " " << imodule << " " << iside << " " << idie << " " << isensor << endl;
2680  if (x <= 0) { // check left border
2681  if ((x + maps_width - maps_passive_left * 2.) >= 0)
2682  result_x = true;
2683  } else { // check right border
2684  if ((x - maps_width + maps_passive_right * 2.) <= 0)
2685  result_x = true;
2686  }
2687  if (y <= 0) { // check bottom border
2688  if ((y + maps_height - maps_passive_bottom * 2.) >= 0)
2689  result_y = true;
2690  } else { // check top border
2691  if ((y - maps_height + maps_passive_top * 2.) <= 0)
2692  result_y = true;
2693  }
2694  //cout << " result is " << result_x << " " << result_y << " therefore " << (result_x && result_y) << endl << endl;
2695  // both must be set to true to be true
2696  return (result_x && result_y);
2697 }
2698 
2699 void PndLmdDim::Transform_global_to_lmd_local(double& x, double& y, double& z,
2700  bool aligned) {
2701  const TGeoHMatrix& matrix = Get_transformation_global_to_lmd_local(aligned);
2702  double from[3] = { x, y, z };
2703  //cout << x << " " << y << " " << z << endl;
2704  double to[3];
2705  matrix.MasterToLocal(from, to);
2706  x = to[0];
2707  y = to[1];
2708  z = to[2];
2709  //cout << x << " " << y << " " << z << endl << endl;
2710 }
2711 
2712 void PndLmdDim::Transform_global_to_lmd_local_vect(double& x, double& y,
2713  double& z, bool aligned) {
2714  const TGeoHMatrix& matrix = Get_transformation_global_to_lmd_local(aligned);
2715  double from[3] = { x, y, z };
2716  //cout << x << " " << y << " " << z << endl;
2717  double to[3];
2718  matrix.MasterToLocalVect(from, to);
2719  x = to[0];
2720  y = to[1];
2721  z = to[2];
2722  //cout << x << " " << y << " " << z << endl << endl;
2723 }
2724 
2725 map<Tkey, TGeoMatrix*>* PndLmdDim::Get_matrices(bool aligned) {
2726  map<Tkey, TGeoMatrix*>* matrices = NULL;
2727  if (aligned) {
2728  matrices = &transformation_matrices_aligned;
2729  } else {
2730  matrices = &transformation_matrices;
2731  }
2732  if (matrices->size() == 0) {
2733  cout
2734  << " Warning in PndLmdDim::Get_matrices: No transformation matrices loaded! => trying to load default ones. "
2735  << endl;
2736  Read_transformation_matrices("", aligned);
2737  }
2738  return matrices;
2739 }
2740 
2741 TGeoMatrix* PndLmdDim::Get_matrix(int ihalf, int iplane, int imodule, int iside,
2742  int idie, int isensor, bool aligned) {
2743  map<Tkey, TGeoMatrix*>* matrices = Get_matrices(aligned);
2744  if (!matrices)
2745  return NULL;
2746  it_transformation_matrices = matrices->find(
2747  Tkey(ihalf, iplane, imodule, iside, idie, isensor));
2748  if (it_transformation_matrices == matrices->end()) {
2749  //cout << " Warning in PndLmdDim::Get_matrix: Transformation matrix not existent! " << endl;
2750  //cout << ihalf << " " << iplane << " " << imodule << " " << iside << " " << idie << " " << isensor << endl;
2751  return NULL;
2752  } else {
2753  return it_transformation_matrices->second;
2754  }
2755 }
2756 
2757 bool PndLmdDim::Get_matrix_difference(int ihalf, int iplane, int imodule,
2758  int iside, int idie, int isensor, double& dx, double& dy, double& dz,
2759  double& dphi, double& dtheta, double& dpsi) {
2760  bool result = true;
2761  dx = 0;
2762  dy = 0;
2763  dz = 0;
2764  dphi = 0;
2765  dtheta = 0;
2766  dpsi = 0;
2767  // try to retrieve the matrix
2768  TGeoMatrix* matrix = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
2769  false);
2770  TGeoMatrix* matrix_aligned = Get_matrix(ihalf, iplane, imodule, iside, idie,
2771  isensor, true);
2772  if (!matrix || !matrix_aligned)
2773  return false;
2774  TGeoCombiTrans combi_trans_matrix(*matrix);
2775  TGeoCombiTrans combi_trans_matrix_aligned(*matrix_aligned);
2776  if (combi_trans_matrix.IsTranslation()
2777  && combi_trans_matrix_aligned.IsTranslation()) {
2778  const double* pos = combi_trans_matrix.GetTranslation();
2779  const double* pos_aligned = combi_trans_matrix_aligned.GetTranslation();
2780  dx = pos[0] - pos_aligned[0];
2781  dy = pos[1] - pos_aligned[1];
2782  dz = pos[2] - pos_aligned[2];
2783  }
2784  if (combi_trans_matrix.IsRotation()
2785  && combi_trans_matrix_aligned.IsRotation()) {
2786  double rot[3], rot_aligned[3];
2787  combi_trans_matrix.GetRotation()->GetAngles(rot[0], rot[1], rot[2]);
2788  combi_trans_matrix_aligned.GetRotation()->GetAngles(rot_aligned[0],
2789  rot_aligned[1], rot_aligned[2]);
2790  dphi = rot[0] - rot_aligned[0];
2791  dtheta = rot[1] - rot_aligned[1];
2792  dpsi = rot[2] - rot_aligned[2];
2793  }
2794  /*if (ihalf==-1){
2795  dx = pos[0];
2796  dy = pos[1];
2797  dz = pos[2];
2798  dphi = rot[0];
2799  dtheta = rot[1];
2800  dpsi = rot[2];
2801  }*/
2802  // format the values to something meaning full for the lumi
2803  dx *= 10.e3; // cm -> µm
2804  dy *= 10.e3; // cm -> µm
2805  dz *= 10.e3; // cm -> µm
2806 
2807  dphi = dphi / 180. * pi * 1.e6; // deg -> µrad
2808  dtheta = dtheta / 180. * pi * 1.e6; // deg -> µrad
2809  dpsi = dpsi / 180. * pi * 1.e6; // deg -> µrad
2810  return result;
2811 }
2812 
2813 #include<TCanvas.h>
2814 void PndLmdDim::Calc_matrix_offsets() {
2815  // go through all available matrices and calculate the differences between in terms of
2816  // displacement and rotation
2817  int ihalf(-1), iplane(-1), imodule(-1), iside(-1), idie(-1), isensor(-1);
2818  double dx(0), dy(0), dz(0), dphi(0), dtheta(0), dpsi(0);
2819  TH1D hist_dx_sensors_local("hist_dx_sensors_local",
2820  "ref local;#Deltax [#mum]; entries", 100, -1e3, 1e3);
2821  TH1D hist_dx_sensors_global("hist_dx_sensors_lmd ",
2822  "ref lmd;#Deltax [#mum]; entries", 100, -1e3, 1e3);
2823  TH1D hist_dy_sensors_local("hist_dy_sensors_local",
2824  "ref local;#Deltay [#mum]; entries", 100, -1e3, 1e3);
2825  TH1D hist_dy_sensors_global("hist_dy_sensors_lmd",
2826  "ref lmd;#Deltay [#mum]; entries", 100, -1e3, 1e3);
2827  TH1D hist_dz_sensors_local("hist_dz_sensors_local",
2828  "ref local;#Deltaz [#mum]; entries", 100, -1e3, 1e3);
2829  TH1D hist_dz_sensors_global("hist_dz_sensors_lmd",
2830  "ref lmd;#Deltaz [#mum]; entries", 100, -1e3, 1e3);
2831  TH1D hist_dphi_sensors_local("hist_dphi_sensors_local",
2832  "ref local;#Delta#phi [#murad]; entries", 100, -1e3, 1e3);
2833  TH1D hist_dphi_sensors_global("hist_dphi_sensors_lmd",
2834  "ref lmd;#Delta#phi [#murad]; entries", 100, -1e3, 1e3);
2835  TH1D hist_dtheta_sensors_local("hist_dtheta_sensors_local",
2836  "ref local;#Delta#theta [#murad]; entries", 100, -1e3, 1e3);
2837  TH1D hist_dtheta_sensors_global("hist_dtheta_sensors_lmd",
2838  "ref lmd;#Delta#theta [#murad]; entries", 100, -1e3, 1e3);
2839  TH1D hist_dpsi_sensors_local("hist_dpsi_sensors_local",
2840  "ref local;#Delta#psi [#murad]; entries", 100, -1e3, 1e3);
2841  TH1D hist_dpsi_sensors_global("hist_dpsi_sensors_lmd",
2842  "ref lmd;#Delta#psi [#murad]; entries", 100, -1e3, 1e3);
2843 
2844  ofstream outfile("offsets.txt");
2845  if (outfile.is_open()) {
2846  outfile.setf(std::ios::fixed, std::ios::floatfield); // floatfield not set
2847  outfile.precision(7);
2848 
2849  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie, isensor, dx,
2850  dy, dz, dphi, dtheta, dpsi)) {
2851  outfile << " coordinate offset of the lmd reference system \n";
2852  outfile
2853  << " \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2854  outfile << dx << "\t" << dy << "\t" << dz << "\t" << dphi << "\t"
2855  << dtheta << "\t" << dpsi << "\n";
2856  outfile << endl;
2857  }
2858  for (ihalf = 0; ihalf < 2; ihalf++) {
2859  iplane = -1;
2860  imodule = -1;
2861  iside = -1;
2862  idie = -1;
2863  isensor = -1;
2864  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie, isensor,
2865  dx, dy, dz, dphi, dtheta, dpsi)) {
2866  outfile << "\t coordinate offset of the lmd half " << ihalf << " \n";
2867  outfile
2868  << "\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2869  outfile << "\t" << dx << "\t" << dy << "\t" << dz << "\t" << dphi
2870  << "\t" << dtheta << "\t" << dpsi << "\n";
2871  outfile << endl;
2872  }
2873  for (iplane = 0; iplane < (int)n_planes; iplane++) {
2874  imodule = -1;
2875  iside = -1;
2876  idie = -1;
2877  isensor = -1;
2878  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie, isensor,
2879  dx, dy, dz, dphi, dtheta, dpsi)) {
2880  outfile << "\t\t coordinate offset of the lmd plane " << iplane
2881  << " \n";
2882  outfile
2883  << "\t\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2884  outfile << "\t\t" << dx << "\t" << dy << "\t" << dz << "\t" << dphi
2885  << "\t" << dtheta << "\t" << dpsi << "\n";
2886  outfile << endl;
2887  }
2888  for (imodule = 0; imodule < n_cvd_discs / 2; imodule++) {
2889  iside = -1;
2890  idie = -1;
2891  isensor = -1;
2892  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie,
2893  isensor, dx, dy, dz, dphi, dtheta, dpsi)) {
2894  outfile << "\t\t\t coordinate offset of the lmd module " << imodule
2895  << " \n";
2896  outfile
2897  << "\t\t\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2898  outfile << "\t\t\t" << dx << "\t" << dy << "\t" << dz << "\t"
2899  << dphi << "\t" << dtheta << "\t" << dpsi << "\n";
2900  outfile << endl;
2901  }
2902  for (iside = 0; iside < 2; iside++) {
2903  idie = -1;
2904  isensor = -1;
2905  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie,
2906  isensor, dx, dy, dz, dphi, dtheta, dpsi)) {
2907  outfile << "\t\t\t\t coordinate offset of the lmd module side "
2908  << iside << " \n";
2909  outfile
2910  << "\t\t\t\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2911  outfile << "\t\t\t\t" << dx << "\t" << dy << "\t" << dz << "\t"
2912  << dphi << "\t" << dtheta << "\t" << dpsi << "\n";
2913  outfile << endl;
2914  }
2915  for (idie = 0; idie < 2; idie++) {
2916  isensor = -1;
2917  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie,
2918  isensor, dx, dy, dz, dphi, dtheta, dpsi)) {
2919  outfile << "\t\t\t\t\t coordinate offset of the lmd die "
2920  << idie << " \n";
2921  outfile
2922  << "\t\t\t\t\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2923  outfile << "\t\t\t\t\t" << dx << "\t" << dy << "\t" << dz
2924  << "\t" << dphi << "\t" << dtheta << "\t" << dpsi << "\n";
2925  outfile << endl;
2926  }
2927  for (isensor = 0; isensor < 3; isensor++) {
2928  if (Get_matrix_difference(ihalf, iplane, imodule, iside, idie,
2929  isensor, dx, dy, dz, dphi, dtheta, dpsi)) {
2930  outfile << "\t\t\t\t\t\t coordinate offset of the lmd sensor "
2931  << isensor << " \n";
2932  outfile
2933  << "\t\t\t\t\t\t \u0394x [µm] \t \u0394y [µm] \t \u0394z [µm] \t \u0394\u03C6 [µrad] \t \u0394\u03D1 [µrad] \t \u0394\u03C8 [µrad] \t \n";
2934  outfile << "\t\t\t\t\t\t" << dx << "\t" << dy << "\t" << dz
2935  << "\t" << dphi << "\t" << dtheta << "\t" << dpsi << "\n";
2936  outfile << endl;
2937  hist_dx_sensors_local.Fill(dx);
2938  hist_dy_sensors_local.Fill(dy);
2939  hist_dz_sensors_local.Fill(dz);
2940  hist_dphi_sensors_local.Fill(dphi);
2941  hist_dtheta_sensors_local.Fill(dtheta);
2942  hist_dpsi_sensors_local.Fill(dpsi);
2943  }
2944  } // loop over the sensor of one die
2945  } // loop over the die of one side module
2946  } // loop over the sides of one module
2947  } // loop over the modules of one plane half
2948  } // loop over the planes of the detector
2949  } // loop over the halves of the detector
2950  outfile.close();
2951  TCanvas outcanvas("canvas", "canvas", 600, 600);
2952  hist_dx_sensors_local.Draw();
2953  outcanvas.Print("offsets.pdf(");
2954  hist_dy_sensors_local.Draw();
2955  outcanvas.Print("offsets.pdf(");
2956  hist_dz_sensors_local.Draw();
2957  outcanvas.Print("offsets.pdf(");
2958  hist_dphi_sensors_local.Draw();
2959  outcanvas.Print("offsets.pdf(");
2960  hist_dtheta_sensors_local.Draw();
2961  outcanvas.Print("offsets.pdf(");
2962  hist_dpsi_sensors_local.Draw();
2963  outcanvas.Print("offsets.pdf)");
2964  } else {
2965  cout << " sorry, could not write file into the current directory. " << endl;
2966  // if file is open
2967  }
2968 }
2969 
2970 TGeoMatrix* PndLmdDim::Get_matrix_global_to_lmd_local(bool aligned) {
2971  return Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
2972 }
2973 
2974 // TGeoMatrix* PndLmdDim::Get_matrix_global_to_sector_local(int ihalf, int iplane, int imodule, bool aligned){
2975 // return Get_matrix(ihalf, iplane, imodule,-1,-1,-1, aligned);
2976 // }
2977 
2978 TGeoMatrix* PndLmdDim::Get_matrix_lmd_local_to_module_side(int ihalf,
2979  int iplane, int imodule, int iside, bool aligned) {
2980  return Get_matrix(ihalf, iplane, imodule, iside, -1, -1, aligned);
2981 }
2982 
2983 TGeoMatrix* PndLmdDim::Get_matrix_module_side_to_sensor(int ihalf, int iplane,
2984  int imodule, int iside, int idie, int isensor, bool aligned) {
2985  return Get_matrix(ihalf, iplane, imodule, iside, idie, isensor, aligned);
2986 }
2987 
2988 TGeoHMatrix PndLmdDim::Get_transformation_global_to_lmd_local(bool aligned) {
2989  TGeoMatrix *matrix = Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
2990  if (!matrix)
2991  return TGeoHMatrix();
2992  return TGeoHMatrix(*matrix);
2993 }
2994 
2995 TGeoHMatrix PndLmdDim::Get_transformation_lmd_local_to_module_side(int ihalf,
2996  int iplane, int imodule, int iside, bool aligned) {
2997  TGeoMatrix* matrix = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
2998  aligned);
2999  if (!matrix)
3000  return TGeoHMatrix();
3001  return TGeoHMatrix(*matrix);
3002 }
3003 
3004 TGeoHMatrix PndLmdDim::Get_transformation_module_side_to_sensor(int ihalf,
3005  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3006  TGeoMatrix* matrix = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3007  aligned);
3008  if (!matrix)
3009  TGeoHMatrix();
3010  return TGeoHMatrix(*matrix);
3011 }
3012 
3013 TGeoHMatrix PndLmdDim::Get_transformation_global_to_sensor(int ihalf,
3014  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3015  // instead of using the 3 methods to transform between the reference frames
3016  // pointers to the matrices are retrieved directly in order to gain a little
3017  // bit of performance (no construction and deletion of new matrices needed)
3018  TGeoMatrix* matrix1 = Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
3019  TGeoMatrix* matrix2 = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
3020  aligned);
3021  TGeoMatrix* matrix3 = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3022  aligned);
3023  if (!matrix1 || !matrix2 || !matrix3)
3024  return TGeoHMatrix();
3025  //return TGeoHMatrix(((*matrix1) * (*matrix2)) * (*matrix3)); //FIXME: I removed that line
3026  return TGeoHMatrix((*matrix1) * (*matrix2) * (*matrix3));
3027 }
3028 
3029 TGeoHMatrix PndLmdDim::Get_transformation_lmd_local_to_sensor(int ihalf,
3030  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3031  TGeoMatrix* matrix2 = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
3032  aligned);
3033  TGeoMatrix* matrix3 = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3034  aligned);
3035  if (!matrix2 || !matrix3)
3036  return TGeoHMatrix();
3037  return TGeoHMatrix((*matrix2) * (*matrix3));
3038 }
3039 
3040 TGeoHMatrix PndLmdDim::Get_transformation_lmd_local_to_global(bool aligned) {
3041  TGeoMatrix* matrix = Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
3042  if (!matrix)
3043  return TGeoHMatrix();
3044  return TGeoHMatrix(matrix->Inverse());
3045 }
3046 
3047 TGeoHMatrix PndLmdDim::Get_transformation_module_side_to_lmd_local(int ihalf,
3048  int iplane, int imodule, int iside, bool aligned) {
3049  TGeoMatrix* matrix = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
3050  aligned);
3051  if (!matrix)
3052  TGeoHMatrix();
3053  return TGeoHMatrix(matrix->Inverse());
3054 }
3055 
3056 TGeoHMatrix PndLmdDim::Get_transformation_sensor_to_module_side(int ihalf,
3057  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3058  TGeoMatrix* matrix = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3059  aligned);
3060  if (!matrix)
3061  return TGeoHMatrix();
3062  return TGeoHMatrix(matrix->Inverse());
3063 }
3064 
3065 TGeoHMatrix PndLmdDim::Get_transformation_sensor_to_global(int ihalf,
3066  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3067  // instead of using the 3 methods to transform between the reference frames
3068  // pointers to the matrices are retrieved directly in order to gain a little
3069  // bit of performance (no construction and deletion of new matrices needed)
3070  TGeoMatrix* matrix1 = Get_matrix(-1, -1, -1, -1, -1, -1, aligned);
3071  TGeoMatrix* matrix2 = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
3072  aligned);
3073  TGeoMatrix* matrix3 = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3074  aligned);
3075  if (!matrix1 || !matrix2 || !matrix3)
3076  return TGeoHMatrix();
3077  return TGeoHMatrix(((*matrix1) * (*matrix2) * (*matrix3)).Inverse());
3078 }
3079 
3080 //FIXME: above and below a return is missing and was added;
3081 
3082 TGeoHMatrix PndLmdDim::Get_transformation_sensor_to_lmd_local(int ihalf,
3083  int iplane, int imodule, int iside, int idie, int isensor, bool aligned) {
3084  TGeoMatrix* matrix2 = Get_matrix(ihalf, iplane, imodule, iside, -1, -1,
3085  aligned);
3086  TGeoMatrix* matrix3 = Get_matrix(ihalf, iplane, imodule, iside, idie, isensor,
3087  aligned);
3088  if (!matrix2 || !matrix3)
3089  return TGeoHMatrix();
3090  return TGeoHMatrix(((*matrix2) * (*matrix3)).Inverse());
3091 }
3092 
3093 TGeoHMatrix PndLmdDim::Get_transformation_sensor_to_sensor_aligned(int ihalf,
3094  int iplane, int imodule, int iside, int idie, int isensor) {
3095  // to do: optimize like in Get_transformation_sensor_to_global
3096  const TGeoHMatrix& matrix = Get_transformation_sensor_to_global(ihalf, iplane,
3097  imodule, iside, idie, isensor, false);
3098  const TGeoHMatrix& matrix_aligned = Get_transformation_global_to_sensor(ihalf,
3099  iplane, imodule, iside, idie, isensor, true);
3100  return TGeoHMatrix(matrix * matrix_aligned);
3101 }
3102 
3103 TGeoHMatrix PndLmdDim::Get_transformation_sensor_aligned_to_sensor(int ihalf,
3104  int iplane, int imodule, int iside, int idie, int isensor) {
3105  // to do: optimize like in Get_transformation_sensor_to_global
3106  const TGeoHMatrix& matrix_aligned = Get_transformation_sensor_to_global(ihalf,
3107  iplane, imodule, iside, idie, isensor, true);
3108  const TGeoHMatrix& matrix = Get_transformation_global_to_sensor(ihalf, iplane,
3109  imodule, iside, idie, isensor, false);
3110  return TGeoHMatrix(matrix_aligned * matrix);
3111 }
3112 
3113 TVector3 PndLmdDim::Transform_global_to_lmd_local(const TVector3& point,
3114  bool isvector, bool aligned) {
3115  const TGeoHMatrix& matrix = Get_transformation_global_to_lmd_local(aligned);
3116  double master[3];
3117  point.GetXYZ(master);
3118  double local[3];
3119  if (isvector)
3120  matrix.MasterToLocalVect(master, local);
3121  else
3122  matrix.MasterToLocal(master, local);
3123  return TVector3(local);
3124 }
3125 
3126 TVector3 PndLmdDim::Transform_lmd_local_to_module_side(const TVector3& point,
3127  int ihalf, int iplane, int imodule, int iside, bool isvector,
3128  bool aligned) {
3129  const TGeoHMatrix& matrix = Get_transformation_lmd_local_to_module_side(ihalf,
3130  iplane, imodule, iside, aligned);
3131  double master[3];
3132  point.GetXYZ(master);
3133  double local[3];
3134  if (isvector)
3135  matrix.MasterToLocalVect(master, local);
3136  else
3137  matrix.MasterToLocal(master, local);
3138  return TVector3(local);
3139 }
3140 
3141 TVector3 PndLmdDim::Transform_module_side_to_sensor(const TVector3& point,
3142  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3143  bool isvector, bool aligned) {
3144  const TGeoHMatrix& matrix = Get_transformation_module_side_to_sensor(ihalf,
3145  iplane, imodule, iside, idie, isensor, aligned);
3146  double master[3];
3147  point.GetXYZ(master);
3148  double local[3];
3149  if (isvector)
3150  matrix.MasterToLocalVect(master, local);
3151  else
3152  matrix.MasterToLocal(master, local);
3153  return TVector3(local);
3154 }
3155 
3156 TVector3 PndLmdDim::Transform_global_to_sensor(const TVector3& point, int ihalf,
3157  int iplane, int imodule, int iside, int idie, int isensor, bool isvector,
3158  bool aligned) {
3159  const TGeoHMatrix& matrix = Get_transformation_global_to_sensor(ihalf, iplane,
3160  imodule, iside, idie, isensor, aligned);
3161  double master[3];
3162  point.GetXYZ(master);
3163  double local[3];
3164  if (isvector)
3165  matrix.MasterToLocalVect(master, local);
3166  else
3167  matrix.MasterToLocal(master, local);
3168  return TVector3(local);
3169 }
3170 
3171 TVector3 PndLmdDim::Transform_lmd_local_to_global(const TVector3& point,
3172  bool isvector, bool aligned) {
3173  // I think Local to Master calculation is faster than the getter of the inverse matrix
3174  const TGeoHMatrix& matrix = Get_transformation_global_to_lmd_local(aligned);
3175  double local[3];
3176  point.GetXYZ(local);
3177  double master[3];
3178  if (isvector)
3179  matrix.LocalToMasterVect(local, master);
3180  else
3181  matrix.LocalToMaster(local, master);
3182  return TVector3(master);
3183 }
3184 
3185 TVector3 PndLmdDim::Transform_lmd_local_to_sensor(const TVector3& point,
3186  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3187  bool isvector, bool aligned) {
3188  const TGeoHMatrix& matrix = Get_transformation_lmd_local_to_sensor(ihalf,
3189  iplane, imodule, iside, idie, isensor, aligned);
3190  double master[3];
3191  point.GetXYZ(master);
3192  double local[3];
3193  if (isvector)
3194  matrix.MasterToLocalVect(master, local);
3195  else
3196  matrix.MasterToLocal(master, local);
3197  return TVector3(local);
3198 }
3199 
3200 TVector3 PndLmdDim::Transform_module_side_to_lmd_local(const TVector3& point,
3201  int ihalf, int iplane, int imodule, int iside, bool isvector,
3202  bool aligned) {
3203  // I think Local to Master calculation is faster than the getter of the inverse matrix
3204  const TGeoHMatrix& matrix = Get_transformation_lmd_local_to_module_side(ihalf,
3205  iplane, imodule, iside, aligned);
3206  double local[3];
3207  point.GetXYZ(local);
3208  double master[3];
3209  if (isvector)
3210  matrix.LocalToMasterVect(local, master);
3211  else
3212  matrix.LocalToMaster(local, master);
3213  return TVector3(master);
3214 }
3215 
3216 TVector3 PndLmdDim::Transform_sensor_to_module_side(const TVector3& point,
3217  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3218  bool isvector, bool aligned) {
3219  // I think Local to Master calculation is faster than the getter of the inverse matrix
3220  const TGeoHMatrix& matrix = Get_transformation_module_side_to_sensor(ihalf,
3221  iplane, imodule, iside, idie, isensor, aligned);
3222  double local[3];
3223  point.GetXYZ(local);
3224  double master[3];
3225  if (isvector)
3226  matrix.LocalToMasterVect(local, master);
3227  else
3228  matrix.LocalToMaster(local, master);
3229  return TVector3(master);
3230 }
3231 
3232 TVector3 PndLmdDim::Transform_sensor_to_global(const TVector3& point, int ihalf,
3233  int iplane, int imodule, int iside, int idie, int isensor, bool isvector,
3234  bool aligned) {
3235  // I think Local to Master calculation is faster than the getter of the inverse matrix
3236  const TGeoHMatrix& matrix = Get_transformation_global_to_sensor(ihalf, iplane,
3237  imodule, iside, idie, isensor, aligned);
3238  double local[3];
3239  point.GetXYZ(local);
3240  double master[3];
3241  if (isvector)
3242  matrix.LocalToMasterVect(local, master);
3243  else
3244  matrix.LocalToMaster(local, master);
3245  return TVector3(master);
3246 }
3247 
3248 TVector3 PndLmdDim::Transform_sensor_to_lmd_local(const TVector3& point,
3249  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3250  bool isvector, bool aligned) {
3251  const TGeoHMatrix& matrix = Get_transformation_lmd_local_to_sensor(ihalf,
3252  iplane, imodule, iside, idie, isensor, aligned);
3253  double local[3];
3254  point.GetXYZ(local);
3255  double master[3];
3256  if (isvector)
3257  matrix.LocalToMasterVect(local, master);
3258  else
3259  matrix.LocalToMaster(local, master);
3260  return TVector3(master);
3261 }
3262 
3263 TVector3 PndLmdDim::Transform_sensor_to_sensor_aligned(const TVector3& point,
3264  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3265  bool isvector) {
3266  const TGeoHMatrix& matrix = Get_transformation_sensor_to_sensor_aligned(ihalf,
3267  iplane, imodule, iside, idie, isensor);
3268  double master[3];
3269  point.GetXYZ(master);
3270  double local[3];
3271  if (isvector)
3272  matrix.MasterToLocalVect(master, local);
3273  else
3274  matrix.MasterToLocal(master, local);
3275  return TVector3(local);
3276 }
3277 
3278 TVector3 PndLmdDim::Transform_sensor_aligned_to_sensor(const TVector3& point,
3279  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3280  bool isvector) {
3281  const TGeoHMatrix& matrix = Get_transformation_sensor_aligned_to_sensor(ihalf,
3282  iplane, imodule, iside, idie, isensor);
3283  double master[3];
3284  point.GetXYZ(master);
3285  double local[3];
3286  if (isvector)
3287  matrix.MasterToLocalVect(master, local);
3288  else
3289  matrix.MasterToLocal(master, local);
3290  return TVector3(local);
3291 }
3292 
3293 TVector3 PndLmdDim::Transform_sensor_to_sensor(const TVector3& point,
3294  int ihalf_from, int iplane_from, int imodule_from, int iside_from,
3295  int idie_from, int isensor_from, int ihalf_to, int iplane_to,
3296  int imodule_to, int iside_to, int idie_to, int isensor_to, bool isvector,
3297  bool aligned) {
3298  const TGeoHMatrix& matrix_from = Get_transformation_sensor_to_lmd_local(
3299  ihalf_from, iplane_from, imodule_from, iside_from, idie_from,
3300  isensor_from, aligned);
3301  const TGeoHMatrix& matrix_to = Get_transformation_lmd_local_to_sensor(
3302  ihalf_to, iplane_to, imodule_to, iside_to, idie_to, isensor_to, aligned);
3303  TGeoHMatrix matrix(matrix_from * matrix_to);
3304  double master[3];
3305  point.GetXYZ(master);
3306  double local[3];
3307  if (isvector)
3308  matrix.MasterToLocalVect(master, local);
3309  else
3310  matrix.MasterToLocal(master, local);
3311  return TVector3(local);
3312 }
3313 
3314 TMatrixD PndLmdDim::Transform_global_to_lmd_local(const TMatrixD& matrix,
3315  bool aligned) {
3316  TMatrixD rotmatrix(3, 3,
3317  Get_transformation_global_to_lmd_local(aligned).GetRotationMatrix());
3318  return TMatrixD(
3319  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3320 }
3321 
3322 TMatrixD PndLmdDim::Transform_lmd_local_to_module_side(const TMatrixD& matrix,
3323  int ihalf, int iplane, int imodule, int iside, bool aligned) {
3324  TMatrixD rotmatrix(3, 3,
3325  Get_transformation_lmd_local_to_module_side(ihalf, iplane, imodule, iside,
3326  aligned).GetRotationMatrix());
3327  return TMatrixD(
3328  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3329 }
3330 
3331 TMatrixD PndLmdDim::Transform_module_side_to_sensor(const TMatrixD& matrix,
3332  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3333  bool aligned) {
3334  TMatrixD rotmatrix(3, 3,
3335  Get_transformation_module_side_to_sensor(ihalf, iplane, imodule, iside,
3336  idie, isensor, aligned).GetRotationMatrix());
3337  return TMatrixD(
3338  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3339 }
3340 
3341 TMatrixD PndLmdDim::Transform_global_to_sensor(const TMatrixD& matrix,
3342  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3343  bool aligned) {
3344  TMatrixD rotmatrix(3, 3,
3345  Get_transformation_global_to_sensor(ihalf, iplane, imodule, iside, idie,
3346  isensor, aligned).GetRotationMatrix());
3347  return TMatrixD(
3348  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3349 }
3350 
3351 TMatrixD PndLmdDim::Transform_lmd_local_to_sensor(const TMatrixD& matrix,
3352  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3353  bool aligned) {
3354  TMatrixD rotmatrix(3, 3,
3355  Get_transformation_lmd_local_to_sensor(ihalf, iplane, imodule, iside,
3356  idie, isensor, aligned).GetRotationMatrix());
3357  return TMatrixD(
3358  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3359 }
3360 
3361 TMatrixD PndLmdDim::Transform_lmd_local_to_global(const TMatrixD& matrix,
3362  bool aligned) {
3363  TMatrixD rotmatrix(3, 3,
3364  Get_transformation_lmd_local_to_global(aligned).GetRotationMatrix());
3365  return TMatrixD(
3366  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3367 }
3368 
3369 TMatrixD PndLmdDim::Transform_module_side_to_lmd_local(const TMatrixD& matrix,
3370  int ihalf, int iplane, int imodule, int iside, bool aligned) {
3371  TMatrixD rotmatrix(3, 3,
3372  Get_transformation_module_side_to_lmd_local(ihalf, iplane, imodule, iside,
3373  aligned).GetRotationMatrix());
3374  return TMatrixD(
3375  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3376 }
3377 
3378 TMatrixD PndLmdDim::Transform_sensor_to_module_side(const TMatrixD& matrix,
3379  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3380  bool aligned) {
3381  TMatrixD rotmatrix(3, 3,
3382  Get_transformation_sensor_to_module_side(ihalf, iplane, imodule, iside,
3383  idie, isensor, aligned).GetRotationMatrix());
3384  return TMatrixD(
3385  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3386 }
3387 
3388 TMatrixD PndLmdDim::Transform_sensor_to_global(const TMatrixD& matrix,
3389  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3390  bool aligned) {
3391  TMatrixD rotmatrix(3, 3,
3392  Get_transformation_sensor_to_global(ihalf, iplane, imodule, iside, idie,
3393  isensor, aligned).GetRotationMatrix());
3394  return TMatrixD(
3395  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3396 }
3397 
3398 TMatrixD PndLmdDim::Transform_sensor_to_lmd_local(const TMatrixD& matrix,
3399  int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
3400  bool aligned) {
3401  TMatrixD rotmatrix(3, 3,
3402  Get_transformation_sensor_to_lmd_local(ihalf, iplane, imodule, iside,
3403  idie, isensor, aligned).GetRotationMatrix());
3404  return TMatrixD(
3405  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3406 }
3407 
3408 TMatrixD PndLmdDim::Transform_sensor_to_sensor_aligned(const TMatrixD& matrix,
3409  int ihalf, int iplane, int imodule, int iside, int idie, int isensor) {
3410  TMatrixD rotmatrix(3, 3,
3411  Get_transformation_sensor_to_sensor_aligned(ihalf, iplane, imodule, iside,
3412  idie, isensor).GetRotationMatrix());
3413  return TMatrixD(
3414  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3415 }
3416 
3417 TMatrixD PndLmdDim::Transform_sensor_aligned_to_sensor(const TMatrixD& matrix,
3418  int ihalf, int iplane, int imodule, int iside, int idie, int isensor) {
3419  TMatrixD rotmatrix(3, 3,
3420  Get_transformation_sensor_aligned_to_sensor(ihalf, iplane, imodule, iside,
3421  idie, isensor).GetRotationMatrix());
3422  return TMatrixD(
3423  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3424 }
3425 
3426 TMatrixD PndLmdDim::Transform_sensor_to_sensor(const TMatrixD& matrix,
3427  int ihalf_from, int iplane_from, int imodule_from, int iside_from,
3428  int idie_from, int isensor_from, int ihalf_to, int iplane_to,
3429  int imodule_to, int iside_to, int idie_to, int isensor_to, bool aligned) {
3430  TMatrixD rotmatrix(3, 3,
3431  (Get_transformation_sensor_to_lmd_local(ihalf_from, iplane_from,
3432  imodule_from, iside_from, idie_from, isensor_from, aligned)
3433  * Get_transformation_lmd_local_to_sensor(ihalf_to, iplane_to,
3434  imodule_to, iside_to, idie_to, isensor_to, aligned)).GetRotationMatrix());
3435  return TMatrixD(
3436  rotmatrix * TMatrixD(matrix, TMatrixD::kMultTranspose, rotmatrix));
3437 }
3438 
3439 void PndLmdDim::Test_matrices() {
3440  //TVector3 from(1.,2.,3.);
3441  //TVector3 to = Transform_global_to_lmd_local(from);
3442  //to = Transform_lmd_local_to_global(to);
3443  //if ((from-to).Mag() > 1e7) cout << " error " << endl;
3444  for (unsigned int ihalf = 0; ihalf < 2; ihalf++)
3445  for (unsigned int iplane = 0; iplane < n_planes; iplane++)
3446  for (unsigned int imodule = 0; imodule < nmodules; imodule++)
3447  for (unsigned int iside = 0; iside < 2; iside++)
3448  for (unsigned int idie = 0; idie < 2; idie++)
3449  for (unsigned int isensor = 0; isensor < 3; isensor++) {
3450  if (idie == 1 && isensor == 0)
3451  continue;
3452  TVector3 from_sens(1., 2., 3.);
3453  TVector3 to_sens = Transform_global_to_sensor(from_sens, ihalf,
3454  iplane, imodule, iside, idie, isensor);
3455  to_sens = Transform_sensor_to_global(to_sens, ihalf, iplane,
3456  imodule, iside, idie, isensor);
3457  if ((from_sens - to_sens).Mag() > 1e7)
3458  cout << " error " << endl;
3459  }
3460 }
3461 
3462 void PndLmdDim::Draw_Sensors(int iplane, bool aligned, bool lmd_frame,
3463  int glside) {
3464 
3465  for (unsigned int ihalf = 0; ihalf < 2; ihalf++)
3466  //for (unsigned int iplane = 0; iplane < n_planes; iplane++)
3467  for (unsigned int imodule = 0; imodule < nmodules; imodule++) {
3468  if (glside > 1) { //default
3469  for (unsigned int iside = 0; iside < 2; iside++) {
3470  for (unsigned int idie = 0; idie < 2; idie++) {
3471  for (unsigned int isensor = 0; isensor < 3; isensor++)
3472  //for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
3473  {
3474  //cout << " sensID " << sensorID << endl;
3475  //int ihalf, _iplane, imodule, iside, idie, isensor;
3476  //Get_sensor_by_id(sensorID, ihalf, _iplane, imodule, iside, idie, isensor);
3477  //cout /*<< sensorID*/ << " " << ihalf << " " << iplane << " " << imodule << " " << iside << " " << idie << " " << isensor << endl;
3478  //if (iplane != _iplane) continue;
3479  //cout << " " << ihalf << " " << iplane << endl;
3480  if (idie == 1 && isensor == 0)
3481  continue;
3482  TPolyLine* sensor_shape = Get_Sensor_Shape(ihalf, iplane, imodule,
3483  iside, idie, isensor, aligned, lmd_frame);
3484  if (iside == 1)
3485  sensor_shape->SetLineColor(17);
3486  else
3487  sensor_shape->SetLineColor(13);
3488  sensor_shape->Draw();
3489  }
3490  }
3491  }
3492  } //default
3493  else {
3494  int iside = glside;
3495  for (unsigned int idie = 0; idie < 2; idie++) {
3496  for (unsigned int isensor = 0; isensor < 3; isensor++) {
3497  if (idie == 1 && isensor == 0)
3498  continue;
3499  TPolyLine* sensor_shape = Get_Sensor_Shape(ihalf, iplane, imodule,
3500  iside, idie, isensor, aligned, lmd_frame);
3501  if (iside == 1)
3502  sensor_shape->SetLineColor(17);
3503  else
3504  sensor_shape->SetLineColor(13);
3505  sensor_shape->Draw();
3506  }
3507  }
3508  }
3509  }
3510 
3511 }
3512 
3513 TPolyLine* PndLmdDim::Get_Sensor_Shape(int ihalf, int iplane, int imodule,
3514  int iside, int idie, int isensor, bool aligned, bool lmd_frame) {
3515  Double_t x[5] = { -maps_width + 2 * maps_passive_left, maps_width
3516  - 2 * maps_passive_right, maps_width - 2 * maps_passive_right, -maps_width
3517  + 2 * maps_passive_left, -maps_width + 2 * maps_passive_right };
3518  Double_t y[5] = { -maps_height + 2 * maps_passive_bottom, -maps_height
3519  + 2 * maps_passive_bottom, maps_height - 2 * maps_passive_top, maps_height
3520  - 2 * maps_passive_top, -maps_height + 2 * maps_passive_bottom };
3521  for (unsigned int ipoint = 0; ipoint < 5; ipoint++) {
3522  TVector3 point(x[ipoint], y[ipoint], 0);
3523  TVector3 point_master;
3524  if (lmd_frame) {
3525  point_master = Transform_sensor_to_lmd_local(point, ihalf, iplane,
3526  imodule, iside, idie, isensor, false, aligned);
3527  } else {
3528  point_master = Transform_sensor_to_global(point, ihalf, iplane, imodule,
3529  iside, idie, isensor, false, aligned);
3530  }
3531  x[ipoint] = point_master.X();
3532  y[ipoint] = point_master.Y();
3533  }
3534  TPolyLine* pline = new TPolyLine(5, x, y);
3535  //pline->SetFillColor(38);
3536  pline->SetLineColor(2);
3537  pline->SetLineWidth(1);
3538  //pline->Draw("f");
3539  //pline->Draw();
3540  //gPad->Update();
3541  return pline;
3542 }
3543 
3544 vector<TGraph*> PndLmdDim::Get_Sensor_Graph(int ihalf, int iplane, int imodule,
3545  int iside, int idie, int isensor, bool aligned, bool lmd_frame,
3546  bool pixel_subdivision) {
3547  //stringstream name;
3548  //name << "multigraph_half_" << ihalf << "_plane_" << iplane << "_module_" << imodule << "_side_" << iside << "_die_" << idie << "_sensor_" << isensor;
3549  //stringstream title;
3550  //title << "multi graph for sensor " << Get_sensor_id(ihalf, iplane, imodule, iside, idie, isensor);
3551  vector<TGraph*> result; // = new TMultiGraph(name.str().c_str(), title.str().c_str());
3552  vector<double*> xs;
3553  vector<double*> ys;
3554  // create the passive area
3555  // top
3556  Double_t x_top[5] = { -maps_width, maps_width, maps_width, -maps_width,
3557  -maps_width };
3558  Double_t y_top[5] = { maps_height - maps_passive_top * 2., maps_height
3559  - maps_passive_top * 2., maps_height, maps_height, maps_height
3560  - maps_passive_top * 2. };
3561  xs.push_back(x_top);
3562  ys.push_back(y_top);
3563  // bottom
3564  Double_t x_bottom[5] = { -maps_width, maps_width, maps_width, -maps_width,
3565  -maps_width };
3566  Double_t y_bottom[5] = { -maps_height + maps_passive_bottom * 2., -maps_height
3567  + maps_passive_bottom * 2., -maps_height, -maps_height, -maps_height
3568  + maps_passive_bottom * 2. };
3569  xs.push_back(x_bottom);
3570  ys.push_back(y_bottom);
3571  // left
3572  Double_t x_left[5] = { -maps_width, -maps_width + maps_passive_left * 2,
3573  -maps_width + maps_passive_left * 2, -maps_width, -maps_width };
3574  Double_t y_left[5] = { -maps_height + maps_passive_bottom * 2., -maps_height
3575  + maps_passive_bottom * 2., maps_height - maps_passive_top * 2.,
3576  maps_height - maps_passive_top * 2., -maps_height
3577  + maps_passive_bottom * 2. };
3578  xs.push_back(x_left);
3579  ys.push_back(y_left);
3580  // right
3581  Double_t x_right[5] = { +maps_width - maps_passive_right * 2., +maps_width,
3582  +maps_width, +maps_width - maps_passive_right * 2., +maps_width
3583  - maps_passive_right * 2. };
3584  Double_t y_right[5] = { -maps_height + maps_passive_bottom * 2., -maps_height
3585  + maps_passive_bottom * 2., maps_height - maps_passive_top * 2.,
3586  maps_height - maps_passive_top * 2., -maps_height
3587  + maps_passive_bottom * 2. };
3588  xs.push_back(x_right);
3589  ys.push_back(y_right);
3590  // the pixels
3591  int ncols = (int) floor(maps_active_width * 2. / maps_active_pixel_size);
3592  int nrows = (int) floor(maps_active_height * 2. / maps_active_pixel_size);
3593  if (!pixel_subdivision) {
3594  // create only one bin for the active area
3595  TPolyLine* shape = Get_Sensor_Shape(ihalf, iplane, imodule, iside, idie,
3596  isensor, aligned, lmd_frame);
3597  double* x = shape->GetX();
3598  double* y = shape->GetY();
3599  double n = shape->GetN();
3600  result.push_back(new TGraph(n, x, y));
3601 
3602  } else {
3603  cout << "creating " << nrows << " rows x " << ncols << " cols of pixels "
3604  << endl;
3605  for (int irow = 0; irow < nrows; irow++) {
3606  //cout << ".";
3607  //flush(cout);
3608  for (int icol = 0; icol < ncols; icol++) {
3609  // put it on the stack instead of the heap
3610  // to have it still outside this loop scope
3611  double* x = new double[5];
3612  double* y = new double[5];
3613  x[0] = -maps_width + maps_passive_left * 2.
3614  + icol * maps_active_pixel_size;
3615  x[1] = x[0] + maps_active_pixel_size;
3616  x[2] = x[1];
3617  x[3] = x[0];
3618  x[4] = x[0];
3619  xs.push_back(x);
3620  y[0] = -maps_height + maps_passive_bottom * 2.
3621  + irow * maps_active_pixel_size;
3622  y[1] = y[0];
3623  y[2] = y[0] + maps_active_pixel_size;
3624  y[3] = y[2];
3625  y[4] = y[0];
3626  ys.push_back(y);
3627  }
3628  }
3629  }
3630  // now transform (if requested) into the corresponding reference frame
3631  //cout << " transforming " << endl;
3632  for (unsigned int ipoints = 0; ipoints < xs.size(); ipoints++) {
3633  //cout << ".";
3634  //flush(cout);
3635  for (unsigned int ipoint = 0; ipoint < 5; ipoint++) {
3636  TVector3 point(xs[ipoints][ipoint], ys[ipoints][ipoint], 0);
3637  TVector3 point_master;
3638  if (lmd_frame) {
3639  point_master = Transform_sensor_to_lmd_local(point, ihalf, iplane,
3640  imodule, iside, idie, isensor, false, aligned);
3641  } else {
3642  point_master = Transform_sensor_to_global(point, ihalf, iplane, imodule,
3643  iside, idie, isensor, false, aligned);
3644  }
3645  xs[ipoints][ipoint] = point_master.X();
3646  ys[ipoints][ipoint] = point_master.Y();
3647  }
3648  result.push_back(new TGraph(5, xs[ipoints], ys[ipoints]));
3649  // clean up
3650  if (ipoints > 3) {
3651  delete[] xs[ipoints];
3652  delete[] ys[ipoints];
3653  }
3654  }
3655  return result;
3656 }
3657 
3658 TH2Poly* PndLmdDim::Get_histogram_Plane(int iplane, int iside, bool aligned,
3659  bool lmd_frame, bool pixel_subdivision) {
3660  stringstream name;
3661  name << "histpoly_plane_" << iplane << "_side_" << iside;
3662  stringstream title;
3663  title << "plane " << iplane << " side " << iside;
3664  TH2Poly* result = new TH2Poly();
3665  result->SetNameTitle(name.str().c_str(), title.str().c_str());
3666  result->SetFloat();
3667  result->SetXTitle("x [cm]");
3668  result->SetYTitle("y [cm]");
3669 
3670  for (unsigned int ihalf = 0; ihalf < 2; ihalf++)
3671  //for (unsigned int iplane = 0; iplane < n_planes; iplane++)
3672  for (unsigned int imodule = 0; imodule < nmodules; imodule++)
3673  //for (unsigned int iside = 0; iside < 2; iside++)
3674  for (unsigned int idie = 0; idie < 2; idie++)
3675  for (unsigned int isensor = 0; isensor < 3; isensor++) {
3676  if (idie == 1 && isensor == 0)
3677  continue;
3678  vector<TGraph*> sensor_graph = Get_Sensor_Graph(ihalf, iplane,
3679  imodule, iside, idie, isensor, aligned, lmd_frame,
3680  pixel_subdivision);
3681  for (size_t igraph = 0; igraph < sensor_graph.size(); igraph++) {
3682  result->AddBin(sensor_graph[igraph]);
3683  //delete sensor_graph[igraph]; is owned by a PolyBin, so don't delete it
3684  }
3685  }
3686  return result;
3687 }
3688 
3689 TH2Poly* PndLmdDim::Get_histogram_Moduleside(int ihalf, int iplane, int imodule,
3690  int iside, bool aligned, bool lmd_frame, bool pixel_subdivision) {
3691  stringstream name;
3692  name << "histpoly_half_" << ihalf << "_plane_" << iplane << "_module_"
3693  << imodule << "_side_" << iside;
3694  stringstream title;
3695  title << "half " << ihalf << " plane " << iplane << " module " << imodule
3696  << " side " << iside;
3697  TH2Poly* result = new TH2Poly();
3698  result->SetNameTitle(name.str().c_str(), title.str().c_str());
3699  result->SetFloat();
3700  result->SetXTitle("x [cm]");
3701  result->SetYTitle("y [cm]");
3702 
3703  //for (unsigned int ihalf = 0; ihalf < 2; ihalf++)
3704  //for (unsigned int iplane = 0; iplane < n_planes; iplane++)
3705  //for (unsigned int imodule = 0; imodule < nmodules; imodule++)
3706  //for (unsigned int iside = 0; iside < 2; iside++)
3707  for (unsigned int idie = 0; idie < 2; idie++)
3708  for (unsigned int isensor = 0; isensor < 3; isensor++) {
3709  if (idie == 1 && isensor == 0)
3710  continue;
3711  vector<TGraph*> sensor_graph = Get_Sensor_Graph(ihalf, iplane, imodule,
3712  iside, idie, isensor, aligned, lmd_frame, pixel_subdivision);
3713  for (size_t igraph = 0; igraph < sensor_graph.size(); igraph++) {
3714  result->AddBin(sensor_graph[igraph]);
3715  //delete sensor_graph[igraph]; is owned by a PolyBin, so don't delete it
3716  }
3717  }
3718  return result;
3719 }
3720 
3721 TH2Poly* PndLmdDim::Get_histogram_Sensor(int ihalf, int iplane, int imodule,
3722  int iside, int idie, int isensor, bool aligned, bool lmd_frame) {
3723  stringstream name;
3724  name << "histpoly_sensor_half_" << ihalf << "_plane_" << iplane << "_module_"
3725  << imodule << "_side_" << iside << "_die_" << idie << "_sensor_"
3726  << isensor;
3727  stringstream title;
3728  title << "sensor "
3729  << Get_sensor_id(ihalf, iplane, imodule, iside, idie, isensor);
3730  TH2Poly* result = new TH2Poly();
3731  result->SetNameTitle(name.str().c_str(), title.str().c_str());
3732  result->SetFloat();
3733  result->SetXTitle("x [cm]");
3734  result->SetYTitle("y [cm]");
3735 
3736  vector<TGraph*> sensor_graph = Get_Sensor_Graph(ihalf, iplane, imodule, iside,
3737  idie, isensor, aligned, lmd_frame, true);
3738  for (size_t igraph = 0; igraph < sensor_graph.size(); igraph++) {
3739  //cout << " adding bins " << endl;
3740  result->AddBin(sensor_graph[igraph]);
3741  //delete sensor_graph[igraph]; is owned by a PolyBin, so don't delete it
3742  }
3743  return result;
3744 }
3745 
3746 int PndLmdDim::makeOverlapID(int firstSensorId, int secondSensorId) {
3747 
3748  int fhalf, fplane, fmodule, fside, fdie, fsensor;
3749  int bhalf, bplane, bmodule, bside, bdie, bsensor;
3750 
3751  int smalloverlap = 0;
3752 
3753  Get_sensor_by_id(firstSensorId, fhalf, fplane, fmodule, fside, fdie, fsensor);
3754  Get_sensor_by_id(secondSensorId, bhalf, bplane, bmodule, bside, bdie,
3755  bsensor);
3756 
3757  if (fside == bside) {
3758  return -1;
3759  }
3760 
3761  //should return the same id for (id1, id2) as for (id2, id1)
3762  if (bside < fside) {
3763  std::swap(firstSensorId, secondSensorId);
3764  Get_sensor_by_id(firstSensorId, fhalf, fplane, fmodule, fside, fdie,
3765  fsensor);
3766  Get_sensor_by_id(secondSensorId, bhalf, bplane, bmodule, bside, bdie,
3767  bsensor);
3768  }
3769 
3770  if (bhalf != fhalf) {
3771  return -1;
3772  }
3773  if (bplane != fplane) {
3774  return -1;
3775  }
3776  if (bmodule != fmodule) {
3777  return -1;
3778  }
3779 
3780  //0to5
3781  if (fdie == 0 && fsensor == 0 && bdie == 0 && bsensor == 0) {
3782  smalloverlap = 0;
3783  }
3784  //3to8
3785  else if (fdie == 1 && fsensor == 1 && bdie == 1 && bsensor == 1) {
3786  smalloverlap = 1;
3787  }
3788  //4to9
3789  else if (fdie == 1 && fsensor == 2 && bdie == 1 && bsensor == 2) {
3790  smalloverlap = 2;
3791  }
3792  //3to6
3793  else if (fdie == 1 && fsensor == 1 && bdie == 0 && bsensor == 1) {
3794  smalloverlap = 3;
3795  }
3796  //1to8
3797  else if (fdie == 0 && fsensor == 1 && bdie == 1 && bsensor == 1) {
3798  smalloverlap = 4;
3799  }
3800  //2to8
3801  else if (fdie == 0 && fsensor == 2 && bdie == 1 && bsensor == 1) {
3802  smalloverlap = 5;
3803  }
3804  //2to9
3805  else if (fdie == 0 && fsensor == 2 && bdie == 1 && bsensor == 2) {
3806  smalloverlap = 6;
3807  }
3808  //3to7
3809  else if (fdie == 1 && fsensor == 1 && bdie == 0 && bsensor == 2) {
3810  smalloverlap = 7;
3811  }
3812  //4to7
3813  else if (fdie == 1 && fsensor == 2 && bdie == 0 && bsensor == 2) {
3814  smalloverlap = 8;
3815  } else {
3816  //all other checks are negative? then the sensors don't overlap!
3817  return -1;
3818  }
3819  return 1000 * fhalf + 100 * fplane + 10 * fmodule + smalloverlap;
3820 }
3821 
3822 int PndLmdDim::getID1fromOverlapID(int overlapID) {
3823 
3824  int fhalf, fplane, fmodule, fside, fdie, fsensor;
3825  //int bhalf, bplane, bmodule, bside, bdie, bsensor; //[R.K.03/2017] unused variable
3826 
3827  //get info from overlapID
3828  //pad overlapID to 4 digits
3829 
3830  fhalf = std::floor(overlapID / 1000); // = bhalf //[R.K.03/2017] unused variable
3831  fplane = std::floor((overlapID % 1000) / 100); // = bplane //[R.K.03/2017] unused variable
3832  fmodule = std::floor((overlapID % 100) / 10);// = bmodule //[R.K.03/2017] unused variable
3833  int smalloverlap = std::floor(overlapID % 10);
3834 
3835  fside = 0;
3836  //bside = 1; //[R.K.03/2017] unused variable
3837 
3838  if (smalloverlap == 0) {
3839  fdie = 0;
3840  fsensor = 0;
3841  //bdie = 0; //[R.K.03/2017] unused variable
3842  //bsensor = 0; //[R.K.03/2017] unused variable
3843  } else if (smalloverlap == 1) {
3844  fdie = 1;
3845  fsensor = 1;
3846  //bdie = 1; //[R.K.03/2017] unused variable
3847  //bsensor = 1; //[R.K.03/2017] unused variable
3848  } else if (smalloverlap == 2) {
3849  fdie = 1;
3850  fsensor = 2;
3851  //bdie = 1; //[R.K.03/2017] unused variable
3852  //bsensor = 2; //[R.K.03/2017] unused variable
3853  } else if (smalloverlap == 3) {
3854  fdie = 1;
3855  fsensor = 1;
3856  //bdie = 0; //[R.K.03/2017] unused variable
3857  //bsensor = 1; //[R.K.03/2017] unused variable
3858  } else if (smalloverlap == 4) {
3859  fdie = 0;
3860  fsensor = 1;
3861  //bdie = 1; //[R.K.03/2017] unused variable
3862  //bsensor = 1; //[R.K.03/2017] unused variable
3863  } else if (smalloverlap == 5) {
3864  fdie = 0;
3865  fsensor = 2;
3866  //bdie = 1; //[R.K.03/2017] unused variable
3867  //bsensor = 1; //[R.K.03/2017] unused variable
3868  } else if (smalloverlap == 6) {
3869  fdie = 0;
3870  fsensor = 2;
3871  //bdie = 1; //[R.K.03/2017] unused variable
3872  //bsensor = 2; //[R.K.03/2017] unused variable
3873  } else if (smalloverlap == 7) {
3874  fdie = 1;
3875  fsensor = 1;
3876  //bdie = 0; //[R.K.03/2017] unused variable
3877  //bsensor = 2; //[R.K.03/2017] unused variable
3878  } else if (smalloverlap == 8) {
3879  fdie = 1;
3880  fsensor = 2;
3881  //bdie = 0; //[R.K.03/2017] unused variable
3882  //bsensor = 2; //[R.K.03/2017] unused variable
3883  } else {
3884  return -1;
3885  }
3886 
3887  int id = Get_sensor_id(fhalf, fplane, fmodule, fside, fdie, fsensor);
3888  return id;
3889 }
3890 
3891 int PndLmdDim::getID2fromOverlapID(int overlapID) {
3892  //int fhalf, fplane, fmodule, fside, fdie, fsensor; //[R.K.03/2017] unused variable
3893  int bhalf, bplane, bmodule, bside, bdie, bsensor;
3894 
3895  //get info from overlapID
3896  //pad overlapID to 4 digits
3897 
3898  bhalf = std::floor(overlapID / 1000); //fhalf = //[R.K.03/2017] unused variable
3899  bplane = std::floor((overlapID % 1000) / 100); //fplane = //[R.K.03/2017] unused variable
3900  bmodule = std::floor((overlapID % 100) / 10); //fmodule = //[R.K.03/2017] unused variable
3901  int smalloverlap = std::floor(overlapID % 10);
3902 
3903  //fside = 0; //[R.K.03/2017] unused variable
3904  bside = 1;
3905 
3906  if (smalloverlap == 0) {
3907  //fdie = 0; //[R.K.03/2017] unused variable
3908  //fsensor = 0; //[R.K.03/2017] unused variable
3909  bdie = 0;
3910  bsensor = 0;
3911  } else if (smalloverlap == 1) {
3912  //fdie = 1; //[R.K.03/2017] unused variable
3913  //fsensor = 1; //[R.K.03/2017] unused variable
3914  bdie = 1;
3915  bsensor = 1;
3916  } else if (smalloverlap == 2) {
3917  //fdie = 1; //[R.K.03/2017] unused variable
3918  //fsensor = 2; //[R.K.03/2017] unused variable
3919  bdie = 1;
3920  bsensor = 2;
3921  } else if (smalloverlap == 3) {
3922  //fdie = 1; //[R.K.03/2017] unused variable
3923  //fsensor = 1; //[R.K.03/2017] unused variable
3924  bdie = 0;
3925  bsensor = 1;
3926  } else if (smalloverlap == 4) {
3927  //fdie = 0; //[R.K.03/2017] unused variable
3928  //fsensor = 1; //[R.K.03/2017] unused variable
3929  bdie = 1;
3930  bsensor = 1;
3931  } else if (smalloverlap == 5) {
3932  //fdie = 0; //[R.K.03/2017] unused variable
3933  //fsensor = 2; //[R.K.03/2017] unused variable
3934  bdie = 1;
3935  bsensor = 1;
3936  } else if (smalloverlap == 6) {
3937  //fdie = 0; //[R.K.03/2017] unused variable
3938  //fsensor = 2; //[R.K.03/2017] unused variable
3939  bdie = 1;
3940  bsensor = 2;
3941  } else if (smalloverlap == 7) {
3942  //fdie = 1; //[R.K.03/2017] unused variable
3943  //fsensor = 1; //[R.K.03/2017] unused variable
3944  bdie = 0;
3945  bsensor = 2;
3946  } else if (smalloverlap == 8) {
3947  //fdie = 1; //[R.K.03/2017] unused variable
3948  //fsensor = 2; //[R.K.03/2017] unused variable
3949  bdie = 0;
3950  bsensor = 2;
3951  } else {
3952  return -1;
3953  }
3954 
3955  int id = Get_sensor_id(bhalf, bplane, bmodule, bside, bdie, bsensor);
3956  return id;
3957 }
3958 
3959 int PndLmdDim::makeModuleID(int overlapID) {
3960  int moduleID = std::floor(overlapID / 10.0);
3961  return moduleID;
3962 }
3963 
3964 std::vector<int> PndLmdDim::getAvailableOverlapIDs() {
3965  std::vector<int> result;
3966  int overlapID;
3967 
3968  for (int iHalf = 0; iHalf < 2; iHalf++) {
3969  for (int iPlane = 0; iPlane < n_planes; iPlane++) {
3970  for (int iModule = 0; iModule < nmodules; iModule++) {
3971  for (int iOverlap = 0; iOverlap < 9; iOverlap++) {
3972  overlapID = 1000 * iHalf + 100 * iPlane + 10 * iModule + iOverlap;
3973  result.push_back(overlapID);
3974  }
3975  }
3976  }
3977  }
3978  return result;
3979 }
3980 
3981 //
3982 
row
int row
Definition: anaLmdDigi.C:67
PndLmdDim::maps_active_pixel_size
double maps_active_pixel_size
Definition: PndLmdDim.h:406
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
PndLmdDim::sensor_tilt_psi
double sensor_tilt_psi
Definition: PndLmdDim.h:452
PndLmdDim::Draw_Sensors
void Draw_Sensors(int iplane, bool aligned=true, bool lmd_frame=true, int glside=2)
Definition: PndLmdDim.cxx:3462
PndLmdDim::rot_z
double rot_z
Definition: PndLmdDim.h:498
PndLmdDim::die_tilt_psi
double die_tilt_psi
Definition: PndLmdDim.h:439
PndLmdDim::Get_offset
void Get_offset(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, double &x, double &y, double &z, double &rotphi, double &rottheta, double &rotpsi, bool random=false)
Definition: PndLmdDim.cxx:2352
PndLmdDim::Test_List_of_Sensors
bool Test_List_of_Sensors(vector< string > list_of_sensors, int &offset)
Definition: PndLmdDim.cxx:2108
PndLmdDim::pos_y
double pos_y
Definition: PndLmdDim.h:495
PndLmdDim::r_bend
double r_bend
Definition: PndLmdDim.h:484
PndLmdDim::itoffset
map< Tkey, vector< double > >::iterator itoffset
Definition: PndLmdDim.h:559
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndLmdDim::Transform_global_to_lmd_local_vect
void Transform_global_to_lmd_local_vect(double &x, double &y, double &z, bool aligned=true)
Definition: PndLmdDim.cxx:2712
PndLmdDim::Get_pos_lmd_global
void Get_pos_lmd_global(double &x, double &y, double &z, double &rotx, double &roty, double &rotz, bool misaligned=false)
Definition: PndLmdDim.h:726
PndLmdDim::Get_histogram_Plane
TH2Poly * Get_histogram_Plane(int iplane, int iside, bool aligned=true, bool lmd_frame=true, bool pixel_subdivision=false)
Definition: PndLmdDim.cxx:3658
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
PndLmdDim::makeModuleID
int makeModuleID(int overlapID)
Definition: PndLmdDim.cxx:3959
PndLmdDim::Propagate_fast_ip_to_lmd
void Propagate_fast_ip_to_lmd(TVector3 &pos, TVector3 &mom, double pbeam)
Definition: PndLmdDim.cxx:2547
PndLmdDim::half_tilt_phi
double half_tilt_phi
Definition: PndLmdDim.h:275
FairMediumVacuum
FairGeoMedium * FairMediumVacuum
Definition: createRootGeoFile.C:86
PndLmdDim::plane_half_tilt_theta
double plane_half_tilt_theta
Definition: PndLmdDim.h:290
PndLmdDim::Generate_key
string Generate_key(int ihalf, int iplane, int imodule, int iside, int idie, int isensor)
Definition: PndLmdDim.h:596
PndLmdDim::half_offset_z
double half_offset_z
Definition: PndLmdDim.h:273
PndLmdDim::Get_matrix_global_to_lmd_local
TGeoMatrix * Get_matrix_global_to_lmd_local(bool aligned=true)
Definition: PndLmdDim.cxx:2970
PndLmdDim::Get_matrix_lmd_local_to_module_side
TGeoMatrix * Get_matrix_lmd_local_to_module_side(int ihalf, int iplane, int imodule, int iside, bool aligned=true)
Definition: PndLmdDim.cxx:2978
PndLmdDim::Get_histogram_Sensor
TH2Poly * Get_histogram_Sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true, bool lmd_frame=true)
Definition: PndLmdDim.cxx:3721
PndLmdDim::pos_z
double pos_z
Definition: PndLmdDim.h:490
PndLmdDim::Get_Sensor_Shape
TPolyLine * Get_Sensor_Shape(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true, bool lmd_frame=true)
Definition: PndLmdDim.cxx:3513
PndLmdDim::plane_half_offset_y
double plane_half_offset_y
Definition: PndLmdDim.h:284
geoLoad
FairGeoLoader * geoLoad
Definition: createRootGeoFile.C:16
i
Int_t i
Definition: run_full.C:25
Media
FairGeoMedia * Media
Definition: createRootFscGeometryFile.C:68
file
TFile * file
Definition: drc/old/draw_geom.C:12
PndLmdDim::Get_matrix_difference
bool Get_matrix_difference(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, double &dx, double &dy, double &dz, double &dphi, double &dtheta, double &dpsi)
Definition: PndLmdDim.cxx:2757
gGeoMan
TGeoManager * gGeoMan
Definition: createRootFscGeometryFile.C:85
PndLmdDim::cvd_offset_y
double cvd_offset_y
Definition: PndLmdDim.h:329
PndLmdDim::die_tilt_theta
double die_tilt_theta
Definition: PndLmdDim.h:437
PndLmdAlignPar::GetRotateZ
Double_t GetRotateZ(int j) const
Definition: PndLmdAlignPar.h:35
PndLmdDim::die_offset_y
double die_offset_y
Definition: PndLmdDim.h:431
PndLmdDim::delta_phi
double delta_phi
Definition: PndLmdDim.h:310
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndLmdDim::die_tilt_phi
double die_tilt_phi
Definition: PndLmdDim.h:435
PndLmdDim::Get_matrices
map< Tkey, TGeoMatrix * > * Get_matrices(bool aligned=true)
Definition: PndLmdDim.cxx:2725
PndLmdDim::Instance
static PndLmdDim * Instance()
Definition: PndLmdDim.cxx:249
PndLmdDim::Transform_sensor_aligned_to_sensor
TVector3 Transform_sensor_aligned_to_sensor(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false)
Definition: PndLmdDim.cxx:3278
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndLmdDim::end_seg_upstream
double end_seg_upstream
Definition: PndLmdDim.h:482
PndLmdDim::plane_half_tilt_phi
double plane_half_tilt_phi
Definition: PndLmdDim.h:288
PndLmdDim::Transform_global_to_lmd_local
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
n
int n
Definition: HypStatDecay.cxx:21
PndLmdDim::cvd_disc_rad
double cvd_disc_rad
Definition: PndLmdDim.h:303
PndLmdDim::Get_transformation_sensor_to_sensor_aligned
TGeoHMatrix Get_transformation_sensor_to_sensor_aligned(int ihalf, int iplane, int imodule, int iside, int idie, int isensor)
Definition: PndLmdDim.cxx:3093
PndLmdDim::phi_bend
double phi_bend
Definition: PndLmdDim.h:486
rotation
TGeoRotation rotation
Definition: createRootGeoFile.C:81
PndLmdDim::Transform_sensor_to_module_side
TVector3 Transform_sensor_to_module_side(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3216
PndLmdDim::maps_passive_right
double maps_passive_right
Definition: PndLmdDim.h:402
PndLmdDim::maps_die_height
double maps_die_height
Definition: PndLmdDim.h:418
PndLmdDim::plane_pos_z
double * plane_pos_z
Definition: PndLmdDim.h:266
PndLmdDim::Get_List_of_Sensors
string Get_List_of_Sensors(vector< string > &list_of_sensors, bool found_lmd=false, bool first_call=true)
Definition: PndLmdDim.cxx:2074
PndLmdDim::cvd_tilt_theta
double cvd_tilt_theta
Definition: PndLmdDim.h:335
PndLmdDim::rad_entrance
double rad_entrance
Definition: PndLmdDim.h:467
PndLmdDim::Transform_module_side_to_lmd_local
TVector3 Transform_module_side_to_lmd_local(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3200
PndLmdDim::sensor_offset_x
double sensor_offset_x
Definition: PndLmdDim.h:442
PndLmdDim::sensor_tilt_phi
double sensor_tilt_phi
Definition: PndLmdDim.h:448
mom
Double_t mom
Definition: plot_dirc.C:14
offset
TVector3 offset(2, 0, 0)
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
PndLmdDim::plane_half_tilt_psi
double plane_half_tilt_psi
Definition: PndLmdDim.h:292
PndLmdDim::side_offset_x
double side_offset_x
Definition: PndLmdDim.h:344
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
PndLmdDim::outer_rad
double outer_rad
Definition: PndLmdDim.h:299
PndLmdDim::box_size_y
double box_size_y
Definition: PndLmdDim.h:459
PndLmdDim::Read_transformation_matrices_from_geometry
bool Read_transformation_matrices_from_geometry(bool aligned=true)
Definition: PndLmdDim.cxx:1714
PndLmdDim::Calc_matrix_offsets
void Calc_matrix_offsets()
Definition: PndLmdDim.cxx:2814
PndLmdDim::Transform_sensor_to_sensor_aligned
TVector3 Transform_sensor_to_sensor_aligned(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false)
Definition: PndLmdDim.cxx:3263
PndLmdDim::side_tilt_psi
double side_tilt_psi
Definition: PndLmdDim.h:354
PndLmdDim::Get_sensor_id
int Get_sensor_id(int ihalf, int iplane, int imodule, int iside, int idie, int isensor)
Definition: PndLmdDim.h:519
PndLmdAlignPar::GetRotateY
Double_t GetRotateY(int j) const
Definition: PndLmdAlignPar.h:34
PndLmdDim::maps_passive_bottom
double maps_passive_bottom
Definition: PndLmdDim.h:400
PndLmdDim::Set_offset
void Set_offset(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, double x, double y, double z, double rotphi, double rottheta, double rotpsi)
Definition: PndLmdDim.cxx:2331
PndLmdDim::pinstance
static PndLmdDim * pinstance
Definition: PndLmdDim.h:227
PndLmdDim::sensIDoffset
unsigned int sensIDoffset
Definition: PndLmdDim.h:250
PndLmdDim::pol_side_dist_min
double pol_side_dist_min
Definition: PndLmdDim.h:315
PndLmdDim::sensor_offset_y
double sensor_offset_y
Definition: PndLmdDim.h:444
PndLmdDim::pos_plane_0
double pos_plane_0
Definition: PndLmdDim.h:479
PndLmdDim::maps_width
double maps_width
Definition: PndLmdDim.h:408
PndLmdDim::side_offset_y
double side_offset_y
Definition: PndLmdDim.h:346
PndLmdDim::Transform_lmd_local_to_global
TVector3 Transform_lmd_local_to_global(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3171
Tkey::die
signed char die
Definition: PndLmdDim.h:148
PndLmdDim::cvd_disc_thick_half
double cvd_disc_thick_half
Definition: PndLmdDim.h:305
PndLmdDim::Get_transformation_lmd_local_to_module_side
TGeoHMatrix Get_transformation_lmd_local_to_module_side(int ihalf, int iplane, int imodule, int iside, bool aligned=true)
Definition: PndLmdDim.cxx:2995
FairMediumAir
FairGeoMedium * FairMediumAir
Definition: createRootGeoFile.C:85
geobuild
FairGeoBuilder * geobuild
Definition: createRootFscGeometryFile.C:69
shape
TGeoShape * shape
Definition: createDummySensor.C:18
PndLmdDim::Retrieve_version_number
bool Retrieve_version_number()
Definition: PndLmdDim.cxx:286
PndLmdDim::Generate_rootgeom
void Generate_rootgeom(TGeoVolume &mothervol, bool misaligned=false)
Definition: PndLmdDim.cxx:331
PndLmdAlignPar
Alignment Parameter Class for LMD.
Definition: PndLmdAlignPar.h:17
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndLmdDim::maps_active_offset_x
double maps_active_offset_x
Definition: PndLmdDim.h:411
PndLmdAlignPar::GetRotateX
Double_t GetRotateX(int j) const
Definition: PndLmdAlignPar.h:33
PndLmdDim::Get_transformation_module_side_to_sensor
TGeoHMatrix Get_transformation_module_side_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3004
PndLmdDim::Test_matrices
void Test_matrices()
Definition: PndLmdDim.cxx:3439
PndLmdDim::Get_sensor_by_id
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
Definition: PndLmdDim.h:526
PndLmdDim::kapton_disc_thick_half
double kapton_disc_thick_half
Definition: PndLmdDim.h:341
PndLmdDim::makeOverlapID
int makeOverlapID(int firstSensorId, int secondSensorId)
Definition: PndLmdDim.cxx:3746
PndLmdDim::Get_matrix
TGeoMatrix * Get_matrix(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:2741
dtheta
double dtheta
Definition: anaLmdCluster.C:54
PndLmdDim::rot_y
double rot_y
Definition: PndLmdDim.h:497
PndLmdDim::rad_pipe
double rad_pipe
Definition: PndLmdDim.h:471
PndLmdDim::Transform_lmd_local_to_module_side
TVector3 Transform_lmd_local_to_module_side(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3126
PndLmdDim::end_seg_bend
double end_seg_bend
Definition: PndLmdDim.h:491
PndLmdDim::transformation_matrices_aligned
map< Tkey, TGeoMatrix * > transformation_matrices_aligned
Definition: PndLmdDim.h:636
PndLmdDim::box_thickness
double box_thickness
Definition: PndLmdDim.h:463
PndLmdDim::pol_side_lg_half
double pol_side_lg_half
Definition: PndLmdDim.h:313
PndLmdDim::getID2fromOverlapID
int getID2fromOverlapID(int overlapID)
Definition: PndLmdDim.cxx:3891
PndLmdDim::pos_rib
double pos_rib
Definition: PndLmdDim.h:465
PndLmdDim::transformation_matrices
map< Tkey, TGeoMatrix * > transformation_matrices
Definition: PndLmdDim.h:635
PndLmdDim::Get_transformation_sensor_to_module_side
TGeoHMatrix Get_transformation_sensor_to_module_side(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3056
PndLmdDim::half_tilt_theta
double half_tilt_theta
Definition: PndLmdDim.h:277
z
Double_t z
Definition: createRootFscGeometryFile.C:146
PndLmdDim::Set_sensIDoffset
bool Set_sensIDoffset(int offset=-1)
Definition: PndLmdDim.cxx:2161
PndLmdDim::Get_transformation_lmd_local_to_sensor
TGeoHMatrix Get_transformation_lmd_local_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3029
PndLmdDim::cvd_offset_x
double cvd_offset_x
Definition: PndLmdDim.h:327
PndLmdDim::Get_transformation_sensor_to_global
TGeoHMatrix Get_transformation_sensor_to_global(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3065
PndLmdDim::maps_passive_left
double maps_passive_left
Definition: PndLmdDim.h:401
PndLmdDim::Transform_global_to_sensor
TVector3 Transform_global_to_sensor(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3156
PndLmdDim::maps_active_offset_y
double maps_active_offset_y
Definition: PndLmdDim.h:412
PndLmdDim::cvd_tilt_phi
double cvd_tilt_phi
Definition: PndLmdDim.h:333
PndLmdDim::getID1fromOverlapID
int getID1fromOverlapID(int overlapID)
Definition: PndLmdDim.cxx:3822
PndLmdDim::box_size_z
double box_size_z
Definition: PndLmdDim.h:461
PndLmdDim::Get_transformation_global_to_sensor
TGeoHMatrix Get_transformation_global_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3013
name
TString name
Definition: createRootGeoFile.C:82
PndLmdDim::Get_histogram_Moduleside
TH2Poly * Get_histogram_Moduleside(int ihalf, int iplane, int imodule, int iside, bool aligned=true, bool lmd_frame=true, bool pixel_subdivision=true)
Definition: PndLmdDim.cxx:3689
PndLmdDim::die_offset_x
double die_offset_x
Definition: PndLmdDim.h:429
Tkey::sensor
signed char sensor
Definition: PndLmdDim.h:149
PndLmdDim::Read_DB_offsets
void Read_DB_offsets(PndLmdAlignPar *lmdalignpar)
Definition: PndLmdDim.cxx:2249
PndLmdDim::cvd_disc_even_odd_offset
double cvd_disc_even_odd_offset
Definition: PndLmdDim.h:308
PndLmdDim::maps_n_col
int maps_n_col
Definition: PndLmdDim.h:390
PndLmdDim::Write_transformation_matrices
void Write_transformation_matrices(string filename, bool aligned=true, int version_number=geometry_version)
Definition: PndLmdDim.cxx:2179
PndLmdDim::reCreate_transformation_matrices
void reCreate_transformation_matrices()
Definition: PndLmdDim.cxx:1324
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
PndLmdDim::pos_rot_z
double pos_rot_z
Definition: PndLmdDim.h:488
PndLmdDim::Read_transformation_matrices
void Read_transformation_matrices(string filename="", bool aligned=true, int version_number=geometry_version)
Definition: PndLmdDim.cxx:1515
PndLmdDim::geometry_version
static int geometry_version
Definition: PndLmdDim.h:225
PndLmdDim::side_tilt_theta
double side_tilt_theta
Definition: PndLmdDim.h:352
PndLmdDim::cvd_disc_dist
double cvd_disc_dist
Definition: PndLmdDim.h:319
PndLmdDim::die_gap
double die_gap
Definition: PndLmdDim.h:415
PndLmdDim::offsets
map< Tkey, vector< double > > offsets
Definition: PndLmdDim.h:558
PndLmdDim::maps_n_row
int maps_n_row
Definition: PndLmdDim.h:391
PndLmdDim::enabled
bool ** enabled
Definition: PndLmdDim.h:393
PndLmdDim::cvd_offset_z
double cvd_offset_z
Definition: PndLmdDim.h:331
PndLmdDim::inner_rad
double inner_rad
Definition: PndLmdDim.h:297
PndLmdDim::Get_transformation_sensor_to_lmd_local
TGeoHMatrix Get_transformation_sensor_to_lmd_local(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:3082
PndLmdDim::fgGeoMan
TGeoManager * fgGeoMan
Definition: PndLmdDim.h:228
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndLmdDim::maps_die_width
double maps_die_width
Definition: PndLmdDim.h:417
PndLmdDim::Cleanup
void Cleanup()
Definition: PndLmdDim.cxx:2234
PndLmdDim::Set_matrix
bool Set_matrix(string path, TGeoHMatrix *matrix, int ihalf=-1, int iplane=-1, int imodule=-1, int iside=-1, int idie=-1, int isensor=-1)
Definition: PndLmdDim.cxx:1670
PndLmdDim::plane_half_offset_x
double plane_half_offset_x
Definition: PndLmdDim.h:282
PndLmdDim::Get_instance
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
PndLmdDim::maps_active_height
double maps_active_height
Definition: PndLmdDim.h:404
PndLmdDim::n_cvd_discs
int n_cvd_discs
Definition: PndLmdDim.h:301
PndLmdDim::it_transformation_matrices
map< Tkey, TGeoMatrix * >::iterator it_transformation_matrices
Definition: PndLmdDim.h:637
rot
TGeoRotation rot
Definition: createRootFscGeometryFile.C:124
PndLmdDim::Get_transformation_global_to_lmd_local
TGeoHMatrix Get_transformation_global_to_lmd_local(bool aligned=true)
Definition: PndLmdDim.cxx:2988
PndLmdDim::Is_on_Sensor
bool Is_on_Sensor(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:2670
PndLmdDim::Transform_sensor_to_lmd_local
TVector3 Transform_sensor_to_lmd_local(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3248
Tkey::module
signed char module
Definition: PndLmdDim.h:146
PndLmdDim::Get_transformation_sensor_aligned_to_sensor
TGeoHMatrix Get_transformation_sensor_aligned_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor)
Definition: PndLmdDim.cxx:3103
PndLmdDim::pi
double pi
Definition: PndLmdDim.h:260
PndLmdDim::Get_overlapping_sensor
bool Get_overlapping_sensor(const TVector3 &point, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor, bool aligned=true)
Definition: PndLmdDim.cxx:2572
PndLmdDim::Get_transformation_lmd_local_to_global
TGeoHMatrix Get_transformation_lmd_local_to_global(bool aligned=true)
Definition: PndLmdDim.cxx:3040
PndLmdDim::Correct_transformation_matrices
void Correct_transformation_matrices()
Definition: PndLmdDim.cxx:1434
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndLmdDim::nmodules
unsigned int nmodules
Definition: PndLmdDim.h:264
geoFace
FairGeoInterface * geoFace
Definition: createRootGeoFile.C:17
PndLmdDim::half_offset_x
double half_offset_x
Definition: PndLmdDim.h:269
y
Double_t y
Definition: createRootFscGeometryFile.C:159
PndLmdDim::plane_half_offset_z
double plane_half_offset_z
Definition: PndLmdDim.h:286
PndLmdDim::half_tilt_psi
double half_tilt_psi
Definition: PndLmdDim.h:279
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
PndLmdDim::half_offset_y
double half_offset_y
Definition: PndLmdDim.h:271
PndLmdDim::rot_x
double rot_x
Definition: PndLmdDim.h:496
PndLmdDim::Write_transformation_matrices_to_geometry
bool Write_transformation_matrices_to_geometry(bool aligned=true)
Definition: PndLmdDim.cxx:1879
PndLmdDim::maps_height
double maps_height
Definition: PndLmdDim.h:409
PndLmdDim::n_sensors
int n_sensors
Definition: PndLmdDim.h:395
Tkey::plane
signed char plane
Definition: PndLmdDim.h:145
PndLmdDim::cvd_tilt_psi
double cvd_tilt_psi
Definition: PndLmdDim.h:337
PndLmdDim::Get_matrix_module_side_to_sensor
TGeoMatrix * Get_matrix_module_side_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true)
Definition: PndLmdDim.cxx:2983
PndLmdDim::nav_paths
vector< string > nav_paths
Definition: PndLmdDim.h:237
PndLmdDim::Get_transformation_module_side_to_lmd_local
TGeoHMatrix Get_transformation_module_side_to_lmd_local(int ihalf, int iplane, int imodule, int iside, bool aligned=true)
Definition: PndLmdDim.cxx:3047
PndLmdDim::side_tilt_phi
double side_tilt_phi
Definition: PndLmdDim.h:350
PndLmdDim::Transform_module_side_to_sensor
TVector3 Transform_module_side_to_sensor(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3141
PndLmdDim::sensor_offset_z
double sensor_offset_z
Definition: PndLmdDim.h:446
PndLmdDim::sensor_tilt_theta
double sensor_tilt_theta
Definition: PndLmdDim.h:450
PndLmdDim::rad_exit
double rad_exit
Definition: PndLmdDim.h:469
PndLmdDim::Transform_sensor_to_sensor
TVector3 Transform_sensor_to_sensor(const TVector3 &point, int ihalf_from, int iplane_from, int imodule_from, int iside_from, int idie_from, int isensor_from, int ihalf_to, int iplane_to, int imodule_to, int iside_to, int idie_to, int isensor_to, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3293
PndLmdAlignPar::GetShiftY
Double_t GetShiftY(int j) const
Definition: PndLmdAlignPar.h:31
PndLmdDim::getAvailableOverlapIDs
std::vector< int > getAvailableOverlapIDs()
Definition: PndLmdDim.cxx:3964
PndLmdDim::maps_thickness
double maps_thickness
Definition: PndLmdDim.h:398
PndLmdDim::propagation_matrices
map< double, PndLmdDimPropMat > propagation_matrices
Definition: PndLmdDim.h:750
Tkey::side
signed char side
Definition: PndLmdDim.h:147
PndLmdAlignPar::GetShiftZ
Double_t GetShiftZ(int j) const
Definition: PndLmdAlignPar.h:32
PndLmdAlignPar::GetShiftX
Double_t GetShiftX(int j) const
Definition: PndLmdAlignPar.h:30
PndLmdDim::maps_active_width
double maps_active_width
Definition: PndLmdDim.h:403
PndLmdDim::Transform_sensor_to_global
TVector3 Transform_sensor_to_global(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3232
PndLmdDim::n_planes
unsigned int n_planes
Definition: PndLmdDim.h:262
PndLmdDim::Transform_lmd_local_to_sensor
TVector3 Transform_lmd_local_to_sensor(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3185
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
PndLmdDim::Get_Sensor_Graph
vector< TGraph * > Get_Sensor_Graph(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned=true, bool lmd_frame=true, bool pixel_subdivision=true)
Definition: PndLmdDim.cxx:3544
PndLmdDim::pos_x
double pos_x
Definition: PndLmdDim.h:493
PndLmdDim::box_size_x
double box_size_x
Definition: PndLmdDim.h:457
PndLmdDim::maps_passive_top
double maps_passive_top
Definition: PndLmdDim.h:399
point
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
filename
const string filename
Definition: Modify_root_geometry.C:23
Tkey::half
signed char half
Definition: PndLmdDim.h:144
PndLmdDim::Decode_hit
TVector3 Decode_hit(const int sensorID, const double column, const double row, const bool aligned=true, bool newVersion=false)
Definition: PndLmdDim.cxx:2503
PndLmdDim::side_offset_z
double side_offset_z
Definition: PndLmdDim.h:348
PndLmdDim::die_offset_z
double die_offset_z
Definition: PndLmdDim.h:433
Generated on Wed Apr 3 2019 08:43:58 for FairRoot/PandaRoot by   doxygen 1.8.5