BmnRoot
Loading...
Searching...
No Matches
BmnTof1GeoUtils.cxx
Go to the documentation of this file.
1#include "BmnTof1GeoUtils.h"
2
3#include "BmnTOF1Point.h"
4#include "BmnTofHit.h"
5#include "FairLogger.h"
6#include "FairRootManager.h"
7#include "TGeoBBox.h"
8#include "TGeoMatrix.h"
9#include "TGeoNode.h"
10
11#include <TGeoManager.h>
12#include <assert.h>
13using namespace std;
14
16
18{
19 size_t NR = 0, NL = 0;
20 const LStrip1* strip2;
21 double distance;
22 for (MStripIT it1 = mStrips.begin(), itEnd1 = mStrips.end(); it1 != itEnd1; it1++) // cycle1 by strips
23 {
24 LStrip1* strip1 = &(it1->second);
25
26 for (MStripCIT it2 = mStrips.begin(), itEnd2 = mStrips.end(); it2 != itEnd2; it2++) // cycle2 by strips
27 {
28 strip2 = &(it2->second);
29
30 // CATION: Ckeck only left and right sides(one row strips NOW)
31 distance = strip1->Distance(LStrip1::kRight, *strip2);
32 if (distance < 0.8) // CAUTION: constant depends on the geometry layout
33 {
34 strip1->neighboring[LStrip1::kRight] = strip2->volumeUID;
35 NR++;
36 }
37
38 distance = strip1->Distance(LStrip1::kLeft, *strip2);
39 if (distance < 0.8) // CAUTION: constant depends on the geometry layout
40 {
41 strip1->neighboring[LStrip1::kLeft] = strip2->volumeUID;
42 NL++;
43 }
44
45 } // cycle2 by strips
46 } // cycle1 by strips
47
48 if (verbose)
49 cout << " [BmnTof1GeoUtils::FindNeighborStrips] Neighbor strips: left = " << NL << ", right = " << NR << endl;
50}
51
53{
54 size_t LO = 0, UP = 0;
55 const LStrip1* strip2;
56 double distance;
57 for (MStripIT it1 = mStrips.begin(), itEnd1 = mStrips.end(); it1 != itEnd1; it1++) // cycle1 by strips
58 {
59 LStrip1* strip1 = &(it1->second);
60
61 for (MStripCIT it2 = mStrips.begin(), itEnd2 = mStrips.end(); it2 != itEnd2; it2++) // cycle2 by strips
62 {
63 strip2 = &(it2->second);
64
65 // CATION: Ckeck only left and right sides(one row strips NOW)
66 distance = strip1->Distance701(LStrip1::kLower, *strip2);
67
68 if (distance < 0.8 && !strip2->GetIsKilled()) // CAUTION: constant depends on the geometry layout
69 {
70 strip1->neighboring[LStrip1::kLower] = strip2->volumeUID;
71 LO++;
72 } else if (distance < 0.8 && strip2->GetIsKilled()) {
73
74 strip1->neighboring[LStrip1::kLower] = LRectangle1::Side_t::kInvalid;
75 }
76 distance = strip1->Distance701(LStrip1::kUpper, *strip2);
77
78 if (distance < 0.8) // CAUTION: constant depends on the geometry layout
79 {
80 strip1->neighboring[LStrip1::kUpper] = strip2->volumeUID;
81 UP++;
82 }
83
84 } // cycle2 by strips
85 } // cycle1 by strips
86 // const LStrip1* pstrip;
87
88 if (verbose)
89 cout << " [BmnTof1GeoUtils::FindNeighborStrips701] Neighbor strips: lower = " << LO << ", upper = " << UP
90 << endl;
91}
92
93Int_t BmnTof1GeoUtils::ParseTGeoManager(bool useMCinput, bool forced, Int_t verbose)
94{
95 assert(gGeoManager);
96 if (!forced && !mStrips.empty())
97 return -1; // already parsed and filled
98 mStrips.clear();
99 // TString stripName, pathTOF = "/cave_1/TOFB1_0";
100 TString stripName, pathTOF = "/cave_1/TOF400_0";
101 gGeoManager->cd(pathTOF);
102 Double_t* X0Y0Z0 = new Double_t[3];
103 X0Y0Z0[0] = X0Y0Z0[1] = X0Y0Z0[2] = 0.; // center of sensetive detector
104 Double_t *local = new Double_t[3], master[3], dX, dY, dZ;
105
106 Int_t volumeUID, detectorID, stripID;
107 Int_t nDetectors = 0, nStrips = 0;
108
109 TObjArray* array = gGeoManager->GetCurrentVolume()->GetNodes();
110 TIterator* it1 = array->MakeIterator();
111
112 TGeoNode *detectorNode, *stripNode;
113 while ((detectorNode = (TGeoNode*)it1->Next())) // detectors
114 {
115 TString PATH1 = pathTOF + "/" + detectorNode->GetName();
116 detectorID = detectorNode->GetNumber();
117 TString detectorNodeName = (TString)detectorNode->GetName();
118 if (!detectorNodeName.Contains("Detector"))
119 continue;
120 nDetectors++;
121 // cout<<"\n DETECTOR: "<<detectorNode->GetName()<<", copy# "<<detectorID<<" path= "<<PATH1.Data();
122
123 TIterator* it2 = detectorNode->GetNodes()->MakeIterator();
124 while ((stripNode = (TGeoNode*)it2->Next())) // strips
125 {
126 stripName = stripNode->GetName();
127 if (!stripName.Contains("StripActiveGas"))
128 continue;
129
130 TString PATH2 = PATH1 + "/" + stripName;
131 stripID = stripNode->GetNumber();
132 nStrips++;
133 // cout<<"\n \tSTRIP: "<<stripNode->GetName()<<", copy# "<<stripID<<" path=
134 // "<<PATH2.Data();
135
136 gGeoManager->cd(PATH2);
137
138 TGeoMatrix* matrix = gGeoManager->GetCurrentMatrix(); // calculate global TGeoHMatrix for current node
139 matrix->LocalToMaster(X0Y0Z0, master); // 0.0.0. --> MRS
140
141 TGeoBBox* box = (TGeoBBox*)gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
142 dX = box->GetDX();
143 dY = box->GetDY();
144 dZ = box->GetDZ();
145 // cout<<"\n center: ("<<master[0]<<", "<<master[1]<<", "<<master[2]<<") d("<<dX<<",
146 // "<<dY<<", "<<dZ<<")";
147
148 volumeUID = BmnTOF1Point::GetVolumeUID(0, detectorID, stripID); // regionID == 0 now
149
150 LStrip1 stripData(volumeUID, 0, 0, detectorID, stripID);
151 stripData.center.SetXYZ(master[0], master[1], master[2]);
152
153 // edges on the front plate of the strips. perp along Z.
154 local[0] = -dX;
155 local[1] = -dY;
156 local[2] = -dZ;
157 matrix->LocalToMaster(local, master);
158 stripData.A.SetXYZ(master[0], master[1], master[2]);
159
160 local[0] = +dX;
161 local[1] = -dY;
162 local[2] = -dZ;
163 matrix->LocalToMaster(local, master);
164 stripData.B.SetXYZ(master[0], master[1], master[2]);
165
166 local[0] = +dX;
167 local[1] = +dY;
168 local[2] = -dZ;
169 matrix->LocalToMaster(local, master);
170 stripData.C.SetXYZ(master[0], master[1], master[2]);
171
172 local[0] = -dX;
173 local[1] = +dY;
174 local[2] = -dZ;
175 matrix->LocalToMaster(local, master);
176 stripData.D.SetXYZ(master[0], master[1], master[2]);
177
178 stripData.InitCenterPerp();
179
180 // stripData.Dump("\n strip:");
181 bool IsUniqueUID = mStrips.insert(make_pair(volumeUID, stripData)).second;
182 assert(IsUniqueUID);
183
184 } // strips
185 delete it2;
186 } // detectors
187
188 if (verbose)
189 LOG(info) << "[BmnTof400HitProducer::ParseTGeoManager] detectors = " << nDetectors << ", strips= " << nStrips
190 << ". ";
191 delete it1;
192 delete[] X0Y0Z0;
193 delete[] local;
194 return nDetectors;
195}
196
197Int_t BmnTof1GeoUtils::ParseTGeoManager701(bool useMCinput, bool forced, Int_t verbose)
198{
199 assert(gGeoManager);
200 if (!forced && !mStrips.empty())
201 return -1; // already parsed and filled
202 mStrips.clear();
203 // TString stripName, pathTOF = "/cave_1/TOFB1_0";
204 TString stripName, gasName, pathTOF = "/cave_1/tof700_0";
205 gGeoManager->cd(pathTOF);
206 Double_t* X0Y0Z0 = new Double_t[3];
207 X0Y0Z0[0] = X0Y0Z0[1] = X0Y0Z0[2] = 0.; // center of sensetive detector
208 Double_t *local = new Double_t[3], master[3], dX, dY, dZ;
209
210 Int_t volumeUID, detectorID, stripID;
211 Int_t nDetectors = 0, nStrips = 0;
212
213 TObjArray* array = gGeoManager->GetCurrentVolume()->GetNodes();
214 TIterator* it1 = array->MakeIterator();
215
216 TGeoNode *detectorNode, *detectorNodeGas, *stripNode;
217 while ((detectorNode = (TGeoNode*)it1->Next())) // detectors
218 {
219 TString PATH1 = pathTOF + "/" + detectorNode->GetName();
220 detectorID = detectorNode->GetNumber();
221 TString detectorNodeName = detectorNode->GetName(); // aka chamName
222 if (!detectorNodeName.Contains("Chamber"))
223 continue;
224 nDetectors++;
225 // cout<<"\n DETECTOR: "<<detectorNode->GetName()<<", copy# "<<detectorID<<" path= "<<PATH1.Data();
226
227 TIterator* it2 = detectorNode->GetNodes()->MakeIterator();
228 while ((detectorNodeGas = (TGeoNode*)it2->Next())) // gas
229 {
230 gasName = detectorNodeGas->GetName();
231 if (!gasName.Contains("Gas"))
232 continue;
233 TString PATHG = PATH1 + "/" + gasName;
234 TObjArray* Nodes = detectorNodeGas->GetNodes();
235 TIterator* it3 = Nodes->MakeIterator();
236
237 while ((stripNode = (TGeoNode*)it3->Next())) {
238
239 stripName = stripNode->GetName();
240 if (!stripName.Contains("ActiveGasStrip"))
241 continue;
242
243 TString PATH2 = PATHG + "/" + stripName;
244 stripID = stripNode->GetNumber();
245 nStrips++;
246 // cout<<"\n \tSTRIP: "<<stripNode->GetName()<<", copy# "<<stripID<<" path=
247 // "<<PATH2.Data();
248
249 gGeoManager->cd(PATH2);
250
251 TGeoMatrix* matrix = gGeoManager->GetCurrentMatrix(); // calculate global TGeoHMatrix for current node
252 matrix->LocalToMaster(X0Y0Z0, master); // 0.0.0. --> MRS
253
254 TGeoBBox* box = (TGeoBBox*)gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
255 dX = box->GetDX();
256 dY = box->GetDY();
257 dZ = box->GetDZ();
258 // cout<<"\n center: ("<<master[0]<<", "<<master[1]<<", "<<master[2]<<") d("<<dX<<", "<<dY<<",
259 // "<<dZ<<")";
260
261 volumeUID = BmnTOF1Point::GetVolumeUID(0, detectorID, stripID); // regionID == 0 now
262 LStrip1 stripData(volumeUID, 0, 0, detectorID, stripID);
263 if (stripID == 1)
264 stripData.SetIsKilled(true);
265 stripData.center.SetXYZ(master[0], master[1], master[2]);
266
267 // edges on the front plate of the strips. perp along Z.
268 local[0] = -dX;
269 local[1] = -dY;
270 local[2] = -dZ;
271 matrix->LocalToMaster(local, master);
272 stripData.A.SetXYZ(master[0], master[1], master[2]);
273
274 local[0] = +dX;
275 local[1] = -dY;
276 local[2] = -dZ;
277 matrix->LocalToMaster(local, master);
278 stripData.B.SetXYZ(master[0], master[1], master[2]);
279
280 local[0] = +dX;
281 local[1] = +dY;
282 local[2] = -dZ;
283 matrix->LocalToMaster(local, master);
284 stripData.C.SetXYZ(master[0], master[1], master[2]);
285
286 local[0] = -dX;
287 local[1] = +dY;
288 local[2] = -dZ;
289 matrix->LocalToMaster(local, master);
290 stripData.D.SetXYZ(master[0], master[1], master[2]);
291
292 stripData.InitCenterPerp();
293 stripData.CalculateLengthnWidth();
294
295 // stripData.Dump("\n strip:");
296 // cout << stripData.GetLength() << endl << endl;
297 bool IsUniqueUID = mStrips.insert(make_pair(volumeUID, stripData)).second;
298 assert(IsUniqueUID);
299 } // strips
300 delete it3;
301 } // gas
302 delete it2;
303 } // detectors
304
305 if (verbose)
306 LOG(info) << "[BmnTof700HitProducer::ParseTGeoManager] detectors = " << nDetectors << ", strips= " << nStrips
307 << ". ";
308 delete it1;
309 delete[] X0Y0Z0;
310 delete[] local;
311 return nDetectors;
312}
313
315{
316
317 MStripCIT cit = mStrips.find(UID);
318 // assert(cit != mStrips.end());
319 if (cit == mStrips.end())
320 return NULL;
321 return &(cit->second);
322
323 // MStripCIT cit = mStrips.find(UID);
324 // assert(cit != mStrips.end());
325 // return &(cit->second);
326}
327
329{
330 MStripCIT cit = mStrips.find(UID);
331 // assert(cit != mStrips.end());
332 if (cit == mStrips.end())
333 return NULL;
334 return &(cit->second);
335}
336
338 const TVector3& a,
339 const TVector3& b,
340 const TVector3& c,
341 const TVector3& d,
342 bool check)
343 : IsInvalid(false)
344 , volumeUID(uid)
345 , A(a)
346 , B(b)
347 , C(c)
348 , D(d)
349{
350 if (check)
351 CheckInValid();
352}
353
354Double_t LRectangle1::DistanceFromPointToLine(const TVector3* pos, const TVector3& P1, const TVector3& P2) const
355{
356 assert(P1 != P2);
357
358 return ((*pos - P1).Cross(*pos - P2)).Mag() / (P2 - P1).Mag();
359}
360
361Double_t LRectangle1::DistanceFromPointToLineSegment(const TVector3* pos, const TVector3& P1, const TVector3& P2) const
362{
363 assert(P1 != P2);
364
365 TVector3 v = P2 - P1;
366 TVector3 w = (*pos) - P1;
367
368 double c1 = w.Dot(v);
369 if (c1 <= 0)
370 return w.Mag();
371
372 double c2 = v.Dot(v);
373 if (c2 <= c1)
374 return ((*pos) - P2).Mag();
375
376 TVector3 Pb = P1 + (c1 / c2) * v;
377 return ((*pos) - Pb).Mag();
378}
379
380Double_t LRectangle1::MinDistanceToEdge(const TVector3* pos, Side_t& side) const
381{
382 double right = DistanceFromPointToLineSegment(pos, A, D);
383 double left = DistanceFromPointToLineSegment(pos, B, C);
384
385 // sorting & return minimal value
386 if (right <= left) {
387 side = LStrip1::kRight;
388 return right;
389 }
390
391 side = LStrip1::kLeft;
392 return left;
393}
394
395Double_t LRectangle1::MinDistanceToEdge701(const TVector3* pos, Side_t& side) const
396{
397 double down = DistanceFromPointToLineSegment(pos, A, B);
398 double up = DistanceFromPointToLineSegment(pos, C, D);
399
400 // sorting & return minimal value
401 if (down <= up) {
402 side = LStrip1::kLower;
403 return down;
404 }
405
406 side = LStrip1::kUpper;
407 return up;
408}
409
410void LRectangle1::Print(ostream& out, const TVector3& point, const char* comment) const
411{
412 if (comment)
413 out << comment;
414 out << " (" << point.X() << ",\t" << point.Y() << ",\t" << point.Z() << ") ";
415}
416
417void LRectangle1::Dump(const char* comment, ostream& out) const
418{
419 if (comment) {
420 out << comment;
421 out << " uid=" << volumeUID << " IsInvalid=" << IsInvalid;
422 }
423 Print(out, A, "\tA:");
424 Print(out, B, "\tB:");
425 Print(out, C, "\tC:");
426 Print(out, D, "\tD:");
427}
428
430 : sectorID(kInvalid)
431 , boxID(kInvalid)
432 , detectorID(kInvalid)
433 , stripID(kInvalid)
434{
437}
438
439LStrip1::LStrip1(Int_t uid, Int_t sector, Int_t box, Int_t detector, Int_t strip)
440 : sectorID(sector)
441 , boxID(box)
442 , detectorID(detector)
443 , stripID(strip)
444{
445 volumeUID = uid;
448}
449
450void LStrip1::Dump(const char* comment, ostream& out) const
451{
452 if (comment)
453 out << comment;
454 out << " ids: " << sectorID << ", " << boxID << ", " << detectorID << ", " << stripID << ", " << GetIsKilled();
455
456 LRectangle1::Dump(nullptr, out);
457}
458
459Double_t LStrip1::Distance(Side_t side, const LStrip1& strip)
460{
461 Double_t value, min1 = 1.e+10, min2 = 1.e+10; // big value
462
463 if ((*this) == strip)
464 return min1 + min2; // same strip
465 if (!IsSameDetector(strip))
466 return min1 + min2; // different modules
467
468 TVector3 *p1, *p2;
469 switch (side) {
470 case kRight:
471 p1 = &A;
472 p2 = &D;
473 break;
474 case kLeft:
475 p1 = &B;
476 p2 = &C;
477 break;
478 case kInvalid:
479 break;
480 }
481
482 value = fabs((*p1 - strip.A).Mag());
483 min1 = (value < min1) ? value : min1;
484 value = fabs((*p1 - strip.B).Mag());
485 min1 = (value < min1) ? value : min1;
486 value = fabs((*p1 - strip.C).Mag());
487 min1 = (value < min1) ? value : min1;
488 value = fabs((*p1 - strip.D).Mag());
489 min1 = (value < min1) ? value : min1;
490
491 value = fabs((*p2 - strip.A).Mag());
492 min2 = (value < min2) ? value : min2;
493 value = fabs((*p2 - strip.B).Mag());
494 min2 = (value < min2) ? value : min2;
495 value = fabs((*p2 - strip.C).Mag());
496 min2 = (value < min2) ? value : min2;
497 value = fabs((*p2 - strip.D).Mag());
498 min2 = (value < min2) ? value : min2;
499
500 return min1 + min2;
501}
502
503Double_t LStrip1::Distance701(Side_t side, const LStrip1& strip)
504{
505 Double_t value, min1 = 1.e+10, min2 = 1.e+10; // big value
506
507 if ((*this) == strip)
508 return min1 + min2; // same strip
509 if (!IsSameDetector(strip))
510 return min1 + min2; // different modules
511
512 TVector3 *p1, *p2;
513 switch (side) {
514 case kUpper:
515 p1 = &C;
516 p2 = &D;
517 break; // ????
518 case kLower:
519 p1 = &B;
520 p2 = &A;
521 break;
522 case kInvalid:
523 break;
524 };
525
526 value = fabs((*p1 - strip.A).Mag());
527 min1 = (value < min1) ? value : min1;
528 value = fabs((*p1 - strip.B).Mag());
529 min1 = (value < min1) ? value : min1;
530 value = fabs((*p1 - strip.C).Mag());
531 min1 = (value < min1) ? value : min1;
532 value = fabs((*p1 - strip.D).Mag());
533 min1 = (value < min1) ? value : min1;
534
535 value = fabs((*p2 - strip.A).Mag());
536 min2 = (value < min2) ? value : min2;
537 value = fabs((*p2 - strip.B).Mag());
538 min2 = (value < min2) ? value : min2;
539 value = fabs((*p2 - strip.C).Mag());
540 min2 = (value < min2) ? value : min2;
541 value = fabs((*p2 - strip.D).Mag());
542 min2 = (value < min2) ? value : min2;
543 // cout << "min1 = " << min1 << "; min2 = " << min2 << endl;
544 return min1 + min2;
545}
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
__m128 v
Definition P4_F32vec4.h:1
Int_t GetVolumeUID() const
const LStrip1 * FindStrip(Int_t UID)
void FindNeighborStrips701(Int_t verbose=0)
const LStrip1 * FindStrip701(Int_t UID)
Int_t ParseTGeoManager701(bool useMCinput, bool forced=false, Int_t verbose=0)
BmnTof1GeoUtils()
indexing strips by detectorUID
Int_t ParseTGeoManager(bool useMCinput, bool forced=false, Int_t verbose=0)
void FindNeighborStrips(Int_t verbose=0)
void CheckInValid()
Double_t MinDistanceToEdge701(const TVector3 *pos, Side_t &side) const
void InitCenterPerp()
Double_t DistanceFromPointToLineSegment(const TVector3 *pos, const TVector3 &P1, const TVector3 &P2) const
Double_t DistanceFromPointToLine(const TVector3 *pos, const TVector3 &P1, const TVector3 &P2) const
void Dump(const char *comment, std::ostream &out=std::cout) const
void Print(std::ostream &out, const TVector3 &point, const char *comment=nullptr) const
void CalculateLengthnWidth()
TVector3 center
Double_t MinDistanceToEdge(const TVector3 *pos, Side_t &side) const
bool IsSameDetector(const LStrip1 &strip) const
Int_t neighboring[2]
void SetIsKilled(bool status)
bool GetIsKilled() const
Double_t Distance701(Side_t side, const LStrip1 &strip)
Int_t detectorID
Double_t Distance(Side_t side, const LStrip1 &strip)
void Dump(const char *comment=nullptr, std::ostream &out=std::cout) const
STL namespace.