BmnRoot 26.10.0-210-gc52d34b0
Software framework of the BM@N experiment
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
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 if (!gGeoManager->cd(pathTOF)) {
102 LOG(error) << "[BmnTof1GeoUtils::ParseTGeoManager] TOF400 volume '" << pathTOF
103 << "' was not found in the geometry, please, check the setup geometry";
104 return -1;
105 }
106 Double_t* X0Y0Z0 = new Double_t[3];
107 X0Y0Z0[0] = X0Y0Z0[1] = X0Y0Z0[2] = 0.; // center of sensetive detector
108 Double_t *local = new Double_t[3], master[3], dX, dY, dZ;
109
110 Int_t volumeUID, detectorID, stripID;
111 Int_t nDetectors = 0, nStrips = 0;
112
113 TObjArray* array = gGeoManager->GetCurrentVolume()->GetNodes();
114 TIterator* it1 = array->MakeIterator();
115
116 TGeoNode *detectorNode, *stripNode;
117 while ((detectorNode = (TGeoNode*)it1->Next())) // detectors
118 {
119 TString PATH1 = pathTOF + "/" + detectorNode->GetName();
120 detectorID = detectorNode->GetNumber();
121 TString detectorNodeName = (TString)detectorNode->GetName();
122 if (!detectorNodeName.Contains("Detector"))
123 continue;
124 nDetectors++;
125 // cout<<"\n DETECTOR: "<<detectorNode->GetName()<<", copy# "<<detectorID<<" path= "<<PATH1.Data();
126
127 TIterator* it2 = detectorNode->GetNodes()->MakeIterator();
128 while ((stripNode = (TGeoNode*)it2->Next())) // strips
129 {
130 stripName = stripNode->GetName();
131 if (!stripName.Contains("StripActiveGas"))
132 continue;
133
134 TString PATH2 = PATH1 + "/" + stripName;
135 stripID = stripNode->GetNumber();
136 nStrips++;
137 // cout<<"\n \tSTRIP: "<<stripNode->GetName()<<", copy# "<<stripID<<" path=
138 // "<<PATH2.Data();
139
140 gGeoManager->cd(PATH2);
141
142 TGeoMatrix* matrix = gGeoManager->GetCurrentMatrix(); // calculate global TGeoHMatrix for current node
143 matrix->LocalToMaster(X0Y0Z0, master); // 0.0.0. --> MRS
144
145 TGeoBBox* box = (TGeoBBox*)gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
146 dX = box->GetDX();
147 dY = box->GetDY();
148 dZ = box->GetDZ();
149 // cout<<"\n center: ("<<master[0]<<", "<<master[1]<<", "<<master[2]<<") d("<<dX<<",
150 // "<<dY<<", "<<dZ<<")";
151
152 volumeUID = BmnTOF1Point::GetVolumeUID(0, detectorID, stripID); // regionID == 0 now
153
154 LStrip1 stripData(volumeUID, 0, 0, detectorID, stripID);
155 stripData.center.SetXYZ(master[0], master[1], master[2]);
156
157 // edges on the front plate of the strips. perp along Z.
158 local[0] = -dX;
159 local[1] = -dY;
160 local[2] = -dZ;
161 matrix->LocalToMaster(local, master);
162 stripData.A.SetXYZ(master[0], master[1], master[2]);
163
164 local[0] = +dX;
165 local[1] = -dY;
166 local[2] = -dZ;
167 matrix->LocalToMaster(local, master);
168 stripData.B.SetXYZ(master[0], master[1], master[2]);
169
170 local[0] = +dX;
171 local[1] = +dY;
172 local[2] = -dZ;
173 matrix->LocalToMaster(local, master);
174 stripData.C.SetXYZ(master[0], master[1], master[2]);
175
176 local[0] = -dX;
177 local[1] = +dY;
178 local[2] = -dZ;
179 matrix->LocalToMaster(local, master);
180 stripData.D.SetXYZ(master[0], master[1], master[2]);
181
182 stripData.InitCenterPerp();
183
184 // stripData.Dump("\n strip:");
185 bool IsUniqueUID = mStrips.insert(make_pair(volumeUID, stripData)).second;
186 assert(IsUniqueUID);
187
188 } // strips
189 delete it2;
190 } // detectors
191
192 if (verbose)
193 LOG(info) << "[BmnTof400HitProducer::ParseTGeoManager] detectors = " << nDetectors << ", strips= " << nStrips
194 << ". ";
195 delete it1;
196 delete[] X0Y0Z0;
197 delete[] local;
198 return nDetectors;
199}
200
201Int_t BmnTof1GeoUtils::ParseTGeoManager701(bool useMCinput, bool forced, Int_t verbose)
202{
203 assert(gGeoManager);
204 if (!forced && !mStrips.empty())
205 return -1; // already parsed and filled
206 mStrips.clear();
207 // TString stripName, pathTOF = "/cave_1/TOFB1_0";
208 TString stripName, gasName, pathTOF = "/cave_1/tof700_0";
209 if (!gGeoManager->cd(pathTOF)) {
210 LOG(error) << "[BmnTof1GeoUtils::ParseTGeoManager701] TOF700 volume '" << pathTOF
211 << "' was not found in the geometry, please, check the setup geometry";
212 return -1;
213 }
214 Double_t* X0Y0Z0 = new Double_t[3];
215 X0Y0Z0[0] = X0Y0Z0[1] = X0Y0Z0[2] = 0.; // center of sensetive detector
216 Double_t *local = new Double_t[3], master[3], dX, dY, dZ;
217
218 Int_t volumeUID, detectorID, stripID;
219 Int_t nDetectors = 0, nStrips = 0;
220
221 TObjArray* array = gGeoManager->GetCurrentVolume()->GetNodes();
222 TIterator* it1 = array->MakeIterator();
223
224 TGeoNode *detectorNode, *detectorNodeGas, *stripNode;
225 while ((detectorNode = (TGeoNode*)it1->Next())) // detectors
226 {
227 TString PATH1 = pathTOF + "/" + detectorNode->GetName();
228 detectorID = detectorNode->GetNumber();
229 TString detectorNodeName = detectorNode->GetName(); // aka chamName
230 if (!detectorNodeName.Contains("Chamber"))
231 continue;
232 nDetectors++;
233 // cout<<"\n DETECTOR: "<<detectorNode->GetName()<<", copy# "<<detectorID<<" path= "<<PATH1.Data();
234
235 TIterator* it2 = detectorNode->GetNodes()->MakeIterator();
236 while ((detectorNodeGas = (TGeoNode*)it2->Next())) // gas
237 {
238 gasName = detectorNodeGas->GetName();
239 if (!gasName.Contains("Gas"))
240 continue;
241 TString PATHG = PATH1 + "/" + gasName;
242 TObjArray* Nodes = detectorNodeGas->GetNodes();
243 TIterator* it3 = Nodes->MakeIterator();
244
245 while ((stripNode = (TGeoNode*)it3->Next())) {
246
247 stripName = stripNode->GetName();
248 if (!stripName.Contains("ActiveGasStrip"))
249 continue;
250
251 TString PATH2 = PATHG + "/" + stripName;
252 stripID = stripNode->GetNumber();
253 nStrips++;
254 // cout<<"\n \tSTRIP: "<<stripNode->GetName()<<", copy# "<<stripID<<" path=
255 // "<<PATH2.Data();
256
257 gGeoManager->cd(PATH2);
258
259 TGeoMatrix* matrix = gGeoManager->GetCurrentMatrix(); // calculate global TGeoHMatrix for current node
260 matrix->LocalToMaster(X0Y0Z0, master); // 0.0.0. --> MRS
261
262 TGeoBBox* box = (TGeoBBox*)gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
263 dX = box->GetDX();
264 dY = box->GetDY();
265 dZ = box->GetDZ();
266 // cout<<"\n center: ("<<master[0]<<", "<<master[1]<<", "<<master[2]<<") d("<<dX<<", "<<dY<<",
267 // "<<dZ<<")";
268
269 volumeUID = BmnTOF1Point::GetVolumeUID(0, detectorID, stripID); // regionID == 0 now
270 LStrip1 stripData(volumeUID, 0, 0, detectorID, stripID);
271 if (stripID == 1)
272 stripData.SetIsKilled(true);
273 stripData.center.SetXYZ(master[0], master[1], master[2]);
274
275 // edges on the front plate of the strips. perp along Z.
276 local[0] = -dX;
277 local[1] = -dY;
278 local[2] = -dZ;
279 matrix->LocalToMaster(local, master);
280 stripData.A.SetXYZ(master[0], master[1], master[2]);
281
282 local[0] = +dX;
283 local[1] = -dY;
284 local[2] = -dZ;
285 matrix->LocalToMaster(local, master);
286 stripData.B.SetXYZ(master[0], master[1], master[2]);
287
288 local[0] = +dX;
289 local[1] = +dY;
290 local[2] = -dZ;
291 matrix->LocalToMaster(local, master);
292 stripData.C.SetXYZ(master[0], master[1], master[2]);
293
294 local[0] = -dX;
295 local[1] = +dY;
296 local[2] = -dZ;
297 matrix->LocalToMaster(local, master);
298 stripData.D.SetXYZ(master[0], master[1], master[2]);
299
300 stripData.InitCenterPerp();
301 stripData.CalculateLengthnWidth();
302
303 // stripData.Dump("\n strip:");
304 // cout << stripData.GetLength() << endl << endl;
305 bool IsUniqueUID = mStrips.insert(make_pair(volumeUID, stripData)).second;
306 assert(IsUniqueUID);
307 } // strips
308 delete it3;
309 } // gas
310 delete it2;
311 } // detectors
312
313 if (verbose)
314 LOG(info) << "[BmnTof700HitProducer::ParseTGeoManager] detectors = " << nDetectors << ", strips= " << nStrips
315 << ". ";
316 delete it1;
317 delete[] X0Y0Z0;
318 delete[] local;
319 return nDetectors;
320}
321
323{
324
325 MStripCIT cit = mStrips.find(UID);
326 // assert(cit != mStrips.end());
327 if (cit == mStrips.end())
328 return NULL;
329 return &(cit->second);
330
331 // MStripCIT cit = mStrips.find(UID);
332 // assert(cit != mStrips.end());
333 // return &(cit->second);
334}
335
337{
338 MStripCIT cit = mStrips.find(UID);
339 // assert(cit != mStrips.end());
340 if (cit == mStrips.end())
341 return NULL;
342 return &(cit->second);
343}
344
346 const TVector3& a,
347 const TVector3& b,
348 const TVector3& c,
349 const TVector3& d,
350 bool check)
351 : IsInvalid(false)
352 , volumeUID(uid)
353 , A(a)
354 , B(b)
355 , C(c)
356 , D(d)
357{
358 if (check)
359 CheckInValid();
360}
361
363{
364 assert(P1 != P2);
365
366 return ((*pos - P1).Cross(*pos - P2)).Mag() / (P2 - P1).Mag();
367}
368
370{
371 assert(P1 != P2);
372
373 TVector3 v = P2 - P1;
374 TVector3 w = (*pos) - P1;
375
376 double c1 = w.Dot(v);
377 if (c1 <= 0)
378 return w.Mag();
379
380 double c2 = v.Dot(v);
381 if (c2 <= c1)
382 return ((*pos) - P2).Mag();
383
384 TVector3 Pb = P1 + (c1 / c2) * v;
385 return ((*pos) - Pb).Mag();
386}
387
389{
390 double right = DistanceFromPointToLineSegment(pos, A, D);
391 double left = DistanceFromPointToLineSegment(pos, B, C);
392
393 // sorting & return minimal value
394 if (right <= left) {
395 side = LStrip1::kRight;
396 return right;
397 }
398
399 side = LStrip1::kLeft;
400 return left;
401}
402
404{
405 double down = DistanceFromPointToLineSegment(pos, A, B);
406 double up = DistanceFromPointToLineSegment(pos, C, D);
407
408 // sorting & return minimal value
409 if (down <= up) {
410 side = LStrip1::kLower;
411 return down;
412 }
413
414 side = LStrip1::kUpper;
415 return up;
416}
417
418void LRectangle1::Print(ostream& out, const TVector3& point, const char* comment) const
419{
420 if (comment)
421 out << comment;
422 out << " (" << point.X() << ",\t" << point.Y() << ",\t" << point.Z() << ") ";
423}
424
425void LRectangle1::Dump(const char* comment, ostream& out) const
426{
427 if (comment) {
428 out << comment;
429 out << " uid=" << volumeUID << " IsInvalid=" << IsInvalid;
430 }
431 Print(out, A, "\tA:");
432 Print(out, B, "\tB:");
433 Print(out, C, "\tC:");
434 Print(out, D, "\tD:");
435}
436
446
448 : sectorID(sector)
449 , boxID(box)
451 , stripID(strip)
452{
453 volumeUID = uid;
456}
457
458void LStrip1::Dump(const char* comment, ostream& out) const
459{
460 if (comment)
461 out << comment;
462 out << " ids: " << sectorID << ", " << boxID << ", " << detectorID << ", " << stripID << ", " << GetIsKilled();
463
464 LRectangle1::Dump(nullptr, out);
465}
466
468{
469 Double_t value, min1 = 1.e+10, min2 = 1.e+10; // big value
470
471 if ((*this) == strip)
472 return min1 + min2; // same strip
473 if (!IsSameDetector(strip))
474 return min1 + min2; // different modules
475
476 TVector3 *p1, *p2;
477 switch (side) {
478 case kRight:
479 p1 = &A;
480 p2 = &D;
481 break;
482 case kLeft:
483 p1 = &B;
484 p2 = &C;
485 break;
486 case kInvalid:
487 break;
488 }
489
490 value = fabs((*p1 - strip.A).Mag());
491 min1 = (value < min1) ? value : min1;
492 value = fabs((*p1 - strip.B).Mag());
493 min1 = (value < min1) ? value : min1;
494 value = fabs((*p1 - strip.C).Mag());
495 min1 = (value < min1) ? value : min1;
496 value = fabs((*p1 - strip.D).Mag());
497 min1 = (value < min1) ? value : min1;
498
499 value = fabs((*p2 - strip.A).Mag());
500 min2 = (value < min2) ? value : min2;
501 value = fabs((*p2 - strip.B).Mag());
502 min2 = (value < min2) ? value : min2;
503 value = fabs((*p2 - strip.C).Mag());
504 min2 = (value < min2) ? value : min2;
505 value = fabs((*p2 - strip.D).Mag());
506 min2 = (value < min2) ? value : min2;
507
508 return min1 + min2;
509}
510
512{
513 Double_t value, min1 = 1.e+10, min2 = 1.e+10; // big value
514
515 if ((*this) == strip)
516 return min1 + min2; // same strip
517 if (!IsSameDetector(strip))
518 return min1 + min2; // different modules
519
520 TVector3 *p1, *p2;
521 switch (side) {
522 case kUpper:
523 p1 = &C;
524 p2 = &D;
525 break; // ????
526 case kLower:
527 p1 = &B;
528 p2 = &A;
529 break;
530 case kInvalid:
531 break;
532 };
533
534 value = fabs((*p1 - strip.A).Mag());
535 min1 = (value < min1) ? value : min1;
536 value = fabs((*p1 - strip.B).Mag());
537 min1 = (value < min1) ? value : min1;
538 value = fabs((*p1 - strip.C).Mag());
539 min1 = (value < min1) ? value : min1;
540 value = fabs((*p1 - strip.D).Mag());
541 min1 = (value < min1) ? value : min1;
542
543 value = fabs((*p2 - strip.A).Mag());
544 min2 = (value < min2) ? value : min2;
545 value = fabs((*p2 - strip.B).Mag());
546 min2 = (value < min2) ? value : min2;
547 value = fabs((*p2 - strip.C).Mag());
548 min2 = (value < min2) ? value : min2;
549 value = fabs((*p2 - strip.D).Mag());
550 min2 = (value < min2) ? value : min2;
551 // cout << "min1 = " << min1 << "; min2 = " << min2 << endl;
552 return min1 + min2;
553}
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
__m128 v
Definition P4_F32vec4.h:1
int Int_t
double Double_t
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
R__EXTERN TGeoManager * gGeoManager
Int_t GetVolumeUID() const
const LStrip1 * FindStrip(Int_t UID)
void FindNeighborStrips701(Int_t verbose=0)
MStripType::iterator MStripIT
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)
MStripType::const_iterator MStripCIT
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 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
void SetIsKilled(bool status)
Int_t neighboring[2]
bool GetIsKilled() const
Double_t Distance701(Side_t side, const LStrip1 &strip)
Int_t detectorID
Double_t Distance(Side_t side, const LStrip1 &strip)
virtual Double_t GetDX() const
virtual Double_t GetDZ() const
virtual Double_t GetDY() const
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
Int_t GetNumber() const
TObjArray * GetNodes() const
virtual TObject * Next()=0
const char * GetName() const override
TIterator * MakeIterator(Bool_t dir=kIterForward) const override
virtual void Dump() const
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Double_t Z() const
void SetXYZ(Double_t x, Double_t y, Double_t z)
Double_t Y() const
Double_t Dot(const TVector3 &) const
Double_t X() const
Double_t Mag() const