Belle II Software  release-05-02-19
ECLGeometryPar.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Poyuan Chen *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <ecl/geometry/ECLGeometryPar.h>
12 #include <framework/logging/Logger.h>
13 
14 #include <G4VTouchable.hh>
15 #include <G4PhysicalVolumeStore.hh>
16 #include <G4NavigationHistory.hh>
17 #include <G4Transform3D.hh>
18 #include <G4Point3D.hh>
19 #include <G4Vector3D.hh>
20 #include <iomanip>
21 
22 using namespace std;
23 using namespace Belle2;
24 using namespace ECL;
25 
26 ECLGeometryPar* ECLGeometryPar::m_B4ECLGeometryParDB = 0;
27 
28 class Mapping_t {
29 public:
30  static void Mapping(int id, int& ThetaId, int& PhiId)
31  {
32  ThetaId = m_Theta[((unsigned int)id) >> 4];
33  PhiId = id - m_dTheta[ThetaId] * 16 - ThetaId * 128;
34  }
35 
36  static void Mapping(int id, int& ThetaId, int& PhiId, int& nrep, int& indx)
37  {
38  Mapping(id, ThetaId, PhiId);
39 
40  int off = m_offsets[ThetaId];
41  int i = m_tbl[ThetaId];
42 
43  int r = m_recip[i];
44  int d = m_denom[i];
45 
46  nrep = (PhiId * r) >> m_RECIPROCAL_SHIFT;
47  indx = off + (PhiId - nrep * d);
48  }
49 
50  static int CellID(int ThetaId, int PhiId)
51  {
52  return PhiId + m_dTheta[ThetaId] * 16 + ThetaId * 128;
53  }
54 
55  static int Offset(int ThetaId)
56  {
57  return m_dTheta[ThetaId] + ThetaId * 8;
58  }
59 
60  static int Indx2ThetaId(int indx)
61  {
62  return m_Theta[indx];
63  }
64 
65  static int ThetaId2NCry(int ThetaId) // Theta Id to the number of crystals @ this Id
66  {
67  return m_denom[m_tbl[ThetaId]];
68  }
69 
70 private:
71  static const char m_dTheta[69];
72  static const unsigned char m_Theta[546], m_tbl[69], m_offsets[69];
73 
74  static const unsigned char m_RECIPROCAL_SHIFT = 16;
75  static const unsigned int m_recip[5], m_denom[5];
76 };
77 
78 const unsigned char Mapping_t::m_Theta[546] = {
79  0, 0, 0, 1, 1, 1,
80  2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4,
81  5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
82  8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10,
83  11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12,
84 
85  13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15,
86  16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18,
87  19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 21, 21, 21, 21,
88  22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24,
89  25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27,
90  28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30,
91  31, 31, 31, 31, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 33, 33, 33, 33,
92  34, 34, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36,
93  37, 37, 37, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39,
94  40, 40, 40, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 42, 42,
95  43, 43, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45,
96  46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 48, 48, 48,
97  49, 49, 49, 49, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 51,
98  52, 52, 52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54,
99  55, 55, 55, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 57, 57, 57,
100  58, 58, 58, 58, 58, 58, 58, 58, 58,
101 
102  59, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 60, 60, 60, 60,
103  61, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 62,
104  63, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 64,
105  65, 65, 65, 65, 65, 65,
106  66, 66, 66, 66, 67, 67, 67, 67, 68, 68, 68, 68
107 };
108 
109 const char Mapping_t::m_dTheta[69] = {
110  0, -5, -10, -14, -18, -22, -24, -26, -28, -30, -32, -34, -33, // forward
111 
112  -32, -31, -30, -29, -28, -27, -26, -25, -24,
113  -23, -22, -21, -20, -19, -18, -17, -16, -15,
114  -14, -13, -12, -11, -10, -9, -8, -7, -6,
115  -5, -4, -3, -2, -1, 0, 1, 2, 3,
116  4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
117 
118  14, 15, 16, 14, 12, 10, 8, 6, 2, -2 // backward
119 };
120 
121 const unsigned char Mapping_t::m_tbl[69] = { // pointer to the denominator/reciprocal arrays or how many crystals in a phi sector
122  1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, // forward
123 
124  0, 0, 0, 0, 0, 0, 0, 0, 0,
125  0, 0, 0, 0, 0, 0, 0, 0, 0,
126  0, 0, 0, 0, 0, 0, 0, 0, 0,
127  0, 0, 0, 0, 0, 0, 0, 0, 0,
128  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
129 
130  4, 4, 3, 3, 3, 3, 3, 2, 2, 2 // backward
131 };
132 
133 const unsigned char Mapping_t::m_offsets[69] = {
134  0, 3, 6, 10, 14, 18, 24, 30, 36, 42, 48, 54, 63, // forward
135 
136  132, 134, 136, 138, 140, 142, 144, 146, 148,
137  150, 152, 154, 156, 158, 160, 162, 164, 166,
138  168, 170, 172, 174, 176, 178, 180, 182, 184,
139  186, 188, 190, 192, 194, 196, 198, 200, 202,
140  204, 206, 208, 210, 212, 214, 216, 218, 220, 222,
141 
142  72, 81, 90, 96, 102, 108, 114, 120, 124, 128 // backward
143 };
144 
145 #define pack(x) ((1<<Mapping_t::m_RECIPROCAL_SHIFT)/x+1)
146 const unsigned int Mapping_t::m_recip[] = {pack(2), pack(3), pack(4), pack(6), pack(9)};
147 const unsigned int Mapping_t::m_denom[] = { (2), (3), (4), (6), (9)};
148 #undef pack
149 
150 
151 ECLGeometryPar* ECLGeometryPar::Instance()
152 {
153  if (!m_B4ECLGeometryParDB) m_B4ECLGeometryParDB = new ECLGeometryPar();
154  return m_B4ECLGeometryParDB;
155 }
156 
157 ECLGeometryPar::ECLGeometryPar()
158 {
159  clear();
160  // delay crystal positions fetching to a moment when it actually needed and geometry is already built
161  // read();
162 }
163 
164 ECLGeometryPar::~ECLGeometryPar()
165 {
166  if (m_B4ECLGeometryParDB) {
167  delete m_B4ECLGeometryParDB;
168  B2DEBUG(150, "m_B4ECLGeometryParDB deleted ");
169  }
170  if (m_ECLForwardGlobalT) delete m_ECLForwardGlobalT;
171  if (m_ECLBarrelGlobalT) delete m_ECLBarrelGlobalT;
172  if (m_ECLBackwardGlobalT) delete m_ECLBackwardGlobalT;
173 }
174 
175 void ECLGeometryPar::clear()
176 {
177  m_ini_cid = -1;
178  mPar_cellID = 0;
179  mPar_thetaID = 0;
180  mPar_phiID = 0;
181 }
182 
183 // There is no way to get world coordinates of a local point of a physical volume in Geant.
184 // The only way to check geometry is to trace particle and check volumes it crosses.
185 // Here particle gun parameters are produced to check crystal positions with the center @ r0 and direction n
186 void ParticleGunParameters(const G4String& comment, const G4ThreeVector& n, const G4ThreeVector& r0, double dphi)
187 {
188  cout << comment << endl;
189  cout << "Center position = " << r0 << ", Direction = " << n << endl;
190  // closest point to z-axis
191  double t = -(n.z() * r0.z() - n * r0) / (n.z() * n.z() - n * n);
192  G4ThreeVector r = n * t + r0;
193  cout << "Closest point to z-axis = " << r << endl; // at the moment I do not see tilt in phi
194  const double r2d = 180 / M_PI;
195  double th = r2d * n.theta();
196  double phi = r2d * n.phi();
197  double z = r.z();
198  dphi *= r2d;
199  cout << " 'thetaParams': [" << th << ", " << th << "]," << endl;
200  cout << " 'phiParams': [" << phi << "+0*" << dphi << ", " << phi << "+0*" << dphi << "]," << endl;
201  cout << " 'zVertexParams': [" << z << ", " << z << "]" << endl;
202 }
203 
204 G4Transform3D getT(const G4VPhysicalVolume& v)
205 {
206  return G4Transform3D(v.GetObjectRotationValue(), v.GetObjectTranslation());
207 }
208 
209 G4Transform3D getT(const G4String& n)
210 {
211  G4PhysicalVolumeStore* store = G4PhysicalVolumeStore::GetInstance();
212  // int N = store->size();
213  // cout<<N<<endl;
214  // for(int i=0;i<N;i++){
215  // G4VPhysicalVolume *vv = (*store)[i];
216  // if(vv->GetName()==n){
217  // cout<<i<<" "<<vv->GetName()<<" "<<vv->GetMultiplicity()<<" "<<vv->GetCopyNo()<<endl;
218  // G4PVReplica *r = dynamic_cast<G4PVReplica*>(vv);
219  // if(r) cout<<"Instance = "<<r->GetInstanceID()<<endl;
220  // }
221  // }
222 
223  G4VPhysicalVolume* v = store->GetVolume(n);
224  return getT(*v);
225 }
226 
227 void print_trans(const G4Transform3D& t)
228 {
229  for (int i = 0; i < 4; i++) {
230  for (int j = 0; j < 3; j++)
231  cout << t(j, i) << " ";
232  cout << endl;
233  }
234 }
235 
236 void ECLGeometryPar::read()
237 {
238  m_crystals.clear();
239  m_crystals.reserve(46 * 2 + 132); // 10752 bytes
240 
241  // Endcap sectors are rotated since behaviour of G4Replica class. It
242  // requires a physical volume during phi replication to be symmetric
243  // around phi=0. So we need to rotate it by -((2*pi)/16)/2 angle at the
244  // geometry description are return back here.
245 
246  G4Point3D p0(0, 0, 0); G4Vector3D v0(0, 0, 1);
247 
248  {
249  m_ECLForwardGlobalT = new G4Transform3D(getT("ECLForwardPhysical"));
250  G4Transform3D T = getT("ECLForwardCrystalSectorPhysical_1");
251  G4String tnamef("ECLForwardWrappedCrystal_Physical_");
252  for (int i = 0; i < 72; i++) {
253  G4String vname(tnamef); vname += to_string(i);
254  G4Transform3D cT = T * getT(vname);
255  CrystalGeom_t c = {cT * p0, cT * v0};
256  m_crystals.push_back(c);
257  // G4cout << i << " " << c.pos << " " << c.dir << G4endl;
258  }
259  }
260 
261  {
262  m_ECLBackwardGlobalT = new G4Transform3D(getT("ECLBackwardPhysical"));
263  G4Transform3D T = getT("ECLBackwardCrystalSectorPhysical_0");
264  G4String tnamef("ECLBackwardWrappedCrystal_Physical_");
265  for (int i = 0; i < 60; i++) {
266  G4String vname(tnamef); vname += to_string(i);
267  G4Transform3D cT = T * getT(vname);
268  CrystalGeom_t c = {cT * p0, cT * v0};
269  m_crystals.push_back(c);
270  // G4cout << i << " " << c.pos << " " << c.dir << G4endl;
271  }
272  }
273 
274  {
275  // get barrel sector (between two septums) transformation
276 
277  // extract global barrel movements assuming we know transformation of the sector in axial symmetric situation
278  m_ECLBarrelGlobalT = new G4Transform3D(getT("ECLBarrelSectorPhysical_0") * G4RotateZ3D(M_PI / 2));
279  G4Transform3D Ts = G4RotateZ3D(-M_PI / 2);
280  // since barrel sector is symmetric around phi=0 we need to
281  // translate crystal with negative phi back to positive rotating
282  // crystal position by (2*M_PI/72) angle
283  G4Transform3D Ts1 = G4RotateZ3D(M_PI / 36) * Ts;
284 
285  G4String tname("ECLBarrelWrappedCrystal_Physical_");
286  for (int i = 0; i < 2 * 46; i++) {
287  G4String vname(tname); vname += to_string(i);
288  G4Transform3D cT = ((i % 2) ? Ts1 : Ts) * getT(vname);
289  CrystalGeom_t c = {cT * p0, cT * v0};
290  m_crystals.push_back(c);
291  // G4cout << i << " " << c.pos << " " <<c.dir<<G4endl;
292  }
293  }
294 
295  B2DEBUG(150, "ECLGeometryPar::read() initialized with " << m_crystals.size() << " crystals.");
296 }
297 
298 template <int n>
299 void sincos(const double* ss, int iphi, double& s, double& c)
300 {
301  int n4 = n / 4;
302  int iq = iphi / n4; // quadrant
303  int is = iphi % n4;
304  int ic = n4 - is;
305  if (iq & 1) {int t = is; is = ic; ic = t;}
306  double ls = ss[is];
307  double lc = ss[ic];
308  // check quadrant and assign sign bit accordingly
309  if ((((0 << 0) + (0 << 1) + (1 << 2) + (1 << 3)) >> iq) & 1) ls = -ls;
310  if ((((0 << 0) + (1 << 1) + (1 << 2) + (0 << 3)) >> iq) & 1) lc = -lc;
311  s = ls;
312  c = lc;
313 }
314 
315 const double ss72[] = {
316  0.0000000000000000000000,
317  0.0871557427476581735581,
318  0.1736481776669303488517,
319  0.2588190451025207623489,
320  0.3420201433256687330441,
321  0.4226182617406994361870,
322  0.5000000000000000000000,
323  0.5735764363510460961080,
324  0.6427876096865393263226,
325  0.7071067811865475244008,
326  0.7660444431189780352024,
327  0.8191520442889917896845,
328  0.8660254037844386467637,
329  0.9063077870366499632426,
330  0.9396926207859083840541,
331  0.9659258262890682867497,
332  0.9848077530122080593667,
333  0.9961946980917455322950,
334  1.0000000000000000000000
335 };
336 
337 const double ss16[] = {
338  0.0000000000000000000000,
339  0.3826834323650897717285,
340  0.7071067811865475244008,
341  0.9238795325112867561282,
342  1.0000000000000000000000
343 };
344 
345 void ECLGeometryPar::InitCrystal(int cid)
346 {
347  if (m_crystals.size() == 0) read();
348  int thetaid, phiid, nreplica, indx;
349  Mapping_t::Mapping(cid, thetaid, phiid, nreplica, indx);
350  // cout<<cid<<" "<<thetaid<<" "<<phiid<<" "<<nreplica<<" "<<indx<<endl;
351  const CrystalGeom_t& t = m_crystals[indx];
352  double s, c;
353  if (indx > 131)
354  sincos<72>(ss72, nreplica, s, c);
355  else
356  sincos<16>(ss16, nreplica, s, c);
357 
358  G4Transform3D* T;
359  if (cid < 1152) {
360  T = m_ECLForwardGlobalT;
361  } else if (cid < 7776) {
362  T = m_ECLBarrelGlobalT;
363  } else {
364  T = m_ECLBackwardGlobalT;
365  c = -c;
366  }
367  double xp = c * t.pos.x() - s * t.pos.y();
368  double yp = s * t.pos.x() + c * t.pos.y();
369  // m_current_crystal.pos.set(xp, yp, t.pos.z());
370 
371  double xv = c * t.dir.x() - s * t.dir.y();
372  double yv = s * t.dir.x() + c * t.dir.y();
373  // m_current_crystal.dir.set(xv, yv, t.dir.z());
374 
375  m_current_crystal.pos = (1 / CLHEP::cm) * (*T * G4Point3D(xp, yp, t.pos.z()));
376  m_current_crystal.dir = *T * G4Vector3D(xv, yv, t.dir.z());
377 
378  // cout<<t.pos<<" "<<t.dir<<" "<<m_current_crystal.pos<<" "<<m_current_crystal.dir<<endl;
379  m_ini_cid = cid;
380 }
381 
382 int ECLGeometryPar::GetCellID(int ThetaId, int PhiId)
383 {
384  return Mapping_t::CellID(ThetaId, PhiId);
385 }
386 
387 void ECLGeometryPar::Mapping(int cid)
388 {
389  mPar_cellID = cid;
390  Mapping_t::Mapping(mPar_cellID, mPar_thetaID, mPar_phiID);
391 }
392 
393 int ECLGeometryPar::TouchableDiodeToCellID(const G4VTouchable* touch)
394 {
395  return TouchableToCellID(touch);
396 }
397 
398 int ECLGeometryPar::TouchableToCellID(const G4VTouchable* touch)
399 {
400  // touch->GetCopyNumber() is a virtual call, avoid it by using
401  // directly G4NavigationHistory so we will have only one virtual
402  // call instead of three here
403  const G4NavigationHistory* h = touch->GetHistory();
404  int hd = h->GetDepth();
405  int i1 = h->GetReplicaNo(hd - 1); // index of each volume is set at physical volume creation
406  int i2 = h->GetReplicaNo(hd - 2); // go up in volume hierarchy
407 
408  int ThetaId = Mapping_t::Indx2ThetaId(i1);
409  int NCryst = Mapping_t::ThetaId2NCry(ThetaId); // the number of crystals in a sector at given ThetaId
410  int Offset = Mapping_t::Offset(ThetaId);
411 
412  int ik = i1 - Offset;
413  int PhiId;
414  if (NCryst == 2) {
415  PhiId = (i2 - (ik % 2)) * NCryst + ik;
416  if (PhiId < 0) PhiId += 144;
417  } else {
418  int i3 = h->GetReplicaNo(hd - 3); // go up in volume hierarchy
419  if (ThetaId < 13)
420  PhiId = (i2 + 2 * i3) * NCryst + ik;
421  else {
422  // int tt[] = {3,2,1,0,7,6,5,4};
423  // if(i3r!=tt[i3]){
424  // cout<<i3<<" "<<i3<<" "<<tt[i3]<<endl;
425  // exit(0);
426  // }
427  int i3r = (3 - (i3 % 4)) + 4 * (i3 / 4);
428  PhiId = (2 * i3r + (1 - i2)) * NCryst + ik;
429  }
430  }
431 
432  int cellID = Offset * 16 + PhiId;
433  // cout<<"ECLGeometryPar::TouchableToCellID "<<h->GetVolume(hd-1)->GetName()<<" "<<i1<<" "<<i2<<" "<<h->GetReplicaNo(hd - 3)<<" "<<ThetaId<<" "<<NCryst<<" "<<Offset<<" "<<PhiId<<endl;
434 
435  // test of the position and direction of crystal
436  if (0) {
437  G4AffineTransform t = h->GetTopTransform();
438  G4ThreeVector o(0, 0, 0), n(0, 0, 1);
439  G4ThreeVector ro = t.Inverse().TransformPoint(o);
440  G4ThreeVector rn = t.Inverse().TransformAxis(n);
441 
442  InitCrystal(cellID);
443  ro *= 1 / CLHEP::cm;
444 
445  G4ThreeVector dr = m_current_crystal.pos - ro, dn = m_current_crystal.dir - rn;
446  if (dr.mag() > 1e-10 || dn.mag() > 1e-10) {
447  cout << "Missmatch " << h->GetVolume(hd - 1)->GetName() << " " << cellID << " " << ThetaId << " " << PhiId << " " << NCryst << " "
448  << hd << " " << i2 << " " << i1 << " " << m_current_crystal.pos << " " << ro << " " << rn << " " << dr << " " << dn << endl;
449 
450  for (int i = 0; i < 144; i++) {
451  int ci = Mapping_t::CellID(ThetaId, i);
452  InitCrystal(ci);
453  dr = m_current_crystal.pos - ro;
454  if (dr.mag() < 1e-10) cout << "best PhiId = " << i << endl;
455  }
456  }
457  }
458  return cellID;
459 }
460 
461 int ECLGeometryPar::ECLVolumeToCellID(const G4VTouchable* touch)
462 {
463  int depth = touch->GetHistoryDepth();
464  if ((depth != 3) && (depth != 5)) {
465  B2WARNING("ECLGeometryPar::ECLVolumeToCellID: History depth = " << depth << " is out of range: should be 3 or 5.");
466  return -1;
467  }
468  const G4String& vname = touch->GetVolume()->GetName();
469  std::size_t pos0 = vname.find("lv_forward_crystal_");
470  std::size_t pos1 = vname.find("lv_barrel_crystal_");
471  std::size_t pos2 = vname.find("lv_backward_crystal_");
472  if (pos0 == string::npos && pos1 == string::npos && pos2 == string::npos) {
473  B2WARNING("ECLGeometryPar::ECLVolumeToCellID: Volume name does not match pattern. NAME=" << vname);
474  return -1;
475  }
476  return TouchableToCellID(touch);
477 }
478 
479 double ECLGeometryPar::time2sensor(int cid, const G4ThreeVector& hit_pos)
480 {
481  if (cid != m_ini_cid) InitCrystal(cid);
482  double z = 15. - (hit_pos - m_current_crystal.pos) * m_current_crystal.dir; // position along the vector of crystal axis
483  double dt = 6.05 + z * (0.0749 - z * 0.00112); // flight time to diode sensor in nanoseconds
484  return dt;
485 }
486 
487 EclNbr::EclNbr() :
488  m_nbrs(*new std::vector< Identifier >)
489 {
490  mNbr_cellID = 0;
491  mNbr_thetaID = 0;
492  mNbr_phiID = 0;
493 }
494 
495 EclNbr::EclNbr(const EclNbr& aNbr) :
496  m_nbrs(*new std::vector< Identifier > (aNbr.m_nbrs)) ,
497  m_nearSize(aNbr.m_nearSize)
498 {
499  mNbr_cellID = 0;
500  mNbr_thetaID = 0;
501  mNbr_phiID = 0;
502 }
503 
504 EclNbr::EclNbr(
505  const std::vector< Identifier >& aNbrs ,
506  const std::vector< Identifier >::size_type aNearSize
507 ) :
508  m_nbrs(*new std::vector< Identifier > (aNbrs)) ,
509  m_nearSize(aNearSize)
510 {
511  // sort vector separately for near, nxt-near nbrs
512  std::sort(m_nbrs.begin() , m_nbrs.begin() + aNearSize , std::less< Identifier >()) ;
513  std::sort(m_nbrs.begin() + aNearSize , m_nbrs.end() , std::less< Identifier >()) ;
514 }
515 
516 EclNbr::~EclNbr()
517 {
518  delete &m_nbrs ;
519 }
520 std::ostream& operator<<(std::ostream& os, const EclNbr& aNbr)
521 {
522  os << "N(" ;
523  unsigned short i(0) ;
524  for (std::vector< EclNbr::Identifier >::const_iterator
525  iNbr(aNbr.nbrs().begin()) ;
526  iNbr != aNbr.nbrs().end() ; ++iNbr) {
527  ++i;
528  if (iNbr != aNbr.nbrs().begin() && i != aNbr.nearSize() + 1) os << "," ;
529  if (i == aNbr.nearSize() + 1) os << "|" ;
530  os << std::setw(4) << (*iNbr) ;
531  }
532  os << ")" ;
533  return os ;
534 }
535 
536 
537 void EclNbr::printNbr()
538 {
539  unsigned short Nri(0) ;
540  cout << "(";
541  for (std::vector< EclNbr::Identifier >::const_iterator
542  iNbr(m_nbrs.begin()) ;
543  iNbr != m_nbrs.end() ; ++iNbr) {
544  ++Nri;
545  if (iNbr != m_nbrs.begin() && Nri != m_nearSize + 1) cout << "," ;
546  if (Nri == m_nearSize + 1) cout << "|" ;
547  cout << std::setw(4) << (*iNbr) ;
548  }
549  cout << ")" << endl;
550 }
551 //
552 // assignment operators
553 
554 
555 EclNbr& EclNbr::operator=(const EclNbr& aNbr)
556 {
557  if (this != &aNbr) {
558  mNbr_cellID = aNbr.mNbr_cellID;
559  mNbr_thetaID = aNbr.mNbr_thetaID;
560  mNbr_phiID = aNbr.mNbr_phiID;
561 
562  m_nbrs = aNbr.m_nbrs ;
563  m_nearSize = aNbr.m_nearSize ;
564  }
565  return *this ;
566 }
567 
568 //
569 // member functions
570 //
571 
572 //
573 // const member functions
574 //
575 
576 const std::vector< EclNbr::Identifier >&
577 EclNbr::nbrs() const
578 {
579  return m_nbrs ;
580 }
581 
582 const std::vector< EclNbr::Identifier >::const_iterator
583 EclNbr::nearBegin() const
584 {
585  return m_nbrs.begin() ;
586 }
587 
588 const std::vector< EclNbr::Identifier >::const_iterator
589 EclNbr::nearEnd() const
590 {
591  return m_nbrs.begin() + m_nearSize ;
592 }
593 
594 const std::vector< EclNbr::Identifier >::const_iterator
595 EclNbr::nextBegin() const
596 {
597  return m_nbrs.begin() + m_nearSize ;
598 }
599 
600 const std::vector< EclNbr::Identifier >::const_iterator
601 EclNbr::nextEnd() const
602 {
603  return m_nbrs.end() ;
604 }
605 
606 std::vector< EclNbr::Identifier >::size_type
607 EclNbr::nearSize() const
608 {
609  return m_nearSize ;
610 }
611 
612 std::vector< EclNbr::Identifier >::size_type
613 EclNbr::nextSize() const
614 {
615  return (m_nbrs.size() - m_nearSize) ;
616 }
617 
618 
619 int EclNbr::GetCellID(int ThetaId, int PhiId)
620 {
621  mNbr_cellID = Mapping_t::CellID(ThetaId, PhiId);
622  mNbr_thetaID = ThetaId;
623  mNbr_phiID = PhiId ;
624  return mNbr_cellID;
625 }
626 
627 void EclNbr::Mapping(int cid)
628 {
629  mNbr_cellID = cid;
630  Mapping_t::Mapping(mNbr_cellID, mNbr_thetaID, mNbr_phiID);
631 }
632 
633 EclNbr
634 EclNbr::getNbr(const Identifier aCellId)
635 {
636  // generate nbr lists. always easier here to work with theta-phi
637 
638  const int cellID = aCellId;
639  Mapping(cellID);
640  const int thetaId = GetThetaID();
641  const int phiId = GetPhiID();
642  std::vector< EclNbr::Identifier >::size_type nearSize(0);
643  std::vector< EclNbr::Identifier > vNbr;
644 
645  vNbr.reserve(24) ; // except for extreme endcaps, always 24
646 
647  int t00 = thetaId;
648  int tm1 = thetaId - 1;
649  int tm2 = thetaId - 2;
650  int tp1 = thetaId + 1;
651  int tp2 = thetaId + 2;
652 
653  if (aCellId > 1151 && aCellId < 7776) {
654  // barrel
655  //
656  // 12 13 14 15 16 ^ theta
657  // 11 2 3 4 17 |
658  // 10 1 0 5 18 +--> phi X--+ view from inside
659  // 9 8 7 6 19 | (foot pointing e- dir)
660  // 24 23 22 21 20 Z
661  int f00 = phiId;
662  int fm1 = (phiId + 143) % 144;
663  int fp1 = (phiId + 1) % 144;
664  int fm2 = (phiId + 142) % 144;
665  int fp2 = (phiId + 2) % 144;
666 
667  vNbr.push_back(GetCellID(t00 , fm1));
668  vNbr.push_back(GetCellID(tp1 , fm1));
669  vNbr.push_back(GetCellID(tp1 , f00));
670  vNbr.push_back(GetCellID(tp1 , fp1));
671  vNbr.push_back(GetCellID(t00 , fp1));
672  vNbr.push_back(GetCellID(tm1 , fp1));
673  vNbr.push_back(GetCellID(tm1 , f00));
674  vNbr.push_back(GetCellID(tm1 , fm1));
675 
676  nearSize = vNbr.size();
677 
678  vNbr.push_back(GetCellID(tm1 , fm2));
679  vNbr.push_back(GetCellID(t00 , fm2));
680  vNbr.push_back(GetCellID(tp1 , fm2));
681  vNbr.push_back(GetCellID(tp2 , fm2));
682  vNbr.push_back(GetCellID(tp2 , fm1));
683  vNbr.push_back(GetCellID(tp2 , f00));
684  vNbr.push_back(GetCellID(tp2 , fp1));
685  vNbr.push_back(GetCellID(tp2 , fp2));
686  vNbr.push_back(GetCellID(tp1 , fp2));
687  vNbr.push_back(GetCellID(t00 , fp2));
688  vNbr.push_back(GetCellID(tm1 , fp2));
689  vNbr.push_back(GetCellID(tm2 , fp2));
690  vNbr.push_back(GetCellID(tm2 , fp1));
691  vNbr.push_back(GetCellID(tm2 , f00));
692  vNbr.push_back(GetCellID(tm2 , fm1));
693  vNbr.push_back(GetCellID(tm2 , fm2));
694  }//if( aCellId > 1152 && aCellId < 7777 )
695  else {
696  // endcap -- not always 24!
697  int n00 = 1000;
698  int np1 = 1000;
699  int np2 = 1000;
700  int nm1 = 1000;
701  int nm2 = 1000;
702  if (aCellId < 1153) { // forward
703  const EclIdentifier mPerRingForward[]
704  = { 48, 48, 64, 64, 64, 96, 96, 96, 96, 96, 96, 144, 144, 144, 144 };
705  if (thetaId > 1) nm2 = mPerRingForward[ thetaId - 2 ];
706  if (thetaId > 0) nm1 = mPerRingForward[ thetaId - 1 ];
707  n00 = mPerRingForward[ thetaId ];
708  np1 = mPerRingForward[ thetaId + 1 ];
709  np2 = mPerRingForward[ thetaId + 2 ];
710  } else { // backward
711  const EclIdentifier mPerRingBackward[]
712  = { 64, 64, 64, 96, 96, 96, 96, 96, 144, 144, 144, 144 };
713  if (thetaId < 67) np2 = mPerRingBackward[ 66 - thetaId ];
714  if (thetaId < 68) np1 = mPerRingBackward[ 67 - thetaId ];
715  n00 = mPerRingBackward[ 68 - thetaId ];
716  nm1 = mPerRingBackward[ 69 - thetaId ];
717  nm2 = mPerRingBackward[ 70 - thetaId ];
718  }
719  // f-- are phi's, t-- are thetas
720  // all calculations should be integer arith - pcs
721  // f(th,phi)
722  // criterion: center -> next bin
723  int f0000 = phiId;
724  int fp100 = (f0000 * np1 + np1 / 2) / n00;
725  int fp200 = (f0000 * np2 + np2 / 2) / n00;
726  int fm100 = (f0000 * nm1 + nm1 / 2) / n00;
727  int fm200 = (f0000 * nm2 + nm2 / 2) / n00;
728 
729  int f00m1 = (f0000 + n00 - 1) % n00; // should be exact
730  int f00m2 = (f0000 + n00 - 2) % n00;
731  int f00p1 = (f0000 + 1) % n00;
732  int f00p2 = (f0000 + 2) % n00;
733 
734  int fp1m1 = (fp100 + np1 - 1) % np1;
735  int fp1m2 = (fp100 + np1 - 2) % np1;
736  int fp1p1 = (fp100 + 1) % np1;
737  int fp1p2 = (fp100 + 2) % np1;
738 
739  int fm1m1 = (fm100 + nm1 - 1) % nm1;
740  int fm1m2 = (fm100 + nm1 - 2) % nm1;
741  int fm1p1 = (fm100 + 1) % nm1;
742  int fm1p2 = (fm100 + 2) % nm1;
743 
744  int fp2m1 = (fp200 + np2 - 1) % np2;
745  int fp2m2 = (fp200 + np2 - 2) % np2;
746  int fp2p1 = (fp200 + 1) % np2;
747  int fp2p2 = (fp200 + 2) % np2;
748 
749  int fm2m1 = (fm200 + nm2 - 1) % nm2;
750  int fm2m2 = (fm200 + nm2 - 2) % nm2;
751  int fm2p1 = (fm200 + 1) % nm2;
752  int fm2p2 = (fm200 + 2) % nm2;
753  int delta = n00 / 16;
754 // int sector = phiId/delta; // 0..15
755  int nth = phiId % delta;
756 
757  switch (thetaId) {
758  case 0:
759  if (nth == 1)
760  fp2p2 = 1000;
761  break;
762  case 1:
763  if (nth == 1) {
764  fp2p2 = 1000;
765  fp1p2 = fp1p1;
766  fp1p1 = 1000;
767  }
768  break;
769  case 2:
770  if ((nth == 0) || (nth == 1)) {
771  fm2p2 = 1000;
772  fm1p2 = fm1p1;
773  fm1p1 = 1000;
774  } else if ((nth == 2) || (nth == 3)) {
775  fm2m2 = 1000;
776  fm1m2 = fm1m1;
777  fm1m1 = 1000;
778  }
779  break;
780  case 3:
781  if ((nth == 0) || (nth == 3)) {
782  fm2p2 = fm2m2 = 1000;
783  } else if (nth == 1) {
784  fm2p2 = 1000;
785  } else if (nth == 2) {
786  fm2m2 = 1000;
787  }
788  break;
789  case 5:
790  if ((nth == 2) || (nth == 5)) {
791  fm2m2 = 1000;
792  fm1m2 = fm1m1;
793  fm1m1 = 1000;
794  } else {
795  fm2p2 = 1000;
796  fm1p2 = fm1p1;
797  fm1p1 = 1000;
798  }
799  break;
800  case 6:
801  fm2p2 = 1000;
802  if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5)) {
803  fm2m2 = 1000;
804  }
805  break;
806  case 11:
807  if ((nth == 2) || (nth == 5) || (nth == 8)) {
808  fm2m2 = 1000;
809  fm1m2 = fm1m1;
810  fm1m1 = 1000;
811  } else {
812  fm2p2 = 1000;
813  fm1p2 = fm1p1;
814  fm1p1 = 1000;
815  }
816  break;
817  case 12:
818  fm2p2 = 1000;
819  if ((nth == 0) || (nth == 2) || (nth == 3)
820  || (nth == 5) || (nth == 6) || (nth == 8))
821  fm2m2 = 1000;
822  break;
823  case 65:
824  if ((nth == 2) || (nth == 5)) {
825  fp2m2 = 1000;
826  fp1m2 = fp1m1;
827  fp1m1 = 1000;
828  } else {
829  fp2p2 = 1000;
830  fp1p2 = fp1p1;
831  fp1p1 = 1000;
832  }
833  break;
834  case 64:
835  fp2p2 = 1000;
836  if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5))
837  fp2m2 = 1000;
838  break;
839  case 60:
840  if ((nth == 2) || (nth == 5) || (nth == 8)) {
841  fp2m2 = 1000;
842  fp1m2 = fp1m1;
843  fp1m1 = 1000;
844  } else {
845  fp2p2 = 1000;
846  fp1p2 = fp1p1;
847  fp1p1 = 1000;
848  }
849  break;
850  case 59:
851  fp2p2 = 1000;
852  if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5)
853  || (nth == 6) || (nth == 8))
854  fp2m2 = 1000;
855  break;
856  }//switch
857 
858  // insert near-nbrs
859  vNbr.push_back(GetCellID(t00, f00m1));
860  vNbr.push_back(GetCellID(t00, f00p1));
861  if (nm1 < 999) {
862  vNbr.push_back(GetCellID(tm1 , fm100));
863  if (fm1m1 < 999)
864  vNbr.push_back(GetCellID(tm1 , fm1m1));
865  if (fm1p1 < 999)
866  vNbr.push_back(GetCellID(tm1 , fm1p1));
867  }
868  if (np1 < 999) {
869  vNbr.push_back(GetCellID(tp1 , fp100));
870  if (fp1m1 < 999)
871  vNbr.push_back(GetCellID(tp1 , fp1m1));
872  if (fp1p1 < 999)
873  vNbr.push_back(GetCellID(tp1 , fp1p1));
874  }
875  nearSize = vNbr.size() ;
876 
877  // now on to next-near neighbors
878  if (nm2 < 999) {
879  vNbr.push_back(GetCellID(tm2 , fm200));
880  if (fm2m1 < 999)
881  vNbr.push_back(GetCellID(tm2 , fm2m1));
882  if (fm2p1 < 999)
883  vNbr.push_back(GetCellID(tm2 , fm2p1));
884  if (fm2m2 < 999)
885  vNbr.push_back(GetCellID(tm2 , fm2m2));
886  if (fm2p2 < 999)
887  vNbr.push_back(GetCellID(tm2 , fm2p2));
888  }
889  if (nm1 < 999) {
890  if (fm1m2 < 999)
891  vNbr.push_back(GetCellID(tm1 , fm1m2));
892  if (fm1p2 < 999)
893  vNbr.push_back(GetCellID(tm1 , fm1p2));
894  }
895  vNbr.push_back(GetCellID(t00 , f00m2));
896  vNbr.push_back(GetCellID(t00 , f00p2));
897  if (np1 < 999) {
898  if (fp1m2 < 999)
899  vNbr.push_back(GetCellID(tp1, fp1m2));
900  if (fp1p2 < 999)
901  vNbr.push_back(GetCellID(tp1, fp1p2));
902  }
903  if (np2 < 999) {
904  vNbr.push_back(GetCellID(tp2, fp200));
905  if (fp2m1 < 999)
906  vNbr.push_back(GetCellID(tp2, fp2m1));
907  if (fp2p1 < 999)
908  vNbr.push_back(GetCellID(tp2, fp2p1));
909  if (fp2m2 < 999)
910  vNbr.push_back(GetCellID(tp2, fp2m2));
911  if (fp2p2 < 999)
912  vNbr.push_back(GetCellID(tp2, fp2p2));
913  }
914  }//else( aCellId > 1152 && aCellId < 7777 )
915  return
916  EclNbr(vNbr, nearSize);
917 }
Belle2::ECL::EclNbr::printNbr
void printNbr()
print crystals nbrs
Definition: ECLGeometryPar.cc:537
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
Definition: IntervalOfValidity.cc:196
Belle2::ECL::EclNbr::m_nearSize
std::vector< Identifier >::size_type m_nearSize
size of near brs
Definition: ECLGeometryPar.h:218
Belle2::ECL::EclNbr::operator=
EclNbr & operator=(const EclNbr &aNbr)
assignment operator(s)
Definition: ECLGeometryPar.cc:555
Belle2::ECL::EclNbr::GetCellID
int GetCellID(int ThetaId, int PhiId)
Get Cell Id.
Definition: ECLGeometryPar.cc:619
Belle2::ECL::EclNbr::getNbr
EclNbr getNbr(const Identifier aCellId)
get crystals nbr
Definition: ECLGeometryPar.cc:634
Belle2::ECL::ECLGeometryPar::CrystalGeom_t
crystal geometry
Definition: ECLGeometryPar.h:128
Belle2::ECL::EclNbr::mNbr_thetaID
int mNbr_thetaID
The Theta ID information.
Definition: ECLGeometryPar.h:215
Belle2::ECL::EclNbr::nearBegin
const std::vector< Identifier >::const_iterator nearBegin() const
get crystals nearBegin
Definition: ECLGeometryPar.cc:583
Belle2::ECL::EclNbr::GetThetaID
int GetThetaID()
Get Cell Id.
Definition: ECLGeometryPar.h:198
Mapping_t
Definition: ECLGeometryPar.cc:28
Belle2::ECL::EclNbr
EclNbr class
Definition: ECLGeometryPar.h:149
Belle2::ECL::EclNbr::nbrs
const std::vector< Identifier > & nbrs() const
get crystals nbrs
Definition: ECLGeometryPar.cc:577
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::EclNbr::Mapping
void Mapping(int cid)
Mapping theta, phi Id.
Definition: ECLGeometryPar.cc:627
Belle2::ECL::EclNbr::GetPhiID
int GetPhiID()
Get Theta Id.
Definition: ECLGeometryPar.h:199
Belle2::ECL::EclNbr::nextSize
std::vector< Identifier >::size_type nextSize() const
get crystals nextSize
Definition: ECLGeometryPar.cc:613
Belle2::ECL::EclNbr::nextBegin
const std::vector< Identifier >::const_iterator nextBegin() const
get crystals nextBegin
Definition: ECLGeometryPar.cc:595
Belle2::rn
TString rn()
Get random string.
Definition: tools.h:41
Belle2::ECL::EclNbr::mNbr_phiID
int mNbr_phiID
The Phi ID information.
Definition: ECLGeometryPar.h:216
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition: ECLGeometryPar.h:45
Belle2::ECL::EclNbr::m_nbrs
std::vector< Identifier > & m_nbrs
id of m_brs
Definition: ECLGeometryPar.h:217
Belle2::ECL::EclNbr::mNbr_cellID
int mNbr_cellID
Get Phi Id.
Definition: ECLGeometryPar.h:199
Belle2::ECL::EclNbr::~EclNbr
virtual ~EclNbr()
destructor
Definition: ECLGeometryPar.cc:516
Belle2::ECL::EclNbr::nearSize
std::vector< Identifier >::size_type nearSize() const
get crystals nearSize
Definition: ECLGeometryPar.cc:607
Belle2::ECL::EclNbr::nearEnd
const std::vector< Identifier >::const_iterator nearEnd() const
get crystals nearEnd
Definition: ECLGeometryPar.cc:589
Belle2::ECL::EclNbr::Identifier
EclIdentifier Identifier
constants, enums and typedefs
Definition: ECLGeometryPar.h:154
Belle2::ECL::EclNbr::EclNbr
EclNbr()
Constructors and destructor.
Definition: ECLGeometryPar.cc:487
Belle2::ECL::EclNbr::nextEnd
const std::vector< Identifier >::const_iterator nextEnd() const
get crystals nextEnd
Definition: ECLGeometryPar.cc:601