BmnRoot
Loading...
Searching...
No Matches
BmnSiBTModule.cxx
Go to the documentation of this file.
1#include "BmnSiBTModule.h"
2
3#include "TF1.h"
4#include "TH1F.h"
5#include "TMath.h"
6#include "TRandom.h"
7
8#include <fstream>
9
10BmnSiBTModule::BmnSiBTModule()
11{
12
13 Verbosity = true;
14
15 ZStartModulePosition = 0.0;
16
17 ModuleThickness = 0.03; // cm
18
19 XMinModule = 0.0;
20 XMaxModule = 0.0;
21 YMinModule = 0.0;
22 YMaxModule = 0.0;
23}
24
25BmnSiBTModule::BmnSiBTModule(Double_t z_start_pos, Double_t module_thick)
26{
27
28 Verbosity = true;
29
30 ZStartModulePosition = z_start_pos;
31
32 ModuleThickness = module_thick; // cm
33
34 XMinModule = 0.0;
35 XMaxModule = 0.0;
36 YMinModule = 0.0;
37 YMaxModule = 0.0;
38}
39
40BmnSiBTModule::~BmnSiBTModule() {}
41
42Double_t BmnSiBTModule::GetZPositionRegistered()
43{
44
45 return ZStartModulePosition + ModuleThickness * 0.5;
46}
47
48void BmnSiBTModule::AddStripLayer(BmnSiBTLayer strip_layer)
49{
50 StripLayers.push_back(strip_layer);
51
52 DefineModuleBorders();
53}
54
55Bool_t BmnSiBTModule::SetStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
56{
57 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
58 return StripLayers[layer_num].SetStripSignal(strip_num, signal);
59
60 return false;
61}
62
63Bool_t BmnSiBTModule::AddStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
64{
65 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
66 return StripLayers[layer_num].AddStripSignal(strip_num, signal);
67
68 return false;
69}
70
71Bool_t BmnSiBTModule::SetStripMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch strip_match)
72{
73 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
74 return StripLayers[layer_num].SetStripMatch(strip_num, strip_match);
75
76 return false;
77}
78
79Double_t BmnSiBTModule::GetStripSignalInLayer(Int_t layer_num, Int_t strip_num)
80{
81 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
82 return StripLayers[layer_num].GetStripSignal(strip_num);
83
84 return -1;
85}
86
87BmnMatch BmnSiBTModule::GetStripMatchInLayer(Int_t layer_num, Int_t strip_num)
88{
89 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
90 return StripLayers[layer_num].GetStripMatch(strip_num);
91
92 return BmnMatch();
93}
94
95void BmnSiBTModule::ResetModuleData()
96{
97 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
98 StripLayers[ilayer].ResetLayer();
99 }
100
101 DefineModuleBorders();
102
103 ResetRealPoints();
104 ResetIntersectionPoints();
105}
106
107Bool_t BmnSiBTModule::IsPointInsideModule(Double_t x, Double_t y)
108{
109 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
110 if (StripLayers[ilayer].IsPointInsideStripLayer(x, y))
111 return true;
112 }
113 return false;
114}
115
116Bool_t BmnSiBTModule::IsPointInsideModule(Double_t x, Double_t y, Double_t z)
117{
118 Double_t coord_eps = 0.01; // 100 um
119 if (((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps))
120 && ((z < (ZStartModulePosition + ModuleThickness))
121 || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps))
122 {
123 if (IsPointInsideModule(x, y))
124 return true;
125 }
126 return false;
127}
128
129Bool_t BmnSiBTModule::IsPointInsideZThickness(Double_t z)
130{
131 Double_t coord_eps = 0.01; // 100 um
132 if (((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps))
133 && ((z < (ZStartModulePosition + ModuleThickness))
134 || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps))
135 {
136 return true;
137 }
138 return false;
139}
140
141// Add single point with Gaussian smearing
142Bool_t BmnSiBTModule::AddRealPointFullOne(Double_t x,
143 Double_t y,
144 Double_t z,
145 Double_t px,
146 Double_t py,
147 Double_t pz,
148 Double_t signal,
149 Int_t refID)
150{
151
152 if (IsPointInsideModule(x, y)) {
153
154 if (signal <= 0.0)
155 signal = 1e-16;
156
157 Double_t AvalancheRadius = 0.05; // cm (avalanche is about 2-3 strips)
158
159 Double_t radius = AvalancheRadius;
160 if (radius <= 0.0)
161 return false;
162
163 Int_t cycle_cnt = 0;
164 while (true) {
165 radius = gRandom->Gaus(AvalancheRadius, 0.025);
166 if (radius > AvalancheRadius / 2.0 && radius < AvalancheRadius * 2.0 && radius > 0.0)
167 break;
168 cycle_cnt++;
169
170 if (cycle_cnt > 5) {
171 radius = AvalancheRadius;
172 break;
173 }
174 }
175
176 // Make strip cluster in each strip layer -------------------------------
177 Int_t n_strip_layers = StripLayers.size();
178 vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
179
180 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
181 if (!StripLayers[ilayer].IsPointInsideStripLayer(x, y))
182 continue;
183 cluster_layers[ilayer] = MakeCluster(ilayer, x, y, signal, radius);
184 }
185 //----------------------------------------------------------------------
186
187 // Add the correct clusters to the strip layers -------------------------
188 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
189 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
190 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
191 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
192 StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
193 }
194 }
195 //----------------------------------------------------------------------
196
197 // Fill strip matches ---------------------------------------------------
198 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
199 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
200 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
201 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
202 StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
203 refID);
204 }
205 }
206 //----------------------------------------------------------------------
207
208 RealPointsX.push_back(x);
209 RealPointsY.push_back(y);
210 RealPointsMC.push_back(refID);
211
212 return true;
213 } else {
214 if (Verbosity)
215 cerr << "WARNING: BmnSiBTModule::AddRealPointFullOne(): Point (" << x << " : " << y << " : " << z
216 << ") is out of the readout plane or inside a dead zone\n";
217 return false;
218 }
219}
220
221// Add single point without smearing and avalanch effects
222Bool_t BmnSiBTModule::AddRealPointSimple(Double_t x,
223 Double_t y,
224 Double_t z,
225 Double_t px,
226 Double_t py,
227 Double_t pz,
228 Double_t signal,
229 Int_t refID)
230{
231
232 if (IsPointInsideModule(x, y)) {
233
234 Int_t n_strip_layers = StripLayers.size();
235
236 // Add a signal to the strips layers ------------------------------------
237 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
238
239 if (!StripLayers[ilayer].IsPointInsideStripLayer(x, y))
240 continue;
241
242 Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
243 StripLayers[ilayer].AddStripSignal((Int_t)strip_pos, signal);
244
245 // strip match
246 StripLayers[ilayer].AddLinkToStripMatch((Int_t)strip_pos, 1.0, refID);
247 }
248 //----------------------------------------------------------------------
249
250 RealPointsX.push_back(x);
251 RealPointsY.push_back(y);
252 RealPointsMC.push_back(refID);
253
254 return true;
255 } else {
256 if (Verbosity)
257 cerr << "WARNING: BmnSiBTModule::AddRealPointSimple(): Point (" << x << " : " << y
258 << ") is out of the readout plane or inside a dead zone\n";
259 return false;
260 }
261}
262
263StripCluster BmnSiBTModule::MakeCluster(Int_t layer_num,
264 Double_t xcoord,
265 Double_t ycoord,
266 Double_t signal,
267 Double_t radius)
268{
269
270 Double_t ClusterDistortion = 0.0; // signal variation (0.1 is 10%)
271
272 Double_t Gain = 1.0;
273
274 StripCluster cluster;
275
276 if (radius <= 0.0)
277 return cluster;
278
279 Double_t Pitch = StripLayers[layer_num].GetPitch();
280
281 Double_t RadiusInZones = radius / Pitch;
282 Double_t Sigma = RadiusInZones / 3.33;
283
284 TF1 gausF("gausF", "gaus", -4 * Sigma, 4 * Sigma);
285 gausF.SetParameter(0, 1.0); // constant (altitude)
286 gausF.SetParameter(1, 0.0); // mean (center position)
287 gausF.SetParameter(2, Sigma); // sigma
288
289 Double_t SCluster = gausF.Integral(-4 * Sigma, 4 * Sigma); // square of the one side distribution (more exactly)
290
291 TRandom rand(0);
292 Double_t var_level = ClusterDistortion; // signal variation (0.1 is 10%)
293
294 Int_t NStripsInLayer = StripLayers[layer_num].GetNStrips();
295
296 if (NStripsInLayer == 0)
297 return cluster;
298
299 Double_t CenterZonePos = StripLayers[layer_num].ConvertPointToStripPosition(xcoord, ycoord);
300
301 if (CenterZonePos < 0.0 || CenterZonePos >= NStripsInLayer)
302 return cluster;
303
304 gausF.SetParameter(1, CenterZonePos);
305 Double_t total_signal = 0.0;
306
307 Double_t LeftZonePos = CenterZonePos - RadiusInZones;
308 Double_t RightZonePos = CenterZonePos + RadiusInZones;
309 Double_t OutLeftBorder = 0.0;
310 Double_t OutRightBorder = 0.0;
311
312 if (LeftZonePos < 0.0) {
313 OutLeftBorder = LeftZonePos;
314 LeftZonePos = 0.0;
315 }
316 if (RightZonePos < 0.0) {
317 OutRightBorder = RightZonePos;
318 RightZonePos = 0.0;
319 }
320 if (LeftZonePos >= NStripsInLayer) {
321 OutLeftBorder = LeftZonePos;
322 LeftZonePos = NStripsInLayer - 0.001;
323 }
324 if (RightZonePos >= NStripsInLayer) {
325 OutRightBorder = RightZonePos;
326 RightZonePos = NStripsInLayer - 0.001;
327 }
328
329 Double_t h = 0.0;
330 Double_t dist = 0.0;
331
332 // avalanche is inside of the one strip
333 if ((Int_t)LeftZonePos == (Int_t)RightZonePos) {
334
335 Int_t NumCurrentZone = (Int_t)LeftZonePos;
336
337 h = RightZonePos - LeftZonePos;
338 if (h < 0.0)
339 h = 0.0;
340
341 Double_t xp = LeftZonePos + dist;
342 Double_t S = gausF.Integral(xp, xp + h);
343 Double_t Energy = (signal * Gain) * (S / SCluster);
344 Energy += rand.Gaus(0.0, var_level * Energy);
345 if (Energy < 0.0)
346 Energy = 0.0;
347
348 if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
349 cluster.AddStrip(NumCurrentZone, Energy);
350 total_signal += Energy;
351 }
352
353 dist += h;
354
355 } else {
356 // left border strip
357 Int_t NumCurrentZone = (Int_t)LeftZonePos;
358
359 h = ((Int_t)LeftZonePos + 1) - LeftZonePos;
360 if (h < 0.0)
361 h = 0.0;
362
363 Double_t xp = LeftZonePos + dist;
364 Double_t S = gausF.Integral(xp, xp + h);
365 Double_t Energy = (signal * Gain) * (S / SCluster);
366 Energy += rand.Gaus(0.0, var_level * Energy);
367 if (Energy < 0.0)
368 Energy = 0.0;
369
370 if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
371 cluster.AddStrip(NumCurrentZone, Energy);
372 total_signal += Energy;
373 }
374
375 dist += h;
376
377 // inner strips
378 for (Int_t i = (Int_t)LeftZonePos + 1; i < (Int_t)RightZonePos; ++i) {
379
380 NumCurrentZone = i;
381
382 h = 1.0;
383
384 xp = LeftZonePos + dist;
385 S = gausF.Integral(xp, xp + h);
386 Energy = (signal * Gain) * (S / SCluster);
387 Energy += rand.Gaus(0.0, var_level * Energy);
388 if (Energy < 0.0)
389 Energy = 0.0;
390
391 if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
392 cluster.AddStrip(NumCurrentZone, Energy);
393 total_signal += Energy;
394 }
395
396 dist += h;
397 }
398 // right border strip
399 NumCurrentZone = (Int_t)RightZonePos;
400
401 h = RightZonePos - (Int_t)RightZonePos;
402 if (h < 0.0)
403 h = 0.0;
404
405 xp = LeftZonePos + dist;
406 S = gausF.Integral(xp, xp + h);
407 Energy = (signal * Gain) * (S / SCluster);
408 Energy += rand.Gaus(0.0, var_level * Energy);
409 if (Energy < 0.0)
410 Energy = 0.0;
411
412 if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
413 cluster.AddStrip(NumCurrentZone, Energy);
414 total_signal += Energy;
415 }
416
417 dist += h;
418 }
419
420 if (cluster.GetClusterSize() <= 0) {
421 return cluster;
422 }
423
424 // find the mean value of the avalanche position (fitting by gaus function)
425 Double_t mean_fit_pos = 0.0;
426
427 Int_t NOutLeftBins = 0;
428 Int_t NOutRightBins = 0;
429 if (OutLeftBorder != 0.0) {
430 NOutLeftBins = (Int_t)(fabs(OutLeftBorder)) + 1;
431 }
432 if (OutRightBorder != 0.0) {
433 NOutRightBins = (Int_t)(OutRightBorder - RightZonePos) + 1;
434 }
435
436 Int_t begin_strip_num = cluster.Strips.at(0);
437 Int_t last_strip_num = cluster.Strips.at(cluster.GetClusterSize() - 1);
438 Int_t nstrips = last_strip_num - begin_strip_num + 1;
439
440 TH1F hist("hist_for_fit", "hist_for_fit", nstrips + NOutLeftBins + NOutRightBins, begin_strip_num - NOutLeftBins,
441 last_strip_num + 1 + NOutRightBins);
442 Int_t hist_index = 0;
443
444 for (Int_t i = 0; i < cluster.GetClusterSize(); ++i) {
445 Double_t value = cluster.Signals.at(i);
446 hist.SetBinContent(hist_index + 1 + NOutLeftBins, value);
447 hist_index++;
448 }
449
450 // on the left border
451 if (NOutLeftBins > 0.0) {
452 Double_t first = OutLeftBorder;
453 Double_t last = (Int_t)OutLeftBorder;
454 Double_t S = gausF.Integral(first, last);
455 Double_t Energy = (signal * Gain) * (S / SCluster);
456 Energy += rand.Gaus(0.0, var_level * Energy);
457 if (Energy < 0.0)
458 Energy = 0.0;
459 Double_t value = Energy;
460 hist.SetBinContent(1, value);
461 dist = 0.0;
462
463 for (Int_t i = 1; i < NOutLeftBins; ++i) {
464 h = 1.0;
465 first = (Int_t)OutLeftBorder + dist;
466 last = first + h;
467 S = gausF.Integral(first, last);
468 Energy = (signal * Gain) * (S / SCluster);
469 Energy += rand.Gaus(0.0, var_level * Energy);
470 if (Energy < 0.0)
471 Energy = 0.0;
472 value = Energy;
473 hist.SetBinContent(1 + i, value);
474 dist += h;
475 }
476 }
477
478 // on the right border
479 if (NOutRightBins > 0.0) {
480 dist = 0.0;
481 for (Int_t i = hist_index; i < hist_index + NOutRightBins - 1; ++i) {
482 h = 1.0;
483 Double_t first = NStripsInLayer + dist;
484 Double_t last = first + h;
485 Double_t S = gausF.Integral(first, last);
486 Double_t Energy = (signal * Gain) * (S / SCluster);
487 Energy += rand.Gaus(0.0, var_level * Energy);
488 if (Energy < 0.0)
489 Energy = 0.0;
490 Double_t value = Energy;
491 hist.SetBinContent(1 + i, value);
492 dist += h;
493 }
494
495 Double_t first = (Int_t)OutRightBorder;
496 Double_t last = first + (OutRightBorder - (Int_t)OutRightBorder);
497 Double_t S = gausF.Integral(first, last);
498 Double_t Energy = (signal * Gain) * (S / SCluster);
499 Energy += rand.Gaus(0.0, var_level * Energy);
500 if (Energy < 0.0)
501 Energy = 0.0;
502 Double_t value = Energy;
503 hist.SetBinContent(hist_index + NOutRightBins, value);
504 }
505
506 TF1* gausFitFunction = 0;
507 TString fit_params = "WQ0";
508
509#ifdef DRAW_REAL_CLUSTER_HISTOGRAMS
510 fit_params = "WQ";
511#endif
512
513 if (nstrips > 1) {
514 hist.Fit("gaus", fit_params); // Q - quit mode (without information on the screen); 0 - not draw
515 gausFitFunction = hist.GetFunction("gaus");
516 if (gausFitFunction) {
517 mean_fit_pos = gausFitFunction->GetParameter(1);
518 } else {
519 mean_fit_pos = hist.GetMean();
520 }
521 } else {
522 mean_fit_pos = hist.GetMean();
523 }
524
525 cluster.MeanPosition = mean_fit_pos;
526 cluster.TotalSignal = total_signal;
527 cluster.PositionResidual = mean_fit_pos - CenterZonePos; // residual between mean and origin positions (lower
528 // cluster): residual = finded(current) - orig(average)
529
530 cluster.IsCorrect = true; // set correct status of the cluster
531 return cluster;
532}
533
534void BmnSiBTModule::CalculateStripHitIntersectionPoints()
535{
536
537 ResetIntersectionPoints();
538
539 Int_t n_strip_layers = StripLayers.size();
540
541 // Find strip clusters and hits in the strip layers -------------------------
542 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
543 StripLayers[ilayer].FindClustersAndStripHits();
544 }
545 //--------------------------------------------------------------------------
546
547 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
548
549 for (Int_t jlayer = ilayer + 1; jlayer < n_strip_layers; ++jlayer) {
550 // cout << "(i-layer : j-layer) = ( " << ilayer << " : " << jlayer << " )\n";
551
552 if (StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
553 continue;
554
555 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
556 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
557
558 if (Abs(iAngleRad - jAngleRad) < 0.01)
559 continue;
560
561 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
562 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
563
564 for (Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit) {
565 for (Int_t j_strip_hit = 0; j_strip_hit < StripLayers[jlayer].GetNStripHits(); ++j_strip_hit) {
566
567 Double_t iStripHitPos = StripLayers[ilayer].GetStripHitPos(i_strip_hit);
568 Double_t jStripHitPos = StripLayers[jlayer].GetStripHitPos(j_strip_hit);
569
570 Double_t iStripHitError = StripLayers[ilayer].GetStripHitError(i_strip_hit);
571 Double_t jStripHitError = StripLayers[jlayer].GetStripHitError(j_strip_hit);
572
573 Double_t xcoord;
574 Double_t ycoord;
575
576 if (Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01) {
577
578 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
579 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
580
581 xcoord = (jB - iB) / (Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
582 ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
583 } else {
584
585 if (Abs(iAngleRad) < 0.01) {
586
587 if (StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight) {
588 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint()
589 + iStripHitPos * StripLayers[ilayer].GetPitch();
590 } else {
591 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint()
592 - iStripHitPos * StripLayers[ilayer].GetPitch();
593 }
594
595 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
596 ycoord = Tan(j_PiOver2MinusAngleRad) * xcoord + jB;
597 }
598
599 if (Abs(jAngleRad) < 0.01) {
600
601 if (StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight) {
602 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint()
603 + jStripHitPos * StripLayers[jlayer].GetPitch();
604 } else {
605 xcoord = StripLayers[jlayer].GetXRightStripBorderPoint()
606 - jStripHitPos * StripLayers[jlayer].GetPitch();
607 }
608
609 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
610 ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
611 }
612 }
613
614 // check: a point belongs to both strip layers together
615 if (StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord)
616 && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord))
617 {
618
619 IntersectionPointsX.push_back(xcoord);
620 IntersectionPointsY.push_back(ycoord);
621
622 // Intersection (x,y)-errors from the strip-errors ------
623 Double_t AngleIntersecRad = Abs(iAngleRad - jAngleRad);
624
625 // a rhomb zone is the intersection of the strip errors
626 Double_t RhombSide1 = (iStripHitError * StripLayers[ilayer].GetPitch()) / Sin(AngleIntersecRad);
627 Double_t RhombSide2 = (jStripHitError * StripLayers[jlayer].GetPitch()) / Sin(AngleIntersecRad);
628
629 // x-strip error
630 Double_t xerr1 = Sin(Abs(jAngleRad)) * RhombSide1;
631 Double_t xerr2 = Sin(Abs(iAngleRad)) * RhombSide2;
632 Double_t xcoord_err = xerr1 + xerr2;
633
634 // y-strip error
635 Double_t yerr1 = Cos(Abs(jAngleRad)) * RhombSide1;
636 Double_t yerr2 = Cos(Abs(iAngleRad)) * RhombSide2;
637 Double_t ycoord_err = yerr1 + yerr2;
638
639 IntersectionPointsXErrors.push_back(xcoord_err);
640 IntersectionPointsYErrors.push_back(ycoord_err);
641 //------------------------------------------------------
642
643 // intersection matching ----------------------------------------
644 BmnMatch iStripMatch = StripLayers[ilayer].GetStripMatch((Int_t)iStripHitPos);
645 BmnMatch jStripMatch = StripLayers[jlayer].GetStripMatch((Int_t)jStripHitPos);
646
647 BmnMatch intersection_match;
648 intersection_match.AddLink(iStripMatch);
649 intersection_match.AddLink(jStripMatch);
650
651 IntersectionPointMatches.push_back(intersection_match);
652 //--------------------------------------------------------------
653
654 // Additional information about the intersection ------------
655 // cluster size (number strips) in both strip layers for each intersection
656 Int_t iLayerClusterSize = StripLayers[ilayer].GetStripHitClusterSize(i_strip_hit);
657 Int_t jLayerClusterSize = StripLayers[jlayer].GetStripHitClusterSize(j_strip_hit);
658
659 // strip position in both strip layers for each intersection
660 Double_t iLayerStripPosition = StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord);
661 Double_t jLayerStripPosition = StripLayers[jlayer].ConvertPointToStripPosition(xcoord, ycoord);
662
663 IntersectionPoints_LowerLayerStripTotalSignal.push_back(
664 StripLayers[ilayer].GetStripHitTotalSignal(i_strip_hit));
665 IntersectionPoints_UpperLayerStripTotalSignal.push_back(
666 StripLayers[jlayer].GetStripHitTotalSignal(j_strip_hit));
667
668 if (StripLayers[ilayer].GetType() == LowerStripLayer) {
669 IntersectionPoints_LowerLayerClusterSize.push_back(iLayerClusterSize);
670 IntersectionPoints_UpperLayerClusterSize.push_back(jLayerClusterSize);
671 IntersectionPoints_LowerLayerStripPosition.push_back(iLayerStripPosition);
672 IntersectionPoints_UpperLayerStripPosition.push_back(jLayerStripPosition);
673 } else {
674 IntersectionPoints_LowerLayerClusterSize.push_back(jLayerClusterSize);
675 IntersectionPoints_UpperLayerClusterSize.push_back(iLayerClusterSize);
676 IntersectionPoints_LowerLayerStripPosition.push_back(jLayerStripPosition);
677 IntersectionPoints_UpperLayerStripPosition.push_back(iLayerStripPosition);
678 }
679 //------------------------------------------------------
680 }
681 }
682 }
683 }
684 }
685}
686
687// need for a separated test (find x,y intersection coords from strip positions)
688Bool_t BmnSiBTModule::SearchIntersectionPoint(Double_t& x,
689 Double_t& y,
690 Double_t strip_pos_layerA,
691 Double_t strip_pos_layerB,
692 Int_t layerA_index,
693 Int_t layerB_index)
694{
695
696 Int_t ilayer = layerA_index;
697 Int_t jlayer = layerB_index;
698
699 if (StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
700 return false;
701
702 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
703 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
704
705 if (Abs(iAngleRad - jAngleRad) < 0.01)
706 return false;
707
708 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
709 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
710
711 Double_t iStripHitPos = strip_pos_layerA;
712 Double_t jStripHitPos = strip_pos_layerB;
713
714 Double_t xcoord;
715 Double_t ycoord;
716
717 if (Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01) {
718
719 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
720 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
721
722 xcoord = (jB - iB) / (Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
723 ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
724 } else {
725
726 if (Abs(iAngleRad) < 0.01) {
727
728 if (StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight) {
729 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + iStripHitPos * StripLayers[ilayer].GetPitch();
730 } else {
731 xcoord =
732 StripLayers[ilayer].GetXRightStripBorderPoint() - iStripHitPos * StripLayers[ilayer].GetPitch();
733 }
734
735 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
736 ycoord = Tan(j_PiOver2MinusAngleRad) * xcoord + jB;
737 }
738
739 if (Abs(jAngleRad) < 0.01) {
740
741 if (StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight) {
742 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + jStripHitPos * StripLayers[jlayer].GetPitch();
743 } else {
744 xcoord =
745 StripLayers[jlayer].GetXRightStripBorderPoint() - jStripHitPos * StripLayers[jlayer].GetPitch();
746 }
747
748 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
749 ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
750 }
751 }
752
753 x = xcoord;
754 y = ycoord;
755
756 // check: a point belongs to both strip layers together
757 if (StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord)
758 && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord))
759 {
760 return true;
761 }
762
763 return false;
764}
765
766void BmnSiBTModule::ResetIntersectionPoints()
767{
768 IntersectionPointsX.clear();
769 IntersectionPointsY.clear();
770 IntersectionPointsXErrors.clear();
771 IntersectionPointsYErrors.clear();
772 IntersectionPoints_LowerLayerClusterSize.clear();
773 IntersectionPoints_UpperLayerClusterSize.clear();
774 IntersectionPoints_LowerLayerStripPosition.clear();
775 IntersectionPoints_UpperLayerStripPosition.clear();
776 IntersectionPoints_LowerLayerStripTotalSignal.clear();
777 IntersectionPoints_UpperLayerStripTotalSignal.clear();
778 IntersectionPointMatches.clear();
779}
780
781void BmnSiBTModule::DefineModuleBorders()
782{
783
784 if (StripLayers.size() == 0)
785 return;
786
787 XMinModule = 1.0E+10;
788 XMaxModule = -1.0E+10;
789 YMinModule = 1.0E+10;
790 YMaxModule = -1.0E+10;
791
792 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
793 if (StripLayers[ilayer].GetXMinLayer() < XMinModule)
794 XMinModule = StripLayers[ilayer].GetXMinLayer();
795 if (StripLayers[ilayer].GetXMaxLayer() > XMaxModule)
796 XMaxModule = StripLayers[ilayer].GetXMaxLayer();
797 if (StripLayers[ilayer].GetYMinLayer() < YMinModule)
798 YMinModule = StripLayers[ilayer].GetYMinLayer();
799 if (StripLayers[ilayer].GetYMaxLayer() > YMaxModule)
800 YMaxModule = StripLayers[ilayer].GetYMaxLayer();
801 }
802 return;
803}
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
Sigma
Double_t Sigma(vector< Double_t > dist, vector< Double_t > w)
Definition BmnMath.cxx:890
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
LeftToRight
@ LeftToRight
Definition BmnStripData.h:17
LowerStripLayer
@ LowerStripLayer
Definition BmnStripData.h:12
BmnMatch::AddLink
void AddLink(const BmnMatch &match)
Definition BmnMatch.cxx:43
BmnSiBTModule::IsPointInsideZThickness
Bool_t IsPointInsideZThickness(Double_t z)
Definition BmnSiBTModule.cxx:129
BmnSiBTModule::MakeCluster
StripCluster MakeCluster(Int_t layer_num, Double_t xcoord, Double_t ycoord, Double_t signal, Double_t radius)
Definition BmnSiBTModule.cxx:263
BmnSiBTModule::IsPointInsideModule
Bool_t IsPointInsideModule(Double_t x, Double_t y)
Definition BmnSiBTModule.cxx:107
BmnSiBTModule::AddRealPointFullOne
Bool_t AddRealPointFullOne(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID)
Definition BmnSiBTModule.cxx:142
BmnSiBTModule::ResetIntersectionPoints
void ResetIntersectionPoints()
Definition BmnSiBTModule.cxx:766
BmnSiBTModule::AddStripSignalInLayer
Bool_t AddStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
Definition BmnSiBTModule.cxx:63
BmnSiBTModule::SetStripMatchInLayer
Bool_t SetStripMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch strip_match)
Definition BmnSiBTModule.cxx:71
BmnSiBTModule::GetZPositionRegistered
Double_t GetZPositionRegistered()
Definition BmnSiBTModule.cxx:42
BmnSiBTModule::GetStripSignalInLayer
Double_t GetStripSignalInLayer(Int_t layer_num, Int_t strip_num)
Definition BmnSiBTModule.cxx:79
BmnSiBTModule::SetStripSignalInLayer
Bool_t SetStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
Definition BmnSiBTModule.cxx:55
BmnSiBTModule::AddRealPointSimple
Bool_t AddRealPointSimple(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID)
Definition BmnSiBTModule.cxx:222
BmnSiBTModule::AddStripLayer
void AddStripLayer(BmnSiBTLayer strip_layer)
Definition BmnSiBTModule.cxx:48
BmnSiBTModule::~BmnSiBTModule
virtual ~BmnSiBTModule()
Definition BmnSiBTModule.cxx:40
BmnSiBTModule::SearchIntersectionPoint
Bool_t SearchIntersectionPoint(Double_t &x, Double_t &y, Double_t strip_pos_layerA, Double_t strip_pos_layerB, Int_t layerA_index, Int_t layerB_index)
Definition BmnSiBTModule.cxx:688
BmnSiBTModule::ResetRealPoints
void ResetRealPoints()
Definition BmnSiBTModule.h:90
BmnSiBTModule::GetStripMatchInLayer
BmnMatch GetStripMatchInLayer(Int_t layer_num, Int_t strip_num)
Definition BmnSiBTModule.cxx:87
BmnSiBTModule::CalculateStripHitIntersectionPoints
void CalculateStripHitIntersectionPoints()
Definition BmnSiBTModule.cxx:534
StripCluster::TotalSignal
Double_t TotalSignal
Definition BmnStripData.h:181
StripCluster::GetClusterSize
Int_t GetClusterSize()
Definition BmnStripData.h:295
StripCluster::MeanPosition
Double_t MeanPosition
Definition BmnStripData.h:180
StripCluster::AddStrip
void AddStrip(Int_t strip_num, Double_t strip_signal)
Definition BmnStripData.h:237
StripCluster::PositionResidual
Double_t PositionResidual
Definition BmnStripData.h:182
StripCluster::Strips
vector< Int_t > Strips
Definition BmnStripData.h:192
StripCluster::Signals
vector< Double_t > Signals
Definition BmnStripData.h:193