347 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
354 if(signal <= 0.0) signal = 1e-16;
356 Double_t AvalancheRadius = 0.01;
358 Double_t radius = AvalancheRadius;
359 if(radius <= 0.0)
return false;
365 radius = gRandom->Gaus(AvalancheRadius, 0.005);
366 if(radius > AvalancheRadius/2.0 && radius < AvalancheRadius*2.0 && radius > 0.0)
break;
370 radius = AvalancheRadius;
376 Int_t n_strip_layers = StripLayers.size();
377 vector<StripCluster> cluster_layers(n_strip_layers,
StripCluster());
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);
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();
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();
407 RealPointsX.push_back(x);
408 RealPointsY.push_back(y);
409 RealPointsMC.push_back(refID);
414 if(Verbosity) cerr <<
"WARNING: BmnVSPModule::AddRealPointFullOne(): Point (" << x <<
" : " << y <<
" : " << z <<
") is out of the readout plane or inside a dead zone\n";
421 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
429 Double_t pz_loc = pz;
434 Double_t z_dist_mid = ModuleThickness*0.5;
435 Double_t vector_coeff_mid =
fabs(z_dist_mid/pz_loc);
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;
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;
449 if(signal <= 0.0) signal = 1e-16;
451 Double_t AvalancheRadius = 0.015;
453 Double_t radius = AvalancheRadius;
454 if(radius <= 0.0)
return false;
460 radius = gRandom->Gaus(AvalancheRadius, 0.005);
461 if(radius > AvalancheRadius/2.0 && radius < AvalancheRadius*2.0 && radius > 0.0)
break;
465 radius = AvalancheRadius;
471 Int_t n_strip_layers = StripLayers.size();
472 vector<StripCluster> cluster_layers(n_strip_layers,
StripCluster());
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);
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();
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();
502 RealPointsX.push_back(x);
503 RealPointsY.push_back(y);
504 RealPointsMC.push_back(refID);
509 if(Verbosity) cerr <<
"WARNING: BmnVSPModule::AddRealPointFullOne(): Point (" << x <<
" : " << y <<
" : " << z <<
") is out of the readout plane or inside a dead zone\n";
515 Double_t px, Double_t py, Double_t pz, Double_t signal, Int_t refID) {
524 Double_t pz_loc = pz;
529 Double_t chargedCarrierCloudRadius = gRandom->Uniform(0.0001, 0.0005);
530 Double_t step = 0.0010;
536 Double_t z_distance_from_in_to_out;
540 z_distance_from_in_to_out = ModuleThickness - (z_loc-ZStartModulePosition);
543 z_distance_from_in_to_out = z_loc - ZStartModulePosition;
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";
553 Double_t vector_coeff =
fabs(z_distance_from_in_to_out/pz_loc);
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;
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;
567 if(x_out < XMinModule || x_out > XMaxModule ) {
568 if(x_out < XMinModule ) x_out = XMinModule;
569 if(x_out > XMaxModule ) x_out = XMaxModule;
571 x_vec_from_in_to_out = x_out - x_loc;
572 vector_coeff = x_vec_from_in_to_out/px_loc;
574 y_vec_from_in_to_out= vector_coeff*py_loc;
575 z_vec_from_in_to_out = vector_coeff*pz_loc;
577 y_out = y_loc + y_vec_from_in_to_out;
581 if(y_out < YMinModule || y_out > YMaxModule ) {
582 if(y_out < YMinModule ) y_out = YMinModule;
583 if(y_out > YMaxModule ) y_out = YMaxModule;
585 y_vec_from_in_to_out = y_out - y_loc;
586 vector_coeff = y_vec_from_in_to_out/py_loc;
588 x_vec_from_in_to_out= vector_coeff*px_loc;
589 z_vec_from_in_to_out = vector_coeff*pz_loc;
591 x_out = x_loc + x_vec_from_in_to_out;
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) );
599 Double_t current_length = 0.0;
600 Double_t current_length_ratio = 0.0;
602 Int_t collisions_cnt = 0;
603 Double_t current_step = 0.0;
606 std::vector<CollPoint> collision_points;
608 while(current_length < track_length) {
612 current_length += current_step;
614 if(current_length > track_length)
break;
616 current_length_ratio = current_length/track_length;
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;
622 collision_points.push_back(CollPoint(current_x, current_y, current_z));
627 if(signal <= 0.0) signal = 0;
629 Int_t n_collisions = collision_points.size();
630 Double_t signal_per_step = signal / n_collisions;
632 for(Int_t icollpoint = 0; icollpoint < n_collisions; ++icollpoint) {
635 Double_t current_radius = chargedCarrierCloudRadius;
636 if(chargedCarrierCloudRadius <= 0.0) {
637 if(Verbosity) cerr <<
"WARNING: BmnVSPModule::AddRealPointFull(): incorrect chargedCarrierCloudRadius!\n";
641 Double_t current_x = collision_points[icollpoint].x;
642 Double_t current_y = collision_points[icollpoint].y;
646 for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
647 if( !StripLayers[ilayer].IsPointInsideStripLayer(current_x, current_y) )
continue;
649 Double_t pitch = StripLayers[ilayer].GetPitch();
650 Double_t radius_in_strips = current_radius/pitch;
651 Int_t zone_id = StripLayers[ilayer].GetZoneID();
653 Int_t first_strip_in_layer = StripLayers[ilayer].GetFirstStripNumber();
654 Int_t last_strip_in_layer = StripLayers[ilayer].GetLastStripNumber();
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;
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;
670 Double_t central_strip_fraction = 1.0;
671 Double_t left_strip_fraction = 0.0;
672 Double_t right_strip_fraction = 0.0;
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;
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;
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);
705 RealPointsX.push_back(x);
706 RealPointsY.push_back(y);
707 RealPointsMC.push_back(refID);
712 if(Verbosity) cerr <<
"WARNING: BmnVSPModule::AddRealPointFull(): Point (" << x <<
" : " << y <<
" : " << z <<
") is out of the readout plane or inside a dead zone\n";
719 Double_t ClusterDistortion = 0.0;
725 if(radius <= 0.0)
return cluster;
727 Double_t Pitch = StripLayers[layer_num].GetPitch();
729 Double_t RadiusInZones = radius/Pitch;
730 Double_t
Sigma = RadiusInZones/3.33;
732 TF1 gausF(
"gausF",
"gaus", -4*
Sigma, 4*
Sigma);
733 gausF.SetParameter(0, 1.0);
734 gausF.SetParameter(1, 0.0);
735 gausF.SetParameter(2,
Sigma);
737 Double_t SCluster = gausF.Integral(-4*
Sigma, 4*
Sigma);
740 Double_t var_level = ClusterDistortion;
742 Int_t NStripsInLayer = StripLayers[layer_num].GetNStrips();
743 Int_t FirstStripInLayer = StripLayers[layer_num].GetFirstStripNumber();
744 Int_t LastStripInLayer = StripLayers[layer_num].GetLastStripNumber();
746 if(NStripsInLayer == 0)
return cluster;
748 Double_t CenterZonePos = StripLayers[layer_num].ConvertPointToStripPosition(xcoord, ycoord);
750 if( CenterZonePos < FirstStripInLayer || CenterZonePos >= LastStripInLayer+1 )
return cluster;
753 gausF.SetParameter(1, CenterZonePos);
754 Double_t total_signal = 0.0;
756 Double_t LeftZonePos = CenterZonePos - RadiusInZones;
757 Double_t RightZonePos = CenterZonePos + RadiusInZones;
758 Double_t OutLeftBorder = 0.0;
759 Double_t OutRightBorder = 0.0;
761 if(LeftZonePos < FirstStripInLayer) {
762 OutLeftBorder = LeftZonePos;
763 LeftZonePos = FirstStripInLayer;
765 if(RightZonePos < FirstStripInLayer) {
766 OutRightBorder = RightZonePos;
767 RightZonePos = FirstStripInLayer;
769 if(LeftZonePos >= LastStripInLayer+1) {
770 OutLeftBorder = LeftZonePos;
771 LeftZonePos = LastStripInLayer+1 - 0.001;
773 if(RightZonePos >= LastStripInLayer+1) {
774 OutRightBorder = RightZonePos;
775 RightZonePos = LastStripInLayer+1 - 0.001;
782 if((Int_t)LeftZonePos == (Int_t)RightZonePos) {
784 Int_t NumCurrentZone = (Int_t) LeftZonePos;
786 h = RightZonePos - LeftZonePos;
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;
795 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
796 cluster.
AddStrip(NumCurrentZone, Energy);
797 total_signal += Energy;
805 Int_t NumCurrentZone = (Int_t) LeftZonePos;
807 h = ((Int_t)LeftZonePos + 1) - LeftZonePos;
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;
816 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
817 cluster.
AddStrip(NumCurrentZone, Energy);
818 total_signal += Energy;
824 for(Int_t
i = (Int_t)LeftZonePos + 1;
i < (Int_t)RightZonePos; ++
i) {
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;
836 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
837 cluster.
AddStrip(NumCurrentZone, Energy);
838 total_signal += Energy;
844 NumCurrentZone = (Int_t)RightZonePos;
846 h = RightZonePos - (Int_t)RightZonePos;
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;
855 if(NumCurrentZone >= FirstStripInLayer && NumCurrentZone < LastStripInLayer+1 && Energy > 0.0) {
856 cluster.
AddStrip(NumCurrentZone, Energy);
857 total_signal += Energy;
868 Double_t mean_fit_pos = 0.0;
870 Int_t NOutLeftBins = 0;
871 Int_t NOutRightBins = 0;
872 if(OutLeftBorder != 0.0) {
873 NOutLeftBins = (Int_t)(
fabs(OutLeftBorder)) + 1;
875 if(OutRightBorder != 0.0) {
876 NOutRightBins = (Int_t)(OutRightBorder - RightZonePos) + 1;
879 Int_t begin_strip_num = cluster.
Strips.at(0);
881 Int_t nstrips = last_strip_num - begin_strip_num + 1;
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;
887 Double_t value = cluster.
Signals.at(
i);
888 hist.SetBinContent(hist_index+1+NOutLeftBins, value);
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);
904 for(Int_t
i = 1;
i < NOutLeftBins; ++
i) {
906 first = (Int_t)OutLeftBorder+
dist;
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;
913 hist.SetBinContent(1+
i, value);
919 if(NOutRightBins > 0.0) {
921 for(Int_t
i = hist_index;
i < hist_index+NOutRightBins-1; ++
i) {
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);
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);
944 TF1* gausFitFunction = 0;
945 TString fit_params =
"WQ0";
947 #ifdef DRAW_REAL_CLUSTER_HISTOGRAMS
952 hist.Fit(
"gaus", fit_params);
953 gausFitFunction = hist.GetFunction(
"gaus");
954 if(gausFitFunction) {
955 mean_fit_pos = gausFitFunction->GetParameter(1);
958 mean_fit_pos = hist.GetMean();
962 mean_fit_pos = hist.GetMean();
978 size_t n_strip_layers = StripLayers.size();
981 for (
size_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
982 StripLayers[ilayer].FindClustersAndStripHits();
984 map<Int_t, Bool_t> UsageStatus_LowerClusters;
985 map<Int_t, Bool_t> UsageStatus_UpperClusters;
988 for (
size_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
990 for (Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit)
992 StripCluster cluster = StripLayers[ilayer].GetStripClusters()[i_strip_hit];
994 UsageStatus_LowerClusters[cluster.GetUniqueID()] =
false;
996 UsageStatus_UpperClusters[cluster.GetUniqueID()] =
false;
1000 for(
size_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
1001 for(
size_t jlayer = ilayer+1; jlayer < n_strip_layers; ++jlayer) {
1004 if(StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
continue;
1006 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
1007 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
1009 if( Abs(iAngleRad - jAngleRad) < 0.01 )
continue;
1011 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
1012 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
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) {
1017 Double_t iStripHitPos = StripLayers[ilayer].GetStripHitPos(i_strip_hit);
1018 Double_t jStripHitPos = StripLayers[jlayer].GetStripHitPos(j_strip_hit);
1020 Double_t iStripHitError = StripLayers[ilayer].GetStripHitError(i_strip_hit);
1021 Double_t jStripHitError = StripLayers[jlayer].GetStripHitError(j_strip_hit);
1026 if( Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01 ) {
1028 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1029 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1031 xcoord = (jB - iB)/(Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
1032 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1036 if( Abs(iAngleRad) < 0.01 ) {
1038 if(StripLayers[ilayer].GetStripNumberingOrder() ==
LeftToRight ) {
1039 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1042 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1045 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1046 ycoord = Tan(j_PiOver2MinusAngleRad)*xcoord + jB;
1049 if( Abs(jAngleRad) < 0.01 ) {
1051 if(StripLayers[jlayer].GetStripNumberingOrder() ==
LeftToRight ) {
1052 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1055 xcoord = StripLayers[jlayer].GetXRightStripBorderPoint() - (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1058 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1059 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1064 if( StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord) && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord) ) {
1069 IntersectionPointsX.push_back(xrot);
1070 IntersectionPointsY.push_back(yrot);
1073 Double_t AngleIntersecRad = Abs(iAngleRad - jAngleRad);
1076 Double_t RhombSide1 = (iStripHitError*StripLayers[ilayer].GetPitch())/Sin(AngleIntersecRad);
1077 Double_t RhombSide2 = (jStripHitError*StripLayers[jlayer].GetPitch())/Sin(AngleIntersecRad);
1080 Double_t xerr1 = Sin(Abs(jAngleRad))*RhombSide1;
1081 Double_t xerr2 = Sin(Abs(iAngleRad))*RhombSide2;
1082 Double_t xcoord_err = xerr1 + xerr2;
1085 Double_t yerr1 = Cos(Abs(jAngleRad))*RhombSide1;
1086 Double_t yerr2 = Cos(Abs(iAngleRad))*RhombSide2;
1087 Double_t ycoord_err = yerr1 + yerr2;
1089 IntersectionPointsXErrors.push_back(xcoord_err);
1090 IntersectionPointsYErrors.push_back(ycoord_err);
1094 BmnMatch iStripMatch = StripLayers[ilayer].GetStripMatch((Int_t)iStripHitPos);
1095 BmnMatch jStripMatch = StripLayers[jlayer].GetStripMatch((Int_t)jStripHitPos);
1098 intersection_match.
AddLink(iStripMatch);
1099 intersection_match.
AddLink(jStripMatch);
1101 IntersectionPointMatches.push_back(intersection_match);
1105 BmnMatch iStripDigitNumberMatch = StripLayers[ilayer].GetStripDigitNumberMatch((Int_t)iStripHitPos);
1106 BmnMatch jStripDigitNumberMatch = StripLayers[jlayer].GetStripDigitNumberMatch((Int_t)jStripHitPos);
1108 BmnMatch intersection_digit_num_match;
1109 intersection_digit_num_match.
AddLink(iStripDigitNumberMatch);
1110 intersection_digit_num_match.
AddLink(jStripDigitNumberMatch);
1112 IntersectionPointDigitNumberMatches.push_back(intersection_digit_num_match);
1117 Int_t iLayerClusterSize = StripLayers[ilayer].GetStripHitClusterSize(i_strip_hit);
1118 Int_t jLayerClusterSize = StripLayers[jlayer].GetStripHitClusterSize(j_strip_hit);
1121 Double_t iLayerStripPosition = StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord);
1122 Double_t jLayerStripPosition = StripLayers[jlayer].ConvertPointToStripPosition(xcoord, ycoord);
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;
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;
1159 for(
auto it : UsageStatus_LowerClusters)
1163 if(it.second ==
false) {
1165 for (
size_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
1169 for (
size_t icluster = 0; icluster < StripLayers[ilayer].GetStripClusters().size(); ++icluster)
1171 StripCluster cluster = StripLayers[ilayer].GetStripClusters()[icluster];
1173 if ((Int_t)cluster.GetUniqueID() != it.first)
continue;
1178 if(StripLayers[ilayer].GetStripNumberingOrder() ==
LeftToRight ) {
1179 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + cluster.
MeanPosition*StripLayers[ilayer].GetPitch();
1182 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - cluster.
MeanPosition*StripLayers[ilayer].GetPitch();
1185 ycoord = (StripLayers[ilayer].GetYMinLayer() + StripLayers[ilayer].GetYMaxLayer())*0.5;
1187 Double_t xerr = cluster.
Error*StripLayers[ilayer].GetPitch();
1188 Double_t yerr = StripLayers[ilayer].GetYSize()*0.5;
1200 PseudoIntersectionsX.push_back(xcoord);
1201 PseudoIntersectionsY.push_back(ycoord);
1203 PseudoIntersectionsXErrors.push_back(xerr);
1204 PseudoIntersectionsYErrors.push_back(yerr);
1206 PseudoIntersections_LowerLayerClusterSize.push_back(cluster.
Width);
1207 PseudoIntersections_UpperLayerClusterSize.push_back(0);
1209 PseudoIntersections_LowerLayerStripPosition.push_back(StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord));
1210 PseudoIntersections_UpperLayerStripPosition.push_back(0);
1212 PseudoIntersections_LowerLayerStripTotalSignal.push_back(cluster.
TotalSignal);
1213 PseudoIntersections_UpperLayerStripTotalSignal.push_back(0);
1215 PseudoIntersectionMatches.push_back(StripLayers[ilayer].GetStripMatch((Int_t)cluster.
MeanPosition));
1216 PseudoIntersectionDigitNumberMatches.push_back(StripLayers[ilayer].GetStripDigitNumberMatch((Int_t)cluster.
MeanPosition));
1218 LowerClusters_PseudoIntersections.push_back(cluster);
1219 UpperClusters_PseudoIntersections.push_back(
StripCluster());
1229 Int_t ilayer = layerA_index;
1230 Int_t jlayer = layerB_index;
1232 if(StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
return false;
1234 Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
1235 Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
1237 if( Abs(iAngleRad - jAngleRad) < 0.01 )
return false;
1239 Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
1240 Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
1242 Double_t iStripHitPos = strip_pos_layerA;
1243 Double_t jStripHitPos = strip_pos_layerB;
1248 if( Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01 ) {
1250 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1251 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1253 xcoord = (jB - iB)/(Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
1254 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1258 if( Abs(iAngleRad) < 0.01 ) {
1260 if(StripLayers[ilayer].GetStripNumberingOrder() ==
LeftToRight ) {
1261 xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1264 xcoord = StripLayers[ilayer].GetXRightStripBorderPoint() - (iStripHitPos-StripLayers[ilayer].GetFirstStripNumber())*StripLayers[ilayer].GetPitch();
1267 Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
1268 ycoord = Tan(j_PiOver2MinusAngleRad)*xcoord + jB;
1271 if( Abs(jAngleRad) < 0.01 ) {
1273 if(StripLayers[jlayer].GetStripNumberingOrder() ==
LeftToRight ) {
1274 xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1277 xcoord = StripLayers[jlayer].GetXRightStripBorderPoint() - (jStripHitPos-StripLayers[jlayer].GetFirstStripNumber())*StripLayers[jlayer].GetPitch();
1280 Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
1281 ycoord = Tan(i_PiOver2MinusAngleRad)*xcoord + iB;
1292 if( StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord) && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord) ) {