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