ATLAS Offline Software
MdtVsTgcRawData_SegmMatching.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Package : MdtVsTgcRawDataValAlg
7 // Author: M.King(Kobe)
8 // Feb. 2011
9 //
10 // DESCRIPTION:
11 // Subject: TGC Efficiency -->TGC Efficiency plots including EIFI by comparing with MDT Segments
13 
14 #include "../MdtVsTgcRawDataValAlg.h"
15 
16 #include "TrkEventPrimitives/TrkEventPrimitivesDict.h"
17 #include "EventPrimitives/EventPrimitivesHelpers.h"
18 #include "EventPrimitives/EventPrimitives.h"
19 #include "GeoPrimitives/GeoPrimitives.h"
20 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
21 
22 #include "TrkRIO_OnTrack/RIO_OnTrack.h"
23 #include "muonEvent/MuonContainer.h"
24 
25 #include <inttypes.h>
26 
27 #include <sstream>
28 #include <algorithm>
29 #include <fstream>
30 
31 // Matches together segments from different stations into a track
32 void
33 MdtVsTgcRawDataValAlg::MatchMDTSegments(std::vector<const Muon::MuonSegment*> (&sortedSegments)[2][4],
34  std::vector<const Muon::MuonSegment*> (&disqualifiedSegments)[2][4],
35  std::vector<SegmTrack> (&matchedSegments)[2]){// Define Cuts:
36  // Cut for matching Segments
37  const float dRhoCutSegmentMatching = 1000;
38  const float dPhiCutSegmentMatching = M_PI/8;
39  const float dPhiCutSegmentDirectionChecking[4][4]={{ 0,M_PI/8, 0.1, 0.01},
40  {M_PI/8, 0,M_PI/8,M_PI/8},
41  { 0.06,M_PI/8, 0, 0.01},
42  { 0.05,M_PI/8, 0.05, 0}};
43  const float dTheCutSegmentDirectionChecking[4][4]={{ 0,M_PI/8, 0.05, 0.06},
44  {M_PI/8, 0,M_PI/8,M_PI/8},
45  { 0.04,M_PI/8, 0, 0.005},
46  { 0.02,M_PI/8, 0.005, 0}};
47  // Loop over sides
48  for(int i=0;i<2;i++){// AC
49  bool skipSegm;int nDisqualifiedSegm; // used when checking the disqualified list for a segment
50  //bool HasStationMatchSegm[4] = {false, false, false, false};// flags for whether the there are segments in each MDT station which were included in tracks
51  //bool HasMatchedTrack = false;// flag for whether a track was found on current side
52 
53  for(int jMDT1=3;jMDT1>=0;jMDT1--){// MDT Stations in reverse]
54  // Get number of segments
55  int nSegm1=sortedSegments[i][jMDT1].size();
56 
57  // Loop over Segments in Station
58  for(int n1=0; n1<nSegm1;n1++){
59  // Get segment
60  const Muon::MuonSegment *segm1 = sortedSegments[i][jMDT1].at(n1);
61 
62  // Check disqualifiedSegments for current segment
63  skipSegm=false;
64  nDisqualifiedSegm=disqualifiedSegments[i][jMDT1].size();
65  for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)if(segm1==disqualifiedSegments[i][jMDT1].at(ndis))skipSegm=true;
66  if(skipSegm)continue;
67 
68  // Get position variables
69  const Amg::Vector3D segm1Pos = segm1->globalPosition();
70 
71  float segm1PosPhi = segm1Pos.phi();
72  float segm1PosZ = segm1Pos.z();
73  if(segm1PosPhi<0)segm1PosPhi+=2*M_PI;
74  Amg::Vector3D segm1PosZunit(segm1Pos/std::abs(segm1PosZ));
75  Amg::Vector3D segm1Dir = segm1->globalDirection();
76 
77  float segm1DirThe = segm1Dir.theta();
78  float segm1DirPhi = segm1Dir.phi();
79  if(segm1DirThe>M_PI/2) segm1DirThe=M_PI-segm1DirThe;
80  if(segm1DirPhi<0) segm1DirPhi+=2*M_PI;
81 
82  // Initialise segment matching variables
83  bool stationMatchFound[4] = {false,false,false,false};
84  std::vector<const Muon::MuonSegment*> matchingSegments[4];
85  for(int jMDT2=0;jMDT2<4;jMDT2++)matchingSegments[jMDT2] = std::vector<const Muon::MuonSegment*>();
86  matchingSegments[jMDT1].push_back(segm1);
87  int nStationMatch=0;
88 
89  for(int jMDT2=3;jMDT2>=0;jMDT2--){// MDT Station in reverse
90  // Get number of segments
91  int nSegm2=sortedSegments[i][jMDT2].size();
92 
93  // Loop over Segments in Station
94  for(int n2=0; n2<nSegm2;n2++){
95  // Get segment
96  const Muon::MuonSegment *segm2 = sortedSegments[i][jMDT2].at(n2);
97  if(segm1==segm2)continue; //do not compare same segments
98 
99  // Check disqualifiedSegments for current segment
100  skipSegm = false;
101  nDisqualifiedSegm=disqualifiedSegments[i][jMDT2].size();
102  for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)if(segm1==disqualifiedSegments[i][jMDT2].at(ndis))skipSegm=true;
103  if(skipSegm)continue;
104 
106  // Position Cut
107  // Fill position variables
108  const Amg::Vector3D segm2Pos = segm2->globalPosition();
109 
110  float segm2PosRho = std::abs(segm2Pos.perp());
111  float segm2PosPhi = segm2Pos.phi();
112  float segm2PosThe = segm2Pos.theta();
113  float segm2PosZ = segm2Pos.z();
114  if(segm2PosThe>M_PI/2) segm2PosThe=M_PI-segm2PosThe;
115  if(segm2PosPhi<0)segm2PosPhi+=2*M_PI;
116  Amg::Vector3D segm2PosZunit(segm2Pos/std::abs(segm2PosZ));
117 
118  // Apply preliminary phi cut between segm1 and segm2 positions
119  float dPhi_Segm1_Segm2 = segm1PosPhi-segm2PosPhi;
120  if(dPhi_Segm1_Segm2<-M_PI)dPhi_Segm1_Segm2+=2*M_PI;
121  if(dPhi_Segm1_Segm2> M_PI)dPhi_Segm1_Segm2-=2*M_PI;
122  if(std::abs(dPhi_Segm1_Segm2)>dPhiCutSegmentMatching)continue;
123 
124  // Extrapolate segm1 position to segm2's Z position
125  float dZ = std::abs(segm2PosZ)-std::abs(segm1PosZ);
126  Amg::Vector3D extrPos(segm1Pos+(segm1PosZunit*dZ));
127 
128  float extrPosRho = std::abs(extrPos.perp());
129  float extrPosThe = extrPos.theta();
130  float extrPosPhi = extrPos.phi();
131  if(extrPosThe>M_PI/2) extrPosThe=M_PI-extrPosThe;
132  if(extrPosPhi<0)extrPosPhi+=2*M_PI;
133 
134  // Get differences between extrapolated and segm2 positions
135  float dRho_Extr_Segm2 = extrPosRho-segm2PosRho;
136  float dPhi_Extr_Segm2 = extrPosPhi-segm2PosPhi;
137  float dThe_Extr_Segm2 = extrPosThe-segm2PosThe;
138  if(dPhi_Extr_Segm2<-M_PI)dPhi_Extr_Segm2+=2*M_PI;
139  if(dPhi_Extr_Segm2> M_PI)dPhi_Extr_Segm2-=2*M_PI;
140 
141  // Fill difference histograms
142  if(m_mdt_segmmatchsag[i][jMDT1][jMDT2][0]) m_mdt_segmmatchsag[i][jMDT1][jMDT2][0]->Fill(dRho_Extr_Segm2);
143  if(m_mdt_segmmatchsag[i][jMDT1][jMDT2][2]) m_mdt_segmmatchsag[i][jMDT1][jMDT2][2]->Fill(dPhi_Extr_Segm2);
144  if(m_mdt_segmmatchsag[i][jMDT1][jMDT2][3]) m_mdt_segmmatchsag[i][jMDT1][jMDT2][3]->Fill(dThe_Extr_Segm2);
145 
146  // Cut segments (segm2) which are outside difference cuts
147  if(std::abs(dPhi_Extr_Segm2)>dPhiCutSegmentMatching)continue;
148  if(std::abs(dRho_Extr_Segm2)>dRhoCutSegmentMatching)continue;
149 
151  // Direction Cut
152  // Fill direction variables
153  //const Trk::GlobalDirection segm2Dir = segm2->globalDirection();
154  const Amg::Vector3D segm2Dir = segm2->globalDirection();
155 
156  float segm2DirThe = segm2Dir.theta();
157  float segm2DirPhi = segm2Dir.phi();
158  if(segm2DirThe>M_PI/2) segm2DirThe=M_PI-segm2DirThe;
159  if(segm2DirPhi<0) segm2DirPhi+=2*M_PI;
160 
161  // Vector Method
162  // Make vector between segment positions
163  //Trk::GlobalDirection segmVector;
164  Amg::Vector3D segmVector;
165 
166  if(jMDT2>jMDT1)segmVector = segm2Pos-segm1Pos;
167  else segmVector = segm1Pos-segm2Pos;
168  float segmVecThe = segmVector.theta();
169  float segmVecPhi = segmVector.phi();
170  if(segmVecThe>M_PI/2) segmVecThe=M_PI-segmVecThe;
171  if(segmVecPhi<0) segmVecPhi+=2*M_PI;
172 
173  // Get difference between vector direction and segment directions
174  float dThe_Vec_Segm1 = segmVecThe-segm1DirThe;
175  float dPhi_Vec_Segm1 = segmVecPhi-segm1DirPhi;
176  if(dPhi_Vec_Segm1<-M_PI)dPhi_Vec_Segm1+=2*M_PI;
177  if(dPhi_Vec_Segm1> M_PI)dPhi_Vec_Segm1-=2*M_PI;
178  float dThe_Vec_Segm2 = segmVecThe-segm2DirThe;
179  float dPhi_Vec_Segm2 = segmVecPhi-segm2DirPhi;
180  if(dPhi_Vec_Segm2<-M_PI)dPhi_Vec_Segm2+=2*M_PI;
181  if(dPhi_Vec_Segm2> M_PI)dPhi_Vec_Segm2-=2*M_PI;
182 
183  // Fill histograms
184  if(jMDT1>jMDT2){
185  if(m_mdt_trackchecksag[i][jMDT2][jMDT1][2][0]) m_mdt_trackchecksag[i][jMDT2][jMDT1][2][0]->Fill(dPhi_Vec_Segm2);
186  if(m_mdt_trackchecksag[i][jMDT2][jMDT1][3][0]) m_mdt_trackchecksag[i][jMDT2][jMDT1][3][0]->Fill(dThe_Vec_Segm2);
187  if(m_mdt_trackchecksag[i][jMDT2][jMDT1][2][1]) m_mdt_trackchecksag[i][jMDT2][jMDT1][2][1]->Fill(dPhi_Vec_Segm1);
188  if(m_mdt_trackchecksag[i][jMDT2][jMDT1][3][1]) m_mdt_trackchecksag[i][jMDT2][jMDT1][3][1]->Fill(dThe_Vec_Segm1);
189  }
190  else if(jMDT1<jMDT2){
191  if(m_mdt_trackchecksag[i][jMDT1][jMDT2][2][0]) m_mdt_trackchecksag[i][jMDT1][jMDT2][2][0]->Fill(dPhi_Vec_Segm1);
192  if(m_mdt_trackchecksag[i][jMDT1][jMDT2][3][0]) m_mdt_trackchecksag[i][jMDT1][jMDT2][3][0]->Fill(dThe_Vec_Segm1);
193  if(m_mdt_trackchecksag[i][jMDT1][jMDT2][2][1]) m_mdt_trackchecksag[i][jMDT1][jMDT2][2][1]->Fill(dPhi_Vec_Segm2);
194  if(m_mdt_trackchecksag[i][jMDT1][jMDT2][3][1]) m_mdt_trackchecksag[i][jMDT1][jMDT2][3][1]->Fill(dThe_Vec_Segm2);
195  }
196 
197  // Direction Comparison Method
198  float dTheDir_Segm1_Segm2 = segm1DirThe-segm2DirThe;
199  float dPhiDir_Segm1_Segm2 = segm1PosPhi-segm2PosPhi;
200  if(dPhiDir_Segm1_Segm2<-M_PI)dPhiDir_Segm1_Segm2+=2*M_PI;
201  if(dPhiDir_Segm1_Segm2> M_PI)dPhiDir_Segm1_Segm2-=2*M_PI;
202  if(m_mdt_trackdirdirsag[i][jMDT1][jMDT2][2]) m_mdt_trackdirdirsag[i][jMDT1][jMDT2][2]->Fill(dTheDir_Segm1_Segm2);
203  if(m_mdt_trackdirdirsag[i][jMDT1][jMDT2][3]) m_mdt_trackdirdirsag[i][jMDT1][jMDT2][3]->Fill(dPhiDir_Segm1_Segm2);
204  // Cut using Vector Method
205  if(dPhi_Vec_Segm1>dPhiCutSegmentDirectionChecking[jMDT1][jMDT2] ||
206  dThe_Vec_Segm1>dTheCutSegmentDirectionChecking[jMDT1][jMDT2] ||
207  dPhi_Vec_Segm2>dPhiCutSegmentDirectionChecking[jMDT2][jMDT1] ||
208  dThe_Vec_Segm2>dTheCutSegmentDirectionChecking[jMDT2][jMDT1]) continue;
209 
210  // Match is found, increment counter and assign segm2 to match found array
211  matchingSegments[jMDT2].push_back(segm2);
212  }// nSegms2
213 
214  if(matchingSegments[jMDT2].size()==1){
215  stationMatchFound[jMDT2]=true;
216  nStationMatch++;
217  }
218  }// Reverse MDT Stations2
219 
220  // If matches found add to matchedSegments and disqualify all segments in array
221  if(nStationMatch>1){
222  //HasMatchedTrack=true;
223  const Muon::MuonSegment *segmArray[4] = {nullptr,nullptr,nullptr,nullptr};
224  for(int jMDT2=0;jMDT2<4;jMDT2++){
225  if(stationMatchFound[jMDT2]){
226  segmArray[jMDT2]=matchingSegments[jMDT2].at(0);
227  //HasStationMatchSegm[jMDT2]=true;
228  disqualifiedSegments[i][jMDT2].push_back(matchingSegments[jMDT2].at(0));
229  }
230  }
231  SegmTrack newTrack(segmArray);
232  matchedSegments[i].push_back(newTrack);
233  }// If Matched Track Found
234 
235  }// nSegms1
236  }// Reverse MDT Stations1
237 
238  }// AC
239 
240  return;
241 }// End of function
MdtVsTgcRawDataValAlg::m_mdt_trackdirdirsag
TH1 * m_mdt_trackdirdirsag[2][4][4][4]
Definition: MdtVsTgcRawDataValAlg.h:152
EventPrimitivesHelpers.h
MuonContainer.h
M_PI
#define M_PI
Definition: ActiveFraction.h:11
MdtVsTgcRawDataValAlg::m_mdt_trackchecksag
TH1 * m_mdt_trackchecksag[2][4][4][4][2]
Definition: MdtVsTgcRawDataValAlg.h:153
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
GeoPrimitives.h
lumiFormat.i
int i
Definition: lumiFormat.py:92
MdtVsTgcRawDataValAlg::MatchMDTSegments
void MatchMDTSegments(std::vector< const Muon::MuonSegment * >(&sortedSegments)[2][4], std::vector< const Muon::MuonSegment * >(&disqualifiedSegments)[2][4], std::vector< SegmTrack >(&matchedSegments)[2])
Definition: MdtVsTgcRawData_SegmMatching.cxx:33
TrkEventPrimitivesDict.h
TH1::Fill
int Fill(double)
Definition: rootspy.cxx:285
EventPrimitives.h
RIO_OnTrack.h
SegmTrack
Definition: SegmTrack.h:22
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Rtt_histogram.n1
n1
Definition: Rtt_histogram.py:21
GeoPrimitivesHelpers.h
MdtVsTgcRawDataValAlg::m_mdt_segmmatchsag
TH1 * m_mdt_segmmatchsag[2][4][4][4]
Definition: MdtVsTgcRawDataValAlg.h:150
Muon::MuonSegment::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
global position
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:157
Muon::MuonSegment
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:45
Muon::MuonSegment::globalDirection
const Amg::Vector3D & globalDirection() const
global direction
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:163
Generated on Tue Jul 2 2024 21:21:27 for ATLAS Offline Software by   doxygen 1.8.18