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