BmnRoot
Loading...
Searching...
No Matches
BmnSiliconModule.cxx
Go to the documentation of this file.
1#include "BmnSiliconModule.h"
2
3#include "TMath.h"
4#include "TRandom.h"
5#include "TF1.h"
6#include "TH1F.h"
7
8#include <iostream>
9
10using namespace TMath;
11using namespace std;
12
14
15 Verbosity = true;
16
17 ZStartModulePosition = 0.0;
18 ModuleThickness = 0.03; //cm
19
20 ModuleRotationAlgleDeg = 0.0;
21 ModuleRotationAlgleRad = ModuleRotationAlgleDeg*TMath::DegToRad();
22 ModuleRotationCenterX = 0.0;
23 ModuleRotationCenterY = 0.0;
24
25 XMinModule = 0.0;
26 XMaxModule = 0.0;
27 YMinModule = 0.0;
28 YMaxModule = 0.0;
29}
30
32
33 Verbosity = true;
34
35 ZStartModulePosition = z_start_pos;
36 ModuleThickness = 0.03; //cm
37
38 ModuleRotationAlgleDeg = 0.0;
39 ModuleRotationAlgleRad = ModuleRotationAlgleDeg*TMath::DegToRad();
40 ModuleRotationCenterX = 0.0;
41 ModuleRotationCenterY = 0.0;
42
43 XMinModule = 0.0;
44 XMaxModule = 0.0;
45 YMinModule = 0.0;
46 YMaxModule = 0.0;
47}
48
52
54 return ZStartModulePosition + ModuleThickness*0.5;
55}
56
57Bool_t BmnSiliconModule::SetModuleRotation(Double_t angleDeg, Double_t rot_center_x, Double_t rot_center_y) {
58 ModuleRotationAlgleDeg = angleDeg;
59 ModuleRotationAlgleRad = ModuleRotationAlgleDeg*TMath::DegToRad();
60 ModuleRotationCenterX = rot_center_x;
61 ModuleRotationCenterY = rot_center_y;
62 return true;
63}
64
66 StripLayers.push_back(strip_layer);
67
68 DefineModuleBorders();
69}
70
71Bool_t BmnSiliconModule::SetStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal) {
72 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
73 return StripLayers[layer_num].SetStripSignal(strip_num, signal);
74
75 return false;
76}
77
78Bool_t BmnSiliconModule::AddStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal) {
79 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
80 return StripLayers[layer_num].AddStripSignal(strip_num, signal);
81
82 return false;
83}
84
85Bool_t BmnSiliconModule::SetStripMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch mc_match) {
86 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
87 return StripLayers[layer_num].SetStripMatch(strip_num, mc_match);
88
89 return false;
90}
91
92Bool_t BmnSiliconModule::SetStripDigitNumberMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch digit_num_match) {
93 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
94 return StripLayers[layer_num].SetStripDigitNumberMatch(strip_num, digit_num_match);
95
96 return false;
97}
98
99Bool_t BmnSiliconModule::SetStripSignalInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t signal) {
100 Int_t n_layers_the_same_id = 0;
101 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
102 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
103 if(StripLayers[ilayer].SetStripSignal(strip_num, signal)) n_layers_the_same_id++;
104 }
105 }
106 return (Bool_t)n_layers_the_same_id;
107}
108
109Bool_t BmnSiliconModule::AddStripSignalInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t signal) {
110 Int_t n_layers_the_same_id = 0;
111 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
112 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
113 if(StripLayers[ilayer].AddStripSignal(strip_num, signal)) n_layers_the_same_id++;
114 }
115 }
116 return (Bool_t)n_layers_the_same_id;
117}
118
119Bool_t BmnSiliconModule::SetStripMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, BmnMatch mc_match) {
120 Int_t n_layers_the_same_id = 0;
121 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
122 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
123 if(StripLayers[ilayer].SetStripMatch(strip_num, mc_match)) n_layers_the_same_id++;
124 }
125 }
126 return (Bool_t)n_layers_the_same_id;
127}
128
129Bool_t BmnSiliconModule::AddStripMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t weight, Int_t mc_num) {
130 Int_t n_layers_the_same_id = 0;
131 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
132 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
133 if(StripLayers[ilayer].AddLinkToStripMatch(strip_num, weight, mc_num)) n_layers_the_same_id++;
134 }
135 }
136 return (Bool_t)n_layers_the_same_id;
137}
138
139Bool_t BmnSiliconModule::SetStripDigitNumberMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, BmnMatch digit_num_match) {
140 Int_t n_layers_the_same_id = 0;
141 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
142 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
143 if(StripLayers[ilayer].SetStripDigitNumberMatch(strip_num, digit_num_match)) n_layers_the_same_id++;
144 }
145 }
146 return (Bool_t)n_layers_the_same_id;
147}
148
149Bool_t BmnSiliconModule::AddStripDigitNumberMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t weight, Int_t digit_num) {
150 Int_t n_layers_the_same_id = 0;
151 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
152 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
153 if(StripLayers[ilayer].AddLinkToStripDigitNumberMatch(strip_num, weight, digit_num)) n_layers_the_same_id++;
154 }
155 }
156 return (Bool_t)n_layers_the_same_id;
157}
158
159Double_t BmnSiliconModule::GetStripSignalInLayer(Int_t layer_num, Int_t strip_num) {
160 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
161 return StripLayers[layer_num].GetStripSignal(strip_num);
162
163 return -1;
164}
165
166BmnMatch BmnSiliconModule::GetStripMatchInLayer(Int_t layer_num, Int_t strip_num) {
167 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
168 return StripLayers[layer_num].GetStripMatch(strip_num);
169
170 return BmnMatch();
171}
172
174 if (layer_num >= 0 && layer_num < (Int_t)StripLayers.size())
175 return StripLayers[layer_num].GetStripDigitNumberMatch(strip_num);
176
177 return BmnMatch();
178}
179
181 Int_t first_strip = 1E6;
182
183 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
184 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
185 Int_t first_strip_curr_layer = StripLayers[ilayer].GetFirstStripNumber();
186 if(first_strip > first_strip_curr_layer) first_strip = first_strip_curr_layer;
187 }
188 }
189 return first_strip;
190}
191
193 Int_t last_strip = 0;
194
195 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
196 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
197 Int_t last_strip_curr_layer = StripLayers[ilayer].GetLastStripNumber();
198 if(last_strip < last_strip_curr_layer) last_strip = last_strip_curr_layer;
199 }
200 }
201 return last_strip;
202}
203
204Double_t BmnSiliconModule::GetStripSignalInZone(Int_t zone_id, Int_t strip_num)
205{
206 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
207 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
208 if( strip_num >= StripLayers[ilayer].GetFirstStripNumber() && strip_num <= StripLayers[ilayer].GetLastStripNumber() ) {
209 return StripLayers[ilayer].GetStripSignal(strip_num);
210 }
211 }
212 }
213 return -1;
214}
215
216BmnMatch BmnSiliconModule::GetStripMatchInZone(Int_t zone_id, Int_t strip_num)
217{
218 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
219 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
220 if( strip_num >= StripLayers[ilayer].GetFirstStripNumber() && strip_num <= StripLayers[ilayer].GetLastStripNumber() ) {
221 return StripLayers[ilayer].GetStripMatch(strip_num);
222 }
223 }
224 }
225 return BmnMatch();
226}
227
229{
230 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
231 if( (StripLayers[ilayer].GetZoneID() == zone_id) ) {
232 if( strip_num >= StripLayers[ilayer].GetFirstStripNumber() && strip_num <= StripLayers[ilayer].GetLastStripNumber() ) {
233 return StripLayers[ilayer].GetStripDigitNumberMatch(strip_num);
234 }
235 }
236 }
237 return BmnMatch();
238}
239
241 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer)
242 StripLayers[ilayer].ResetLayer();
243
244 DefineModuleBorders();
245
249}
250
251Bool_t BmnSiliconModule::IsPointInsideModule(Double_t x, Double_t y, Bool_t isLocal = false) {
252
253 Double_t xloc = x;
254 Double_t yloc = y;
255
256 if(!isLocal) {
257 xloc = ConvertRotatedToLocalX(x, y);
258 yloc = ConvertRotatedToLocalY(x, y);
259 }
260
261 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer)
262 if( StripLayers[ilayer].IsPointInsideStripLayer(xloc, yloc) ) return true;
263
264 return false;
265}
266
267Bool_t BmnSiliconModule::IsPointInsideModule(Double_t x, Double_t y, Double_t z, Bool_t isLocal = false) {
268 Double_t coord_eps = 0.01; //100 um
269 if( ((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps)) &&
270 ((z < (ZStartModulePosition + ModuleThickness)) || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps) ) {
271 if( IsPointInsideModule(x, y, isLocal) ) return true;
272 }
273 return false;
274}
275
277 Double_t coord_eps = 0.01; //100 um
278 if( ((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps)) &&
279 ((z < (ZStartModulePosition + ModuleThickness)) || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps) ) {
280 return true;
281 }
282 return false;
283}
284
285//Add a single point without smearing
286Bool_t BmnSiliconModule::AddRealPointSimple(Double_t x, Double_t y, Double_t z,
287 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
288
289 Double_t x_loc = ConvertRotatedToLocalX(x, y);
290 Double_t y_loc = ConvertRotatedToLocalY(x, y);
291 //Double_t z_loc = z;
292
293 Double_t px_loc = ConvertRotatedToLocalX(px, py);
294 Double_t py_loc = ConvertRotatedToLocalY(px, py);
295 Double_t pz_loc = pz;
296
297 //Effective XY position at the middle of the detector thickness (dz/2) -----
298
299 //Distance from entry point (x,y,z) to the midde of the detector (along z-axis)
300 Double_t z_dist_mid = ModuleThickness*0.5;
301 Double_t vector_coeff_mid = fabs(z_dist_mid/pz_loc); //for middle point
302
303 Double_t x_vec_from_in_to_mid = vector_coeff_mid*px_loc;
304 Double_t y_vec_from_in_to_mid = vector_coeff_mid*py_loc;
305 //Double_t z_vec_from_in_to_mid = vector_coeff_mid*pz_loc;
306
307 //Middle point coordinates (x_mid, y_mid, z_mid)
308 Double_t x_mid = x_loc + x_vec_from_in_to_mid;
309 Double_t y_mid = y_loc + y_vec_from_in_to_mid;
310 //Double_t z_mid = z_loc + z_vec_from_in_to_mid;
311 //--------------------------------------------------------------------------
312
313 if( IsPointInsideModule(x, y) ) {
314
315 Int_t n_strip_layers = StripLayers.size();
316
317 //Add a signal to the strips layers ------------------------------------
318 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
319 if( !StripLayers[ilayer].IsPointInsideStripLayer(x_mid, y_mid) ) continue;
320
321 Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_mid, y_mid);
322 Int_t zone_id = StripLayers[ilayer].GetZoneID();
323
324 AddStripSignalInLayerByZoneId(zone_id, (Int_t)strip_pos, signal);
325
326 //strip match
327 AddStripMatchInLayerByZoneId(zone_id, (Int_t)strip_pos, signal, refID);
328
329 }
330 //----------------------------------------------------------------------
331
332 RealPointsX.push_back(x);
333 RealPointsY.push_back(y);
334 RealPointsMC.push_back(refID);
335
336 return true;
337 }
338 else {
339 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointSimple(): Point (" << x << " : " << y << ") is out of the readout plane or inside a dead zone\n";
340 return false;
341 }
342}
343
344//Add single point with Gaussian smearing
345Bool_t BmnSiliconModule::AddRealPointFullOne(Double_t x, Double_t y, Double_t z,
346 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
347
348 Double_t xloc = ConvertRotatedToLocalX(x, y);
349 Double_t yloc = ConvertRotatedToLocalY(x, y);
350
351 if( IsPointInsideModule(x, y) ) {
352
353 if(signal <= 0.0) signal = 1e-16;
354
355 Double_t AvalancheRadius = 0.01; //cm
356
357 Double_t radius = AvalancheRadius;
358 if(radius <= 0.0) return false;
359
360 //AZ-260322 gRandom->SetSeed(0);
361
362 Int_t cycle_cnt = 0;
363 while(true) {
364 radius = gRandom->Gaus(AvalancheRadius, 0.005);
365 if(radius > AvalancheRadius/2.0 && radius < AvalancheRadius*2.0 && radius > 0.0) break;
366 cycle_cnt++;
367
368 if(cycle_cnt > 5) {
369 radius = AvalancheRadius;
370 break;
371 }
372 }
373
374 //Make strip cluster in each strip layer -------------------------------
375 Int_t n_strip_layers = StripLayers.size();
376 vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
377
378 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
379 if( !StripLayers[ilayer].IsPointInsideStripLayer(xloc, yloc) ) continue;
380 cluster_layers[ilayer] = MakeCluster(ilayer, xloc, yloc, signal, radius);
381 }
382 //----------------------------------------------------------------------
383
384 //Add the correct clusters to the strip layers -------------------------
385 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
386 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
387 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
388 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
389 Int_t zone_id = StripLayers[ilayer].GetZoneID();
390 AddStripSignalInLayerByZoneId(zone_id, strip_num, strip_signal);
391 }
392 }
393 //----------------------------------------------------------------------
394
395 //Fill strip matches ---------------------------------------------------
396 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
397 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
398 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
399 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
400 Int_t zone_id = StripLayers[ilayer].GetZoneID();
401 AddStripMatchInLayerByZoneId(zone_id, strip_num, strip_signal, refID);
402 }
403 }
404 //----------------------------------------------------------------------
405
406 RealPointsX.push_back(x);
407 RealPointsY.push_back(y);
408 RealPointsMC.push_back(refID);
409
410 return true;
411 }
412 else {
413 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointFullOne(): Point (" << x << " : " << y << " : " << z << ") is out of the readout plane or inside a dead zone\n";
414 return false;
415 }
416}
417
418//Add single point with Gaussian smearing + taking into account particle track incline
419Bool_t BmnSiliconModule::AddRealPointFullOne_WithIncline(Double_t x, Double_t y, Double_t z,
420 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
421
422 Double_t x_loc = ConvertRotatedToLocalX(x, y);
423 Double_t y_loc = ConvertRotatedToLocalY(x, y);
424 //Double_t z_loc = z;
425
426 Double_t px_loc = ConvertRotatedToLocalX(px, py);
427 Double_t py_loc = ConvertRotatedToLocalY(px, py);
428 Double_t pz_loc = pz;
429
430 //Effective XY position at the middle of the detector thickness (dz/2) -----
431
432 //Distance from entry point (x,y,z) to the midde of the detector (along z-axis)
433 Double_t z_dist_mid = ModuleThickness*0.5;
434 Double_t vector_coeff_mid = fabs(z_dist_mid/pz_loc); //for middle point
435
436 Double_t x_vec_from_in_to_mid = vector_coeff_mid*px_loc;
437 Double_t y_vec_from_in_to_mid = vector_coeff_mid*py_loc;
438 //Double_t z_vec_from_in_to_mid = vector_coeff_mid*pz_loc;
439
440 //Middle point coordinates (x_mid, y_mid, z_mid)
441 Double_t x_mid = x_loc + x_vec_from_in_to_mid;
442 Double_t y_mid = y_loc + y_vec_from_in_to_mid;
443 //Double_t z_mid = z_loc + z_vec_from_in_to_mid;
444 //--------------------------------------------------------------------------
445
446 if( IsPointInsideModule(x, y) ) {
447
448 if(signal <= 0.0) signal = 1e-16;
449
450 Double_t AvalancheRadius = 0.015; //cm
451
452 Double_t radius = AvalancheRadius;
453 if(radius <= 0.0) return false;
454
455 //AZ-260322 gRandom->SetSeed(0);
456
457 Int_t cycle_cnt = 0;
458 while(true) {
459 radius = gRandom->Gaus(AvalancheRadius, 0.005);
460 if(radius > AvalancheRadius/2.0 && radius < AvalancheRadius*2.0 && radius > 0.0) break;
461 cycle_cnt++;
462
463 if(cycle_cnt > 5) {
464 radius = AvalancheRadius;
465 break;
466 }
467 }
468
469 //Make strip cluster in each strip layer -------------------------------
470 Int_t n_strip_layers = StripLayers.size();
471 vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
472
473 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
474 if( !StripLayers[ilayer].IsPointInsideStripLayer(x_mid, y_mid) ) continue;
475 cluster_layers[ilayer] = MakeCluster(ilayer, x_mid, y_mid, signal, radius);
476 }
477 //----------------------------------------------------------------------
478
479 //Add the correct clusters to the strip layers -------------------------
480 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
481 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
482 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
483 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
484 Int_t zone_id = StripLayers[ilayer].GetZoneID();
485 AddStripSignalInLayerByZoneId(zone_id, strip_num, strip_signal);
486 }
487 }
488 //----------------------------------------------------------------------
489
490 //Fill strip matches ---------------------------------------------------
491 for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
492 for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
493 Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
494 Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
495 Int_t zone_id = StripLayers[ilayer].GetZoneID();
496 AddStripMatchInLayerByZoneId(zone_id, strip_num, strip_signal, refID);
497 }
498 }
499 //----------------------------------------------------------------------
500
501 RealPointsX.push_back(x);
502 RealPointsY.push_back(y);
503 RealPointsMC.push_back(refID);
504
505 return true;
506 }
507 else {
508 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointFullOne(): Point (" << x << " : " << y << " : " << z << ") is out of the readout plane or inside a dead zone\n";
509 return false;
510 }
511}
512
513Bool_t BmnSiliconModule::AddRealPointFull(Double_t x, Double_t y, Double_t z,
514 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
515
516 //Convert coodinates if rotation around Z axis is used ---------------------
517 Double_t x_loc = ConvertRotatedToLocalX(x, y);
518 Double_t y_loc = ConvertRotatedToLocalY(x, y);
519 Double_t z_loc = z;
520
521 Double_t px_loc = ConvertRotatedToLocalX(px, py);
522 Double_t py_loc = ConvertRotatedToLocalY(px, py);
523 Double_t pz_loc = pz;
524 //--------------------------------------------------------------------------
525
526 //Parameters of the signal simulation of one charge carrier (electron or hole)
527 //Double_t chargedCarrierCloudRadius = 0.005; //Signal radius (not diameter!) for ONE STEP (cm) ~50 micron (for large step)
528 Double_t chargedCarrierCloudRadius = gRandom->Uniform(0.0001, 0.0005); //Signal radius (not diameter!) for ONE STEP (cm) 1-5 micron (for small step)
529 Double_t step = 0.0010; //simulation step (is not mean collision dist. for semiconductors!)
530 //--------------------------------------------------------------------------
531
532 if( IsPointInsideModule(x_loc, y_loc, true) ) {
533
534 //Distance from entry point (x,y,z) to exit point (along z-axis)
535 Double_t z_distance_from_in_to_out;
536
537 //Condition: track direction along z-axis
538 if(pz_loc > 0.0) {
539 z_distance_from_in_to_out = ModuleThickness - (z_loc-ZStartModulePosition);
540 }
541 else {
542 z_distance_from_in_to_out = z_loc - ZStartModulePosition;
543 }
544
545 //Condition: track distance along z-axis must be not zero
546 if(z_distance_from_in_to_out <= 0.0) {
547 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointFull(): Point (" << x << " : " << y << " : " << z << ") has a track with an incorrect length\n";
548 return false;
549 }
550
551 //Scale coefficient of the track vector (the path from entry point to exit point)
552 Double_t vector_coeff = fabs(z_distance_from_in_to_out/pz_loc);
553
554 //Components of the track vector
555 Double_t x_vec_from_in_to_out = vector_coeff*px_loc;
556 Double_t y_vec_from_in_to_out = vector_coeff*py_loc;
557 Double_t z_vec_from_in_to_out = vector_coeff*pz_loc;
558
559 //Exit point coordinates (x_out, y_out, z_out)
560 Double_t x_out = x_loc + x_vec_from_in_to_out;
561 Double_t y_out = y_loc + y_vec_from_in_to_out;
562 //Double_t z_out = z_loc + z_vec_from_in_to_out;
563
564 //Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
565 //Condition: x-coordinate of the exit point is inside the module
566 if(x_out < XMinModule || x_out > XMaxModule ) {
567 if(x_out < XMinModule ) x_out = XMinModule;
568 if(x_out > XMaxModule ) x_out = XMaxModule;
569
570 x_vec_from_in_to_out = x_out - x_loc;
571 vector_coeff = x_vec_from_in_to_out/px_loc;
572
573 y_vec_from_in_to_out= vector_coeff*py_loc;
574 z_vec_from_in_to_out = vector_coeff*pz_loc;
575
576 y_out = y_loc + y_vec_from_in_to_out;
577 //z_out = z_loc + z_vec_from_in_to_out;
578 }
579 //Condition: y-coordinate of the exit point is inside the module
580 if(y_out < YMinModule || y_out > YMaxModule ) {
581 if(y_out < YMinModule ) y_out = YMinModule;
582 if(y_out > YMaxModule ) y_out = YMaxModule;
583
584 y_vec_from_in_to_out = y_out - y_loc;
585 vector_coeff = y_vec_from_in_to_out/py_loc;
586
587 x_vec_from_in_to_out= vector_coeff*px_loc;
588 z_vec_from_in_to_out = vector_coeff*pz_loc;
589
590 x_out = x_loc + x_vec_from_in_to_out;
591 //z_out = z_loc + z_vec_from_in_to_out;
592 }
593 //----------------------------------------------------------------------
594
595 //Common track length (Fix: a square root!)
596 Double_t track_length = std::sqrt( (x_vec_from_in_to_out)*(x_vec_from_in_to_out) + (y_vec_from_in_to_out)*(y_vec_from_in_to_out) + (z_vec_from_in_to_out)*(z_vec_from_in_to_out) );
597
598 Double_t current_length = 0.0; //traversed length (counter)
599 Double_t current_length_ratio = 0.0; //ratio of the traversed length to common track length (counter)
600
601 Int_t collisions_cnt = 0; //total collision counter
602 Double_t current_step = 0.0; //current distance between two collision points
603
604 //Collection of collision points
605 std::vector<CollPoint> collision_points;
606
607 while(current_length < track_length) {
608
609 //current_step = rand.Exp(step);
610 current_step = step;
611 current_length += current_step;
612
613 if(current_length > track_length) break;
614
615 current_length_ratio = current_length/track_length;
616
617 Double_t current_x = x_loc + current_length_ratio*x_vec_from_in_to_out;
618 Double_t current_y = y_loc + current_length_ratio*y_vec_from_in_to_out;
619 Double_t current_z = z_loc + current_length_ratio*z_vec_from_in_to_out;
620
621 collision_points.push_back(CollPoint(current_x, current_y, current_z));
622 collisions_cnt++;
623 }
624 //----------------------------------------------------------------------
625
626 if(signal <= 0.0) signal = 0;
627
628 Int_t n_collisions = collision_points.size();
629 Double_t signal_per_step = signal / n_collisions;
630
631 for(Int_t icollpoint = 0; icollpoint < n_collisions; ++icollpoint) {
632
633 //generate the current signal radius on strips for the current charged carrier
634 Double_t current_radius = chargedCarrierCloudRadius;
635 if(chargedCarrierCloudRadius <= 0.0) {
636 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointFull(): incorrect chargedCarrierCloudRadius!\n";
637 return false;
638 }
639
640 Double_t current_x = collision_points[icollpoint].x;
641 Double_t current_y = collision_points[icollpoint].y;
642
643 Int_t n_strip_layers = GetNStripLayers();
644
645 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
646 if( !StripLayers[ilayer].IsPointInsideStripLayer(current_x, current_y) ) continue;
647
648 Double_t pitch = StripLayers[ilayer].GetPitch();
649 Double_t radius_in_strips = current_radius/pitch;
650 Int_t zone_id = StripLayers[ilayer].GetZoneID();
651
652 Int_t first_strip_in_layer = StripLayers[ilayer].GetFirstStripNumber();
653 Int_t last_strip_in_layer = StripLayers[ilayer].GetLastStripNumber();
654
655
656 Double_t central_strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(current_x, current_y);
657 Double_t left_strip_pos = central_strip_pos - radius_in_strips;
658 Double_t right_strip_pos = central_strip_pos + radius_in_strips;
659
660 //cout << " layer: " << ilayer << "\n";
661 //cout << " current_x = " << current_x << "\n";
662 //cout << " current_y = " << current_y << "\n";
663 //cout << " central_strip_pos = " << central_strip_pos << "\n";
664
665 Double_t central_strip_signal = signal_per_step;
666 Double_t left_strip_signal = 0.0;
667 Double_t right_strip_signal = 0.0;
668
669 Double_t central_strip_fraction = 1.0;
670 Double_t left_strip_fraction = 0.0;
671 Double_t right_strip_fraction = 0.0;
672
673 //calculate signal fraction for the left neighbour strip
674 if( (Int_t)left_strip_pos != (Int_t)central_strip_pos ) {
675 if( left_strip_pos >= first_strip_in_layer ) {
676 left_strip_fraction = (Int_t)central_strip_pos - left_strip_pos;
677 if(left_strip_fraction > 1.0) left_strip_fraction = 1.0;
678 }
679 }
680 //calculate signal fraction for the right neighbour strip
681 if( (Int_t)right_strip_pos != (Int_t)central_strip_pos ) {
682 if( right_strip_pos < (last_strip_in_layer+1) ) {
683 right_strip_fraction = right_strip_pos - (Int_t)(central_strip_pos+1);
684 if(right_strip_fraction > 1.0) right_strip_fraction = 1.0;
685 }
686 }
687
688 //signal for the triplet strips (one central and two possible neighbour strips)
689 central_strip_signal = signal_per_step*central_strip_fraction / (central_strip_fraction + left_strip_fraction + right_strip_fraction);
690 left_strip_signal = signal_per_step*left_strip_fraction / (central_strip_fraction + left_strip_fraction + right_strip_fraction);
691 right_strip_signal = signal_per_step*right_strip_fraction / (central_strip_fraction + left_strip_fraction + right_strip_fraction);
692
693 AddStripSignalInLayerByZoneId(zone_id, (Int_t)central_strip_pos, central_strip_signal);
694 if(left_strip_signal > 0.0) AddStripSignalInLayerByZoneId(zone_id, (Int_t)left_strip_pos, left_strip_signal);
695 if(right_strip_signal > 0.0) AddStripSignalInLayerByZoneId(zone_id, (Int_t)right_strip_pos, right_strip_signal);
696
697 //strip match
698 AddStripMatchInLayerByZoneId(zone_id, (Int_t)central_strip_pos, central_strip_signal, refID);
699 if(left_strip_signal > 0.0) AddStripMatchInLayerByZoneId(zone_id, (Int_t)left_strip_pos, left_strip_signal, refID);
700 if(right_strip_signal > 0.0) AddStripMatchInLayerByZoneId(zone_id, (Int_t)right_strip_pos, right_strip_signal, refID);
701 }
702 }
703
704 RealPointsX.push_back(x);
705 RealPointsY.push_back(y);
706 RealPointsMC.push_back(refID);
707
708 return true;
709 }
710 else {
711 if(Verbosity) cerr << "WARNING: BmnSiliconModule::AddRealPointFull(): Point (" << x << " : " << y << " : " << z << ") is out of the readout plane or inside a dead zone\n";
712 return false;
713 }
714}
715
716StripCluster BmnSiliconModule::MakeCluster(Int_t layer_num, Double_t xcoord, Double_t ycoord, Double_t signal, Double_t radius) {
717
718 Double_t ClusterDistortion = 0.0; //signal variation (0.1 is 10%)
719
720 Double_t Gain = 1.0;
721
722 StripCluster cluster;
723
724 if(radius <= 0.0) return cluster;
725
726 Double_t Pitch = StripLayers[layer_num].GetPitch();
727
728 Double_t RadiusInZones = radius/Pitch;
729 Double_t Sigma = RadiusInZones/3.33;
730
731 TF1 gausF("gausF", "gaus", -4*Sigma, 4*Sigma);
732 gausF.SetParameter(0, 1.0); // constant (altitude)
733 gausF.SetParameter(1, 0.0); // mean (center position)
734 gausF.SetParameter(2, Sigma); //sigma
735
736 Double_t SCluster = gausF.Integral(-4*Sigma, 4*Sigma); //square of the one side distribution (more exactly)
737
738 TRandom rand(0);
739 Double_t var_level = ClusterDistortion; //signal variation (0.1 is 10%)
740
741 Int_t NStripsInLayer = StripLayers[layer_num].GetNStrips();
742 Int_t FirstStripInLayer = StripLayers[layer_num].GetFirstStripNumber();
743 Int_t LastStripInLayer = StripLayers[layer_num].GetLastStripNumber();
744
745 if(NStripsInLayer == 0) return cluster;
746
747 Double_t CenterZonePos = StripLayers[layer_num].ConvertPointToStripPosition(xcoord, ycoord);
748
749 if( CenterZonePos < FirstStripInLayer || CenterZonePos >= LastStripInLayer+1 ) return cluster;
750
751
752 gausF.SetParameter(1, CenterZonePos);
753 Double_t total_signal = 0.0;
754
755 Double_t LeftZonePos = CenterZonePos - RadiusInZones;
756 Double_t RightZonePos = CenterZonePos + RadiusInZones;
757 Double_t OutLeftBorder = 0.0;
758 Double_t OutRightBorder = 0.0;
759
760 if(LeftZonePos < FirstStripInLayer) {
761 OutLeftBorder = LeftZonePos;
762 LeftZonePos = FirstStripInLayer;
763 }
764 if(RightZonePos < FirstStripInLayer) {
765 OutRightBorder = RightZonePos;
766 RightZonePos = FirstStripInLayer;
767 }
768 if(LeftZonePos >= LastStripInLayer+1) {
769 OutLeftBorder = LeftZonePos;
770 LeftZonePos = LastStripInLayer+1 - 0.001;
771 }
772 if(RightZonePos >= LastStripInLayer+1) {
773 OutRightBorder = RightZonePos;
774 RightZonePos = LastStripInLayer+1 - 0.001;
775 }
776
777 Double_t h = 0.0;
778 Double_t dist = 0.0;
779
780 //avalanche is inside of the one strip
781 if((Int_t)LeftZonePos == (Int_t)RightZonePos) {
782
783 Int_t NumCurrentZone = (Int_t) LeftZonePos;
784
785 h = RightZonePos - LeftZonePos;
786 if(h < 0.0) h = 0.0;
787
788 Double_t xp = LeftZonePos + dist;
789 Double_t S = gausF.Integral(xp, xp+h);
790 Double_t Energy = (signal*Gain)*(S/SCluster);
791 Energy += rand.Gaus(0.0, var_level*Energy);
792 if(Energy < 0.0) Energy = 0.0;
793
794 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
795 cluster.AddStrip(NumCurrentZone, Energy);
796 total_signal += Energy;
797 }
798
799 dist += h;
800
801 }
802 else {
803 //left border strip
804 Int_t NumCurrentZone = (Int_t) LeftZonePos;
805
806 h = ((Int_t)LeftZonePos + 1) - LeftZonePos;
807 if(h < 0.0) h = 0.0;
808
809 Double_t xp = LeftZonePos + dist;
810 Double_t S = gausF.Integral(xp, xp+h);
811 Double_t Energy = (signal*Gain)*(S/SCluster);
812 Energy += rand.Gaus(0.0, var_level*Energy);
813 if(Energy < 0.0) Energy = 0.0;
814
815 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
816 cluster.AddStrip(NumCurrentZone, Energy);
817 total_signal += Energy;
818 }
819
820 dist += h;
821
822 //inner strips
823 for(Int_t i = (Int_t)LeftZonePos + 1; i < (Int_t)RightZonePos; ++i) {
824
825 NumCurrentZone = i;
826
827 h = 1.0;
828
829 xp = LeftZonePos + dist;
830 S = gausF.Integral(xp, xp+h);
831 Energy = (signal*Gain)*(S/SCluster);
832 Energy += rand.Gaus(0.0, var_level*Energy);
833 if(Energy < 0.0) Energy = 0.0;
834
835 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
836 cluster.AddStrip(NumCurrentZone, Energy);
837 total_signal += Energy;
838 }
839
840 dist += h;
841 }
842 //right border strip
843 NumCurrentZone = (Int_t)RightZonePos;
844
845 h = RightZonePos - (Int_t)RightZonePos;
846 if(h < 0.0) h = 0.0;
847
848 xp = LeftZonePos + dist;
849 S = gausF.Integral(xp, xp+h);
850 Energy = (signal*Gain)*(S/SCluster);
851 Energy += rand.Gaus(0.0, var_level*Energy);
852 if(Energy < 0.0) Energy = 0.0;
853
854 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
855 cluster.AddStrip(NumCurrentZone, Energy);
856 total_signal += Energy;
857 }
858
859 dist += h;
860 }
861
862 if (cluster.GetClusterSize() <= 0) {
863 return cluster;
864 }
865
866 //find the mean value of the avalanche position (fitting by gaus function)
867 Double_t mean_fit_pos = 0.0;
868
869 Int_t NOutLeftBins = 0;
870 Int_t NOutRightBins = 0;
871 if(OutLeftBorder != 0.0) {
872 NOutLeftBins = (Int_t)(fabs(OutLeftBorder)) + 1;
873 }
874 if(OutRightBorder != 0.0) {
875 NOutRightBins = (Int_t)(OutRightBorder - RightZonePos) + 1;
876 }
877
878 Int_t begin_strip_num = cluster.Strips.at(0);
879 Int_t last_strip_num = cluster.Strips.at(cluster.GetClusterSize()-1);
880 Int_t nstrips = last_strip_num - begin_strip_num + 1;
881
882 TH1F hist("hist_for_fit", "hist_for_fit", nstrips+NOutLeftBins+NOutRightBins, begin_strip_num-NOutLeftBins, last_strip_num+1+NOutRightBins);
883 Int_t hist_index = 0;
884
885 for(Int_t i = 0; i < cluster.GetClusterSize(); ++i) {
886 Double_t value = cluster.Signals.at(i);
887 hist.SetBinContent(hist_index+1+NOutLeftBins, value);
888 hist_index++;
889 }
890
891 //on the left border
892 if(NOutLeftBins > 0.0) {
893 Double_t first = OutLeftBorder;
894 Double_t last = (Int_t)OutLeftBorder;
895 Double_t S = gausF.Integral(first, last);
896 Double_t Energy = (signal*Gain)*(S/SCluster);
897 Energy += rand.Gaus(0.0, var_level*Energy);
898 if(Energy < 0.0) Energy = 0.0;
899 Double_t value = Energy;
900 hist.SetBinContent(1, value);
901 dist = 0.0;
902
903 for(Int_t i = 1; i < NOutLeftBins; ++i) {
904 h = 1.0;
905 first = (Int_t)OutLeftBorder+dist;
906 last = first + h;
907 S = gausF.Integral(first, last);
908 Energy = (signal*Gain)*(S/SCluster);
909 Energy += rand.Gaus(0.0, var_level*Energy);
910 if(Energy < 0.0) Energy = 0.0;
911 value = Energy;
912 hist.SetBinContent(1+i, value);
913 dist += h;
914 }
915 }
916
917 //on the right border
918 if(NOutRightBins > 0.0) {
919 dist = 0.0;
920 for(Int_t i = hist_index; i < hist_index+NOutRightBins-1; ++i) {
921 h = 1.0;
922 Double_t first = NStripsInLayer+dist;
923 Double_t last = first + h;
924 Double_t S = gausF.Integral(first, last);
925 Double_t Energy = (signal*Gain)*(S/SCluster);
926 Energy += rand.Gaus(0.0, var_level*Energy);
927 if(Energy < 0.0) Energy = 0.0;
928 Double_t value = Energy;
929 hist.SetBinContent(1+i, value);
930 dist += h;
931 }
932
933 Double_t first = (Int_t)OutRightBorder;
934 Double_t last = first + (OutRightBorder - (Int_t)OutRightBorder);
935 Double_t S = gausF.Integral(first, last);
936 Double_t Energy = (signal*Gain)*(S/SCluster);
937 Energy += rand.Gaus(0.0, var_level*Energy);
938 if(Energy < 0.0) Energy = 0.0;
939 Double_t value = Energy;
940 hist.SetBinContent(hist_index+NOutRightBins, value);
941 }
942
943 TF1* gausFitFunction = 0;
944 TString fit_params = "WQ0";
945
946 #ifdef DRAW_REAL_CLUSTER_HISTOGRAMS
947 fit_params = "WQ";
948 #endif
949
950 if(nstrips > 1) {
951 hist.Fit("gaus", fit_params); //Q - quit mode (without information on the screen); 0 - not draw
952 gausFitFunction = hist.GetFunction("gaus");
953 if(gausFitFunction) {
954 mean_fit_pos = gausFitFunction->GetParameter(1);
955 }
956 else {
957 mean_fit_pos = hist.GetMean();
958 }
959 }
960 else {
961 mean_fit_pos = hist.GetMean();
962 }
963
964 cluster.MeanPosition = mean_fit_pos;
965 cluster.TotalSignal = total_signal;
966 cluster.PositionResidual = mean_fit_pos - CenterZonePos; //residual between mean and origin positions (lower cluster): residual = finded(current) - orig(average)
967
968 cluster.IsCorrect = true; // set correct status of the cluster
969 return cluster;
970}
971
973
976
977 Int_t n_strip_layers = StripLayers.size();
978
979 //Find strip clusters and hits in the strip layers -------------------------
980 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
981 StripLayers[ilayer].FindClustersAndStripHits();
982 }
983 //--------------------------------------------------------------------------
984
985 map<Int_t, Bool_t> UsageStatus_LowerClusters;
986 map<Int_t, Bool_t> UsageStatus_UpperClusters;
987
988 //Saving all unique id marks of clusters in a module
989 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
990 for(Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit) {
991 StripCluster cluster = StripLayers[ilayer].GetStripClusters()[i_strip_hit];
992 if(cluster.GetType() == LowerStripLayer) {
993 UsageStatus_LowerClusters[cluster.GetUniqueID()] = false;
994 }
995 if(cluster.GetType() == UpperStripLayer) {
996 UsageStatus_UpperClusters[cluster.GetUniqueID()] = false;
997 }
998 }
999 }
1000
1001 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
1002 for(Int_t jlayer = ilayer+1; jlayer < n_strip_layers; ++jlayer) {
1003 //cout << "(i-layer : j-layer) = ( " << ilayer << " : " << jlayer << " )\n";
1004
1005 if(StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType()) continue;
1006
1007 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
1008 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
1009
1010 if( Abs(iAngleRad - jAngleRad) < 0.01 ) continue;
1011
1012 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
1013 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
1014
1015 for(Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit) {
1016 for(Int_t j_strip_hit = 0; j_strip_hit < StripLayers[jlayer].GetNStripHits(); ++j_strip_hit) {
1017
1018 Double_t iStripHitPos = StripLayers[ilayer].GetStripHitPos(i_strip_hit);
1019 Double_t jStripHitPos = StripLayers[jlayer].GetStripHitPos(j_strip_hit);
1020
1021 Double_t iStripHitError = StripLayers[ilayer].GetStripHitError(i_strip_hit);
1022 Double_t jStripHitError = StripLayers[jlayer].GetStripHitError(j_strip_hit);
1023
1024 Double_t xcoord;
1025 Double_t ycoord;
1026
1027 if( Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01 ) {
1028
1029 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1030 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1031
1032 xcoord = (jB - iB)/(Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
1033 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1034 }
1035 else {
1036
1037 if( Abs(iAngleRad) < 0.01 ) {
1038
1039 if(StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight ) {
1040 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1041 }
1042 else {
1043 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1044 }
1045
1046 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1047 ycoord = Tan(j_PiOver2MinusAngleRad)*xcoord + jB;
1048 }
1049
1050 if( Abs(jAngleRad) < 0.01 ) {
1051
1052 if(StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight ) {
1053 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1054 }
1055 else {
1056 xcoord = StripLayers[jlayer].GetXRightStripBorderPoint() - (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1057 }
1058
1059 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1060 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1061 }
1062 }
1063
1064 //check: a point belongs to both strip layers together
1065 if( StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord) && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord) ) {
1066
1067 Double_t xrot = ConvertLocalToRotatedX(xcoord, ycoord);
1068 Double_t yrot = ConvertLocalToRotatedY(xcoord, ycoord);
1069
1070 IntersectionPointsX.push_back(xrot);
1071 IntersectionPointsY.push_back(yrot);
1072
1073 //Intersection (x,y)-errors from the strip-errors ------
1074 Double_t AngleIntersecRad = Abs(iAngleRad - jAngleRad);
1075
1076 //a rhomb zone is the intersection of the strip errors
1077 Double_t RhombSide1 = (iStripHitError*StripLayers[ilayer].GetPitch())/Sin(AngleIntersecRad);
1078 Double_t RhombSide2 = (jStripHitError*StripLayers[jlayer].GetPitch())/Sin(AngleIntersecRad);
1079
1080 //x-strip error
1081 Double_t xerr1 = Sin(Abs(jAngleRad))*RhombSide1;
1082 Double_t xerr2 = Sin(Abs(iAngleRad))*RhombSide2;
1083 Double_t xcoord_err = xerr1 + xerr2;
1084
1085 //y-strip error
1086 Double_t yerr1 = Cos(Abs(jAngleRad))*RhombSide1;
1087 Double_t yerr2 = Cos(Abs(iAngleRad))*RhombSide2;
1088 Double_t ycoord_err = yerr1 + yerr2;
1089
1090 IntersectionPointsXErrors.push_back(xcoord_err);
1091 IntersectionPointsYErrors.push_back(ycoord_err);
1092 //------------------------------------------------------
1093
1094 //intersection MC-matching -----------------------------
1095 BmnMatch iStripMatch = StripLayers[ilayer].GetStripMatch((Int_t)iStripHitPos);
1096 BmnMatch jStripMatch = StripLayers[jlayer].GetStripMatch((Int_t)jStripHitPos);
1097
1098 BmnMatch intersection_match;
1099 intersection_match.AddLink(iStripMatch);
1100 intersection_match.AddLink(jStripMatch);
1101
1102 IntersectionPointMatches.push_back(intersection_match);
1103 //------------------------------------------------------
1104
1105 //intersection digit number matching -------------------
1106 BmnMatch iStripDigitNumberMatch = StripLayers[ilayer].GetStripDigitNumberMatch((Int_t)iStripHitPos);
1107 BmnMatch jStripDigitNumberMatch = StripLayers[jlayer].GetStripDigitNumberMatch((Int_t)jStripHitPos);
1108
1109 BmnMatch intersection_digit_num_match;
1110 intersection_digit_num_match.AddLink(iStripDigitNumberMatch);
1111 intersection_digit_num_match.AddLink(jStripDigitNumberMatch);
1112
1113 IntersectionPointDigitNumberMatches.push_back(intersection_digit_num_match);
1114 //------------------------------------------------------
1115
1116 //Additional information about the intersection --------
1117 //cluster size (number strips) in both strip layers for each intersection
1118 Int_t iLayerClusterSize = StripLayers[ilayer].GetStripHitClusterSize(i_strip_hit);
1119 Int_t jLayerClusterSize = StripLayers[jlayer].GetStripHitClusterSize(j_strip_hit);
1120
1121 //strip position in both strip layers for each intersection
1122 Double_t iLayerStripPosition = StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord);
1123 Double_t jLayerStripPosition = StripLayers[jlayer].ConvertPointToStripPosition(xcoord, ycoord);
1124
1125 if(StripLayers[ilayer].GetType() == LowerStripLayer) {
1126 IntersectionPoints_LowerLayerClusterSize.push_back(iLayerClusterSize);
1127 IntersectionPoints_UpperLayerClusterSize.push_back(jLayerClusterSize);
1128 IntersectionPoints_LowerLayerStripPosition.push_back(iLayerStripPosition);
1129 IntersectionPoints_UpperLayerStripPosition.push_back(jLayerStripPosition);
1130 IntersectionPoints_LowerLayerStripTotalSignal.push_back(StripLayers[ilayer].GetStripHitTotalSignal(i_strip_hit));
1131 IntersectionPoints_UpperLayerStripTotalSignal.push_back(StripLayers[jlayer].GetStripHitTotalSignal(j_strip_hit));
1132 LowerClusters.push_back(StripLayers[ilayer].GetStripClusters().at(i_strip_hit));
1133 UpperClusters.push_back(StripLayers[jlayer].GetStripClusters().at(j_strip_hit));
1134 UsageStatus_LowerClusters[StripLayers[ilayer].GetStripClusters()[i_strip_hit].GetUniqueID()] = true;
1135 UsageStatus_UpperClusters[StripLayers[jlayer].GetStripClusters()[j_strip_hit].GetUniqueID()] = true;
1136 }
1137 else {
1138 IntersectionPoints_LowerLayerClusterSize.push_back(jLayerClusterSize);
1139 IntersectionPoints_UpperLayerClusterSize.push_back(iLayerClusterSize);
1140 IntersectionPoints_LowerLayerStripPosition.push_back(jLayerStripPosition);
1141 IntersectionPoints_UpperLayerStripPosition.push_back(iLayerStripPosition);
1142 IntersectionPoints_LowerLayerStripTotalSignal.push_back(StripLayers[jlayer].GetStripHitTotalSignal(j_strip_hit));
1143 IntersectionPoints_UpperLayerStripTotalSignal.push_back(StripLayers[ilayer].GetStripHitTotalSignal(i_strip_hit));
1144 LowerClusters.push_back(StripLayers[jlayer].GetStripClusters().at(j_strip_hit));
1145 UpperClusters.push_back(StripLayers[ilayer].GetStripClusters().at(i_strip_hit));
1146 UsageStatus_LowerClusters[StripLayers[jlayer].GetStripClusters()[j_strip_hit].GetUniqueID()] = true;
1147 UsageStatus_UpperClusters[StripLayers[ilayer].GetStripClusters()[i_strip_hit].GetUniqueID()] = true;
1148 }
1149 //------------------------------------------------------
1150 }
1151 }
1152 }
1153 }
1154 }
1155
1156 /*for (auto it : UsageStatus_UpperClusters)
1157 if(it.second == false) cout << " upper_cluster[" << it.first << "] = " << it.second << "\n";
1158 */
1159
1160 for(auto it : UsageStatus_LowerClusters)
1161 {
1162 //if(it.second == false) cout << " lower_cluster[" << it.first << "] = " << it.second << "\n";
1163 if(it.second == false)
1164 {
1165 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
1166 {
1167 if(StripLayers[ilayer].GetType() != LowerStripLayer) continue;
1168
1169 for (size_t icluster = 0; icluster < StripLayers[ilayer].GetStripClusters().size(); ++icluster) {
1170 StripCluster cluster = StripLayers[ilayer].GetStripClusters()[icluster];
1171
1172 if((int)cluster.GetUniqueID() != it.first) continue;
1173
1174 Double_t xcoord;
1175 Double_t ycoord;
1176
1177 if(StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight ) {
1178 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + cluster.MeanPosition*StripLayers[ilayer].GetPitch();
1179 }
1180 else {
1181 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - cluster.MeanPosition*StripLayers[ilayer].GetPitch();
1182 }
1183
1184 ycoord = (StripLayers[ilayer].GetYMinLayer() + StripLayers[ilayer].GetYMaxLayer())*0.5;
1185
1186 Double_t xerr = cluster.Error*StripLayers[ilayer].GetPitch();
1187 Double_t yerr = StripLayers[ilayer].GetYSize()*0.5;
1188
1189 //cout << " low_cluster[" << icluster << "]:\n";
1190 //cout << " ID = " << cluster.GetUniqueID() << "\n";
1191 //cout << " xcoord = " << xcoord << "\n";
1192 //cout << " xmin = " << StripLayers[ilayer].GetXMinLayer() << ", xmax = " << StripLayers[ilayer].GetXMaxLayer() << "\n";
1193 //cout << " ycoord = " << ycoord << "\n";
1194 //cout << " ymin = " << StripLayers[ilayer].GetYMinLayer() << ", ymax = " << StripLayers[ilayer].GetYMaxLayer() << "\n";
1195 //cout << " xerr = " << xerr << "\n";
1196 //cout << " yerr = " << yerr << "\n";
1197 //cout << "\n";
1198
1199 PseudoIntersectionsX.push_back(xcoord);
1200 PseudoIntersectionsY.push_back(ycoord);
1201
1202 PseudoIntersectionsXErrors.push_back(xerr);
1203 PseudoIntersectionsYErrors.push_back(yerr);
1204
1205 PseudoIntersections_LowerLayerClusterSize.push_back(cluster.Width);
1206 PseudoIntersections_UpperLayerClusterSize.push_back(0);
1207
1208 PseudoIntersections_LowerLayerStripPosition.push_back(StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord));
1209 PseudoIntersections_UpperLayerStripPosition.push_back(0);
1210
1211 PseudoIntersections_LowerLayerStripTotalSignal.push_back(cluster.TotalSignal);
1212 PseudoIntersections_UpperLayerStripTotalSignal.push_back(0);
1213
1214 PseudoIntersectionMatches.push_back(StripLayers[ilayer].GetStripMatch((Int_t)cluster.MeanPosition)); //mc-match
1215 PseudoIntersectionDigitNumberMatches.push_back(StripLayers[ilayer].GetStripDigitNumberMatch((Int_t)cluster.MeanPosition)); //digit number matches
1216
1217 LowerClusters_PseudoIntersections.push_back(cluster);
1218 UpperClusters_PseudoIntersections.push_back(StripCluster()); //empty cluster with default parameters
1219 }
1220 }
1221 }
1222 }
1223}
1224
1225//need for a separated test (find x,y intersection coords from strip positions)
1226Bool_t BmnSiliconModule::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) {
1227
1228 Int_t ilayer = layerA_index;
1229 Int_t jlayer = layerB_index;
1230
1231 if(StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType()) return false;
1232
1233 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
1234 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
1235
1236 if( Abs(iAngleRad - jAngleRad) < 0.01 ) return false;
1237
1238 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
1239 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
1240
1241 Double_t iStripHitPos = strip_pos_layerA;
1242 Double_t jStripHitPos = strip_pos_layerB;
1243
1244 Double_t xcoord;
1245 Double_t ycoord;
1246
1247 if( Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01 ) {
1248
1249 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1250 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1251
1252 xcoord = (jB - iB)/(Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
1253 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1254 }
1255 else {
1256
1257 if( Abs(iAngleRad) < 0.01 ) {
1258
1259 if(StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight ) {
1260 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1261 }
1262 else {
1263 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1264 }
1265
1266 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1267 ycoord = Tan(j_PiOver2MinusAngleRad)*xcoord + jB;
1268 }
1269
1270 if( Abs(jAngleRad) < 0.01 ) {
1271
1272 if(StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight ) {
1273 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1274 }
1275 else {
1276 xcoord = StripLayers[jlayer].GetXRightStripBorderPoint() - (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1277 }
1278
1279 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1280 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1281 }
1282 }
1283
1284 Double_t xrot = ConvertLocalToRotatedX(xcoord, ycoord);
1285 Double_t yrot = ConvertLocalToRotatedY(xcoord, ycoord);
1286
1287 x = xrot;
1288 y = yrot;
1289
1290 //check: a point belongs to both strip layers together
1291 if( StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord) && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord) ) {
1292 return true;
1293 }
1294
1295 return false;
1296}
1297
1298Double_t BmnSiliconModule::ConvertRotatedToLocalX(Double_t xrot, Double_t yrot) {
1299 //Getting local X-coordinate by applying an inverse rotation to rotated coordinates
1300 Double_t xloc = ModuleRotationCenterX + (xrot-ModuleRotationCenterX)*cos(ModuleRotationAlgleRad) - (yrot-ModuleRotationCenterY)*sin(ModuleRotationAlgleRad);
1301 return xloc;
1302}
1303
1304Double_t BmnSiliconModule::ConvertRotatedToLocalY(Double_t xrot, Double_t yrot) {
1305 //Getting local X-coordinate by applying an inverse rotation to rotated coordinates
1306 Double_t yloc = ModuleRotationCenterY + (xrot-ModuleRotationCenterX)*sin(ModuleRotationAlgleRad) + (yrot-ModuleRotationCenterY)*cos(ModuleRotationAlgleRad);
1307 return yloc;
1308}
1309
1310Double_t BmnSiliconModule::ConvertLocalToRotatedX(Double_t xloc, Double_t yloc) {
1311 //Getting rotated X-coordinate by applying a rotation to local (non-rotated) coordinates
1312 Double_t xrot = ModuleRotationCenterX + (xloc-ModuleRotationCenterX)*cos(-ModuleRotationAlgleRad) - (yloc-ModuleRotationCenterY)*sin(-ModuleRotationAlgleRad);
1313 return xrot;
1314}
1315
1316Double_t BmnSiliconModule::ConvertLocalToRotatedY(Double_t xloc, Double_t yloc) {
1317 //Getting rotated Y-coordinate by applying a rotation to local (non-rotated) coordinates
1318 Double_t yrot = ModuleRotationCenterY + (xloc-ModuleRotationCenterX)*sin(-ModuleRotationAlgleRad) + (yloc-ModuleRotationCenterY)*cos(-ModuleRotationAlgleRad);
1319 return yrot;
1320}
1321
1323 IntersectionPointsX.clear();
1324 IntersectionPointsY.clear();
1325 IntersectionPointsXErrors.clear();
1326 IntersectionPointsYErrors.clear();
1327 IntersectionPoints_LowerLayerClusterSize.clear();
1328 IntersectionPoints_UpperLayerClusterSize.clear();
1329 IntersectionPoints_LowerLayerStripPosition.clear();
1330 IntersectionPoints_UpperLayerStripPosition.clear();
1331 IntersectionPoints_LowerLayerStripTotalSignal.clear();
1332 IntersectionPoints_UpperLayerStripTotalSignal.clear();
1333 IntersectionPointMatches.clear();
1334 IntersectionPointDigitNumberMatches.clear();
1335 UpperClusters.clear();
1336 LowerClusters.clear();
1337}
1338
1340 PseudoIntersectionsX.clear();
1341 PseudoIntersectionsY.clear();
1342 PseudoIntersectionsXErrors.clear();
1343 PseudoIntersectionsYErrors.clear();
1344 PseudoIntersections_LowerLayerClusterSize.clear();
1345 PseudoIntersections_UpperLayerClusterSize.clear();
1346 PseudoIntersections_LowerLayerStripPosition.clear();
1347 PseudoIntersections_UpperLayerStripPosition.clear();
1348 PseudoIntersections_LowerLayerStripTotalSignal.clear();
1349 PseudoIntersections_UpperLayerStripTotalSignal.clear();
1350 PseudoIntersectionMatches.clear();
1351 PseudoIntersectionDigitNumberMatches.clear();
1352 UpperClusters_PseudoIntersections.clear();
1353 LowerClusters_PseudoIntersections.clear();
1354 }
1355
1356void BmnSiliconModule::DefineModuleBorders() {
1357
1358 if( StripLayers.size() == 0 ) return;
1359
1360 XMinModule = 1.0E+10;
1361 XMaxModule = -1.0E+10;
1362 YMinModule = 1.0E+10;
1363 YMaxModule = -1.0E+10;
1364
1365 for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
1366 if( StripLayers[ilayer].GetXMinLayer() < XMinModule ) XMinModule = StripLayers[ilayer].GetXMinLayer();
1367 if( StripLayers[ilayer].GetXMaxLayer() > XMaxModule ) XMaxModule = StripLayers[ilayer].GetXMaxLayer();
1368 if( StripLayers[ilayer].GetYMinLayer() < YMinModule ) YMinModule = StripLayers[ilayer].GetYMinLayer();
1369 if( StripLayers[ilayer].GetYMaxLayer() > YMaxModule ) YMaxModule = StripLayers[ilayer].GetYMaxLayer();
1370 }
1371 return;
1372}
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
Double_t Sigma(vector< Double_t > dist, vector< Double_t > w)
Definition BmnMath.cxx:890
friend F32vec4 sin(const F32vec4 &a)
Definition P4_F32vec4.h:124
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
friend F32vec4 cos(const F32vec4 &a)
Definition P4_F32vec4.h:125
int i
Definition P4_F32vec4.h:22
@ LeftToRight
@ UpperStripLayer
@ LowerStripLayer
void AddLink(const BmnMatch &match)
Definition BmnMatch.cxx:43
Bool_t SetStripMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, BmnMatch mc_match)
BmnMatch GetStripMatchInLayer(Int_t layer_num, Int_t strip_num)
Bool_t IsPointInsideZThickness(Double_t z)
Bool_t SetStripDigitNumberMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, BmnMatch digit_num_match)
Double_t GetStripSignalInZone(Int_t zone_id, Int_t strip_num)
Double_t ConvertLocalToRotatedX(Double_t xloc, Double_t yloc)
Bool_t AddRealPointFull(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID)
Int_t GetLastStripInZone(Int_t zone_id)
Bool_t SetStripMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch mc_match)
void AddStripLayer(BmnSiliconLayer strip_layer)
Double_t ConvertLocalToRotatedY(Double_t xloc, Double_t yloc)
BmnMatch GetStripDigitNumberMatchInZone(Int_t zone_id, Int_t strip_num)
Bool_t SetModuleRotation(Double_t angleDeg, Double_t rot_center_x, Double_t rot_center_y)
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)
Double_t GetZPositionRegistered()
BmnMatch GetStripMatchInZone(Int_t zone_id, Int_t strip_num)
Bool_t AddStripMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t weight, Int_t mc_num)
Double_t GetStripSignalInLayer(Int_t layer_num, Int_t strip_num)
Bool_t SetStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
Bool_t SetStripDigitNumberMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch digit_num_match)
Bool_t AddRealPointFullOne_WithIncline(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID)
Bool_t AddStripDigitNumberMatchInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t weight, Int_t digit_num)
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)
Bool_t AddStripSignalInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t signal)
Double_t ConvertRotatedToLocalX(Double_t xrot, Double_t yrot)
Bool_t IsPointInsideModule(Double_t x, Double_t y, Bool_t isLocal)
StripCluster MakeCluster(Int_t layer_num, Double_t xcoord, Double_t ycoord, Double_t signal, Double_t radius)
BmnMatch GetStripDigitNumberMatchInLayer(Int_t layer_num, Int_t strip_num)
void CalculateStripHitIntersectionPoints()
Bool_t SetStripSignalInLayerByZoneId(Int_t zone_id, Int_t strip_num, Double_t signal)
Int_t GetFirstStripInZone(Int_t zone_id)
Bool_t AddStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
Double_t ConvertRotatedToLocalY(Double_t xrot, Double_t yrot)
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)
Double_t TotalSignal
Int_t GetClusterSize()
Int_t GetType()
Double_t MeanPosition
Double_t Error
void AddStrip(Int_t strip_num, Double_t strip_signal)
Double_t PositionResidual
vector< Int_t > Strips
vector< Double_t > Signals
STL namespace.