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