ATLAS Offline Software
MdtVsTgcRawData_MidstationMatching.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 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 // Use Midstation Segments to look for additional Midstation Segments
32 // Only the Midstation covers some trajectories so normal SegmTracks cannot be constructed
33 void
34 MdtVsTgcRawDataValAlg::MidstationOnlyCheck(std::vector<const Muon::MuonSegment*> (&sortedSegments)[2][4],
35  std::vector<const Muon::MuonSegment*> (&disqualifiedSegments)[2][4],
36  const Muon::TgcPrepDataContainer *tgc_prepcontainer){
37  // Define Cuts:
38  // Cut for matching Segments
39  const float dRhoCutSegmentMatching = 1000;
40  const float dPhiCutSegmentMatching = M_PI/8;
41 
42  // Cuts for number of Mdt Measurements required for Midstation only track
43  const int nMeasCutMdtMidstation = 5;
44  const int nMeasCutTGCMidPRD[2] = {2,2};
45 
46  // Loose cut for PRD and Extrapolated
47  const float dPhiCut_Loose = M_PI/8;
48  // Cut for Global Position Efficiencies
49  const float dPhiCutGlobal[2] = {static_cast<float>(M_PI/24),static_cast<float>(M_PI/12)};//[WireStrip]
50  const float dRhoCutGlobal[2] = { 0.08, 0.5};//[WireStrip]
51  // Cut for Sector Efficiencies
52  const float dPhiCutSector[2] = { 0.2, 0.1};//[WireStrip]
53  const float dRhoCutSector[2] = { 300, 3000};//[WireStrip]
54  // Cut for TgcPrepData comparison
55  const float dPhiCutTPD[2] = { 0.15, 0.02};//[WireStrip]
56  const float dRhoCutTPD[2] = { 150, 3000};//[WireStrip]
57 
58  // Loop over sides
59  for(int i=0;i<2;i++){// AC
60  // Number of Segments found which passed all cuts
61  int nValidatedSegm = 0;
62 
63  // Following "Fill" variables are only used if one segment passed all cuts
64  // Flags for whether different stations can be checked by the segment
65  //bool canCheckGlobalFill[4] = {0, 0, 0, 0};
66  bool canCheckSectorFill[4] = {0, 0, 0, 0};
67  // Segment Global Position variables
68  // float posThetaFill = 0;
69  // float posPhiFill = 0;
70  // Segment Sector Position
71  int TGCstation_StationFEFill[4] = {-1,-1,-1,-1};// [TGCStation]
72  int TGCstation_StationEtaFill[4] = { 0, 0, 0, 0};// [TGCStation]
73  int TGCstation_StationPhiFill[4] = { 0, 0, 0, 0};// [TGCStation]
74  // Hit registered arrays
75  // bool hitregisteredFill[9][2] = {{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}};
76  bool sectorhitregisteredFill[9][2] = {{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}};
77 
78  bool skipSegm;int nDisqualifiedSegm; // used when checking the disqualified list for a segment
79 
80  // Make copy of disqualifiedSegment vector array
81  std::vector<const Muon::MuonSegment*> copyDisqualifiedSegments;
82  nDisqualifiedSegm=disqualifiedSegments[i][2].size();
83  copyDisqualifiedSegments.reserve(nDisqualifiedSegm);
84 for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)copyDisqualifiedSegments.push_back(disqualifiedSegments[i][2].at(ndis));
85 
87  // Apply preliminary nMdtMeas Cut
88  int nSegm = sortedSegments[i][2].size();
89  for(int n0=0; n0<nSegm;n0++){
90  // Get segment
91  const Muon::MuonSegment *segm0 = sortedSegments[i][2].at(n0);
92 
93  // Check copyDisqualifiedSegments for current segment
94  skipSegm=false;
95  nDisqualifiedSegm=copyDisqualifiedSegments.size();
96  for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)if(segm0==copyDisqualifiedSegments.at(ndis))skipSegm=true;
97  if(skipSegm)continue;
98 
99  // Apply nMDTMeasurements Cut to Segment Collection
100  int stationName = 0;
101  int nMdtMeas = 0;
102  // Loop through contained ROTs and identify used stations
103  for(unsigned int iROT=0; iROT<segm0->numberOfContainedROTs(); ++iROT) {
104  const Trk::RIO_OnTrack* rio = segm0->rioOnTrack(iROT);
105  if(!rio) continue;
106  Identifier id = rio->identify();
107  stationName = int(m_idHelperSvc->mdtIdHelper().stationName(id));
108 
109  if((stationName==17)||(stationName==18))nMdtMeas++;// MDT
110  }
111  // Cut Segments with insufficient numbers of hits in the stations
112  if(nMdtMeas<nMeasCutMdtMidstation){
113  copyDisqualifiedSegments.push_back(segm0);
114  }
115  }
116 
117  // Loop over Segments in MidStation
118  for(int n1=0; n1<nSegm;n1++){
120  // Start operation on Segment
121  // Get segment
122  const Muon::MuonSegment *segm1 = sortedSegments[i][2].at(n1);
123 
124  // Check copyDisqualifiedSegments for current segment
125  skipSegm=false;
126  nDisqualifiedSegm=copyDisqualifiedSegments.size();
127  for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)if(segm1==copyDisqualifiedSegments.at(ndis))skipSegm=true;
128  if(skipSegm)continue;
129 
130  // Flag for Grouping cut
131  bool failedGroupingCut = false;
132  // Flags for whether different stations can be checked by the segment
133  bool canCheckGlobal[4] = {true, true, true, false};
134  bool canCheckSector[4] = {true, true, true, false};
135 
136  // Get position variables
137  //const Trk::GlobalPosition segm1Pos = segm1->globalPosition();
138  const Amg::Vector3D segm1Pos = segm1->globalPosition();
139 
140  float segm1PosPhi = segm1Pos.phi();
141  float segm1PosThe = segm1Pos.theta();
142  float segm1PosZ = segm1Pos.z();
143  if(segm1PosPhi<0)segm1PosPhi+=2*M_PI;
144  if(segm1PosThe>M_PI/2) segm1PosThe=M_PI-segm1PosThe;
145  Amg::Vector3D segm1PosZunit(segm1Pos/std::abs(segm1PosZ));
146 
147 
149  // Apply nTGCStripMeasurements Cut to canCheck arrays
150  int stationName = 0;
151  int nTGCStrips[4] = { 0, 0, 0, 0};
152  // Loop through contained ROTs and identify used stations
153  for(unsigned int iROT=0; iROT<segm1->numberOfContainedROTs(); ++iROT){
154  const Trk::RIO_OnTrack* rio = segm1->rioOnTrack(iROT);
155  if(!rio) continue;
156  Identifier id = rio->identify();
157  stationName = int(m_idHelperSvc->mdtIdHelper().stationName(id));
158  bool isStrip = m_idHelperSvc->tgcIdHelper().isStrip(id);
159 
160  if(((stationName==41)||(stationName==42))&&isStrip)nTGCStrips[0]++;// TGC
161  if(((stationName==43)||(stationName==44))&&isStrip)nTGCStrips[1]++;// TGC
162  if(((stationName==45)||(stationName==46))&&isStrip)nTGCStrips[2]++;// TGC
163  if(((stationName==47)||(stationName==48))&&isStrip)nTGCStrips[3]++;// TGC
164  }
165 
166  // Don't check mid-stations when there are no strips in other mid-stations
167  if((nTGCStrips[1]==0)&&(nTGCStrips[2]==0)){canCheckSector[0]=false;canCheckGlobal[0]=false;}
168  if((nTGCStrips[0]==0)&&(nTGCStrips[2]==0)){canCheckSector[1]=false;canCheckGlobal[1]=false;}
169  if((nTGCStrips[0]==0)&&(nTGCStrips[1]==0)){canCheckSector[2]=false;canCheckGlobal[2]=false;}
170 
171 
173  // Do Grouping Cut
174  // Ignore any segments which are too close to any other segments not in the original disqualified list
175  for(int n2=0; n2<nSegm;n2++){
176  if(n1==n2)continue; // don't compare segment to itself
177  // Get segment
178  const Muon::MuonSegment *segm2 = sortedSegments[i][2].at(n2);
179 
180  // Check disqualifiedSegments for current segment
181  skipSegm=false;
182  nDisqualifiedSegm=disqualifiedSegments[i][2].size();
183  for(int ndis=0;ndis<nDisqualifiedSegm;ndis++)if(segm1==disqualifiedSegments[i][2].at(ndis))skipSegm=true;
184  if(skipSegm)continue;
185 
186  // Get position variables
187  const Amg::Vector3D segm2Pos = segm2->globalPosition();
188  float segm2PosRho = std::abs(segm2Pos.perp());
189  float segm2PosPhi = segm2Pos.phi();
190  float segm2PosZ = segm2Pos.z();
191  if(segm2PosPhi<0)segm2PosPhi+=2*M_PI;
192  Amg::Vector3D segm2PosZunit(segm2Pos/std::abs(segm2PosZ));
193 
194  // Apply preliminary phi cut between segm1 and segm2 positions
195  float dPhi_Segm1_Segm2 = segm1PosPhi-segm2PosPhi;
196  if(dPhi_Segm1_Segm2<-M_PI)dPhi_Segm1_Segm2+=2*M_PI;
197  if(dPhi_Segm1_Segm2> M_PI)dPhi_Segm1_Segm2-=2*M_PI;
198  if(std::abs(dPhi_Segm1_Segm2)<dPhiCutSegmentMatching){
199  failedGroupingCut=true;
200  break;
201  }
202 
203  // Extrapolate segm1 position to segm2's Z position
204  float dZ = std::abs(segm2PosZ)-std::abs(segm1PosZ);
205  Amg::Vector3D extrPos(segm1Pos+(segm1PosZunit*dZ));
206  float extrPosRho = std::abs(extrPos.perp());
207  float extrPosThe = extrPos.theta();
208  float extrPosPhi = extrPos.phi();
209  if(extrPosThe>M_PI/2) extrPosThe=M_PI-extrPosThe;
210  if(extrPosPhi<0)extrPosPhi+=2*M_PI;
211 
212  // Get differences between extrapolated and segm2 positions
213  float dRho_Extr_Segm2 = extrPosRho-segm2PosRho;
214  float dPhi_Extr_Segm2 = extrPosPhi-segm2PosPhi;
215  if(dPhi_Extr_Segm2<-M_PI)dPhi_Extr_Segm2+=2*M_PI;
216  if(dPhi_Extr_Segm2> M_PI)dPhi_Extr_Segm2-=2*M_PI;
217 
218  // Cut segments (segm1) which are inside difference cuts
219  if((std::abs(dPhi_Extr_Segm2)<dPhiCutSegmentMatching)||
220  (std::abs(dRho_Extr_Segm2)<dRhoCutSegmentMatching)){
221  failedGroupingCut=true;
222  break;
223  }
224  }// nSegm2
225  if(failedGroupingCut)continue;
226 
228  // Check track against TGC Sectors to find which it passes
229  // Initialize variables for TRE array search
230  int TGCStationNames[8] ={41, 42, 43, 44, 45, 46, 47, 48};
231  int TGCstation_StationFE[4] ={-1,-1,-1,-1};// [TGCStation]
232  int TGCstation_StationEta[4]={ 0, 0, 0, 0};// [TGCStation]
233  int TGCstation_StationPhi[4]={ 0, 0, 0, 0};// [TGCStation]
234  int nStationMatch[4] ={ 0, 0, 0, 0};// [TGCStation]
235 
236  // Loop through TRE array, finding sectors which match the track in each layer
237  for(int stationnameindex=0; stationnameindex<6; stationnameindex++){// Station {T1F,T1E,T2F,T2E,T3F,T3E}
238  // Skip stations which don't have sufficient Segments to run efficiency check
239  int stationName = TGCStationNames[stationnameindex];
240  int stationIndex= TGCstationname2stationindex(stationName);
241 
242  // Loop over StationEta&StationPhi
243  for(int stationeta=1; stationeta<=8; stationeta++){// AbsStationEta
244  for(int stationphi=1; stationphi<=48; stationphi++){// StationPhi
245  // Cut Station EtaPhi combinations with no TGC element
246  if(m_TREarray[stationnameindex][i][stationeta][stationphi]==nullptr)continue;
247  const MuonGM::TgcReadoutElement *tre=m_TREarray[stationnameindex][i][stationeta][stationphi];
248 
249  // Extrapolate position from nearest Station's Segment to Sector's Z
250  float sectorZ=tre->globalPosition().z();
251  float dZ_sector=std::abs(sectorZ)-std::abs(segm1PosZ);
252  //Trk::GlobalPosition sectorExtrapolatedPos = segm1Pos+(segm1PosZunit*dZ_sector);
253  Amg::Vector3D sectorExtrapolatedPos = segm1Pos+(segm1PosZunit*dZ_sector);
254 
255  // Convert extrapolated position to local position on the Sector
256  Identifier sector_id=tre->identify();
257  //const HepGeom::Point3D<double> sectorLocalPos3D=tre->globalToLocalCoords(sectorExtrapolatedPos, sector_id);
258  const Amg::Vector3D sectorLocalPos3D=tre->globalToLocalTransf(sector_id) * sectorExtrapolatedPos;
259  //Trk::LocalPosition sectorLocalPos2D(sectorLocalPos3D.y(),sectorLocalPos3D.z());
260  Amg::Vector2D sectorLocalPos2D(sectorLocalPos3D.y(),sectorLocalPos3D.z());
261 
262  double avWidth = (tre->getLongSsize()+tre->getSsize())/2;
263  //double dWidth = (tre->longWidth()-tre->shortWidth());
264  double length = tre->length();
265 
266  // Cut sectors which the track does not go through the central 80% of (avoids double counts)
267  double tol1=-0.1*(avWidth/2);
268  double tol2=-0.1*(length/2);
269 
270  bool insideSectorBounds=tre->bounds().inside(sectorLocalPos2D,tol1,tol2);
271  if(!insideSectorBounds)continue;
272  // Assign values to matching station variables
273  if(stationIndex<0) continue;
274  TGCstation_StationFE[stationIndex]= (tre->isForward()==false);
275  TGCstation_StationEta[stationIndex]=stationeta;
276  TGCstation_StationPhi[stationIndex]=stationphi;
277  nStationMatch[stationIndex]++;
278  }// StationPhi
279  }// StationEta
280  }// StationName
281 
282  // Don't check stations that don't have exactly 1 match for this track
283  for(int jTGC=0;jTGC<4;jTGC++){// TGC Station
284  if(nStationMatch[jTGC]==0){
285  canCheckSector[jTGC]=false;
286  }
287  else if(nStationMatch[jTGC]>1){
288  canCheckSector[jTGC]=false;
289  // Should be impossible, but happens due a problem with the bounds().inside function
290  // commenting out until this bug is resolved
291  //m_log << MSG::WARNING << "MidstationOnly: Number of matches for TGC" << jTGC+1 << " is " << nStationMatch[jTGC] << endl;
292  }
293  }// TGC Station
294 
295  // Cut Segment if no stations can be checked
296  if((!canCheckGlobal[0])&&(!canCheckGlobal[1])&&(!canCheckGlobal[2])&&(!canCheckGlobal[3]))continue;
298  // Check which PRD matches Segm1
299  // Initialise hit registered arrays
300  bool sectorhitregistered[9][2] = {{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}};
301  std::vector<const Muon::TgcPrepData*> tpdVector[2];
302 
303  // Loop over TGC Prep Data container
304  Muon::TgcPrepDataContainer::const_iterator prepit_end=tgc_prepcontainer->end();
305  for( Muon::TgcPrepDataContainer::const_iterator prepit=tgc_prepcontainer->begin();
306  prepit!=prepit_end;
307  ++prepit){
308 
309  //loop over TGC Prep Data collection
310  Muon::TgcPrepDataCollection::const_iterator prepitc_end=(*prepit)->end();
311  for( Muon::TgcPrepDataCollection::const_iterator prepitc=(*prepit)->begin();
312  prepitc!= prepitc_end;
313  ++prepitc){
314  // Get PRD and variables
315  const Muon::TgcPrepData* tpd=*prepitc;
316  const MuonGM::TgcReadoutElement *tre = tpd->detectorElement();
317  const std::string stationType = tre->getStationType();
318 
319  // Get id values
320  Identifier tgcid=(*prepitc)->identify();
321  int tgcAC=(tre->sideA()==false);//isNotAside a:0, c:1
322  int tgcFE=(tre->isForward()==false);//isNotForward f:0, e:1
323  int tgcWS=(m_idHelperSvc->tgcIdHelper().isStrip(tgcid));//isStrip w=0, s=1
324  int stationName = m_idHelperSvc->tgcIdHelper().stationName(tgcid);
325  int stationEta = std::abs(tre->getStationEta());
326  int stationPhi = tre->getStationPhi();
327  int gasGap = m_idHelperSvc->tgcIdHelper().gasGap(tgcid);
328 
329  // Cut hits except those from same side Midstation
330  if(tgcAC!=i) continue;
331  if(stationName>46 || stationName<41) continue;
332 
333  // Get layer number and stationIndex
334  int layer = TGCgetlayer(stationName,gasGap);
335  int stationIndex = TGCstationname2stationindex(stationName);
336  if(stationIndex==3)continue;
337 
338  // Get position variables
339  const Amg::Vector3D prdPos = tpd->globalPosition();
340  float tgcRho = std::abs(prdPos.perp());
341  float tgcPhi = prdPos.phi();
342  float tgcZ = prdPos.z();
343  if(tgcPhi<0)tgcPhi+=2*M_PI;
344 
345  // Extrapolate Segm1 to PRD Z position
346  float dZ = std::abs(tgcZ) - std::abs(segm1PosZ);
347  Amg::Vector3D tgcExtrapolatedPos = ((segm1Pos)+((segm1PosZunit)*dZ));
348 
349  // Get extrapolated variables
350  float tgcExtrRho = std::abs(tgcExtrapolatedPos.perp());
351  float tgcExtrPhi = tgcExtrapolatedPos.phi();
352  if(tgcExtrPhi<0)tgcExtrPhi+=2*M_PI;
353 
354  // Get differences between extrapolated and segm2 positions
355  float dRho = tgcRho-tgcExtrRho;
356  float dPhi = tgcPhi-tgcExtrPhi;
357  if(dPhi<-M_PI)dPhi+=2*M_PI;
358  if(dPhi> M_PI)dPhi-=2*M_PI;
359 
360  // Pass through loose phi cut to eliminate some noise
361  if(std::abs(dPhi)<dPhiCut_Loose){
362  // Fill PRD sagitta histograms
363  if(m_mvt_extrprdsag2[i][stationIndex][tgcFE][tgcWS][0]) m_mvt_extrprdsag2[i][stationIndex][tgcFE][tgcWS][0]->Fill(dRho);
364  if(m_mvt_extrprdsag2[i][stationIndex][tgcFE][tgcWS][2]) m_mvt_extrprdsag2[i][stationIndex][tgcFE][tgcWS][2]->Fill(dPhi);
365 
366  // Do Global check
367  if(canCheckGlobal[stationIndex]){
368  float dRhoCut = dRhoCutGlobal[tgcWS]*tgcExtrRho;
369  if(std::abs(dPhi)<dPhiCutGlobal[tgcWS] && std::abs(dRho)<dRhoCut){
370  }
371  }// global
372 
373  // Add PRD which matches Segm1 position to vector for further analysis
374  if(std::abs(dPhi)<dPhiCutSector[tgcWS] && std::abs(dRho)<dRhoCutSector[tgcWS]){
375  tpdVector[tgcWS].push_back(tpd);
376  }
377  // Do Sector Efficiency Check
378  if(canCheckSector[stationIndex]){
379  // Do check against PRD from matching Sectors
380  if((stationEta==TGCstation_StationEta[stationIndex])&&
381  (stationPhi==TGCstation_StationPhi[stationIndex])&&
382  (tgcFE==TGCstation_StationFE[stationIndex])){
383  // Do check
384  if(std::abs(dPhi)<dPhiCutSector[tgcWS] && std::abs(dRho)<dRhoCutSector[tgcWS]){
385  if(layer>=0)sectorhitregistered[layer][tgcWS]=true;// Sector hit
386  }
387  }
388  }// sector
389 
390  }// dPhi Loose Cut
391  }// TGC PRD Collection
392  }// TGC PRD Container
393 
394  // Cut Segment if no stations can be checked
395  if((!canCheckGlobal[0])&&(!canCheckGlobal[1])&&(!canCheckGlobal[2])&&(!canCheckGlobal[3]))continue;
397  // Do PRD checks
398  // Find vector of PRD which forms a coherent line in the vicinity of Segm1
399 
400  // Variables to hold best PRD matching results
401  std::vector<const Muon::TgcPrepData*> *bestTPDmatches[2];
402  bestTPDmatches[0] = nullptr;
403  bestTPDmatches[1] = nullptr;
404  if(bestTPDmatches[0]->size()>0) bestTPDmatches[0]->clear();
405  if(bestTPDmatches[1]->size()>0) bestTPDmatches[1]->clear();
406  int bestTPDlayerMatches[2][9] = {{0,0,0,0,0,0,0,0,0},
407  {0,0,0,0,0,0,0,0,0}};
408 
409  for(int k=0;k<2;k++){// WireStrip
410  // Variables to record quality of best match found
411  int nPRDMax = 0;
412  int nlayerMax = 0;
413 
414  // Loop over PRD vector
415  int nTPD = tpdVector[k].size();
416  for(int iTPD1=0;iTPD1<nTPD;iTPD1++){
417  // Variables to hold matches found for this PRD
418  std::vector<const Muon::TgcPrepData*> thisTPDmatches;
419  int thisTPDlayerMatches[9] = {0,0,0,0,0,0,0,0,0};
420 
421  // Get position variables
422  const Amg::Vector3D prdPos1 = tpdVector[k].at(iTPD1)->globalPosition();
423 
424  float prd1Phi = prdPos1.phi();
425  float prd1Z = prdPos1.z();
426  if(prd1Phi<0)prd1Phi+=2*M_PI;
427  Amg::Vector3D prd1PosZunit(prdPos1/std::abs(prd1Z));
428 
429  // Get id values
430  Identifier tgcid1=(tpdVector[k].at(iTPD1))->identify();
431  int stationName1 = m_idHelperSvc->tgcIdHelper().stationName(tgcid1);
432  int gasGap1 = m_idHelperSvc->tgcIdHelper().gasGap(tgcid1);
433  int layer1 = TGCgetlayer(stationName1,gasGap1);
434  if(layer1>=0)thisTPDlayerMatches[layer1]++;
435 
436  // Loop over PRD
437  for(int iTPD2=0;iTPD2<nTPD;iTPD2++){
438  if(iTPD2==iTPD1)continue;// do not check PRD against itself
439 
440  // Get position variables
441  const Amg::Vector3D prdPos2 = tpdVector[k].at(iTPD2)->globalPosition();
442  float prd2Rho = std::abs(prdPos2.perp());
443  float prd2Phi = prdPos2.phi();
444  float prd2Z = prdPos2.z();
445  if(prd2Phi<0)prd2Phi+=2*M_PI;
446 
447  // Extrapolate PRD1 to PRD2 Z position
448  float dZ = std::abs(prd2Z)- std::abs(prd1Z);
449  Amg::Vector3D prdExtrPos = ((prdPos1)+((prd1PosZunit)*dZ));
450 
451  // Get extrapolated variables
452  float prdExtrRho = std::abs(prdExtrPos.perp());
453  float prdExtrPhi = prdExtrPos.phi();
454  if(prdExtrPhi<0)prdExtrPhi+=2*M_PI;
455 
456  // Get differences between extrapolated and PRD2 positions
457  float dRho = prd2Rho-prdExtrRho;
458  float dPhi = prd2Phi-prdExtrPhi;
459  if(dPhi<-M_PI)dPhi+=2*M_PI;
460  if(dPhi> M_PI)dPhi-=2*M_PI;
461 
462  // Fill PRD comparison sagitta histograms
463  if(m_tgc_prdcompsag[i][k][0]) m_tgc_prdcompsag[i][k][0]->Fill(dRho);
464  if(m_tgc_prdcompsag[i][k][2]) m_tgc_prdcompsag[i][k][2]->Fill(dPhi);
465 
466  // Do check
467  if(std::abs(dPhi)<dPhiCutTPD[k] && std::abs(dRho)<dRhoCutTPD[k]){
468  // Get id values
469  Identifier tgcid2=(tpdVector[k].at(iTPD2))->identify();
470  int stationName2 = m_idHelperSvc->tgcIdHelper().stationName(tgcid2);
471  int gasGap2 = m_idHelperSvc->tgcIdHelper().gasGap(tgcid2);
472  int layer2 = TGCgetlayer(stationName2,gasGap2);
473 
474  // Add PRD2 to matches for PRD1
475  if(layer2>=0)thisTPDlayerMatches[layer2]++;
476  thisTPDmatches.push_back(tpdVector[k].at(iTPD2));
477  }
478  }// nTPD2
479 
480  // Get quality of matching array produced
481  int nPRDCurrent = 0;
482  int nlayerCurrent = 0;
483  for(int l=0;l<9;l++){
484  if(thisTPDlayerMatches[l]>0)nlayerCurrent++;
485  nPRDCurrent+=thisTPDlayerMatches[l];
486  }
487 
488  // Check quality variables against maximum found
489  if(nlayerMax <= nlayerCurrent){
490  if(nPRDMax < nPRDCurrent){
491  // Set maximum values to current segment's values
492  nlayerMax = nlayerCurrent;
493  nPRDMax = nPRDCurrent;
494  bestTPDmatches[k] = &thisTPDmatches;
495  for(int l=0;l<9;l++){
496  bestTPDlayerMatches[k][l] = thisTPDlayerMatches[l];
497  }
498  }
499  }
500  }// nTPD1
501 
502  // If matching array was somehow empty (should be impossible)
503  if(nlayerMax==0)continue;
504  if(bestTPDmatches[k]->size()==0){
505  ATH_MSG_WARNING( "MidstationOnly: empty bestTPDmatches["<<k<<"] passed" );
506  continue;
507  }
508 
509  // Set canCheck variables based on contents of matched PRD array
510  for(int jTGC1=0;jTGC1<3;jTGC1++){// TGC Stations
511  // Number of matched PRD found in other Stations
512  int nMatchOther = 0;
513  for(int l=0;l<9;l++){
514  int jTGC2 = TGClayer2stationindex(l);
515  if(jTGC1==jTGC2)continue;
516  nMatchOther+=bestTPDlayerMatches[k][l];
517  }
518 
519  // Check against cut
520  if(nMatchOther<nMeasCutTGCMidPRD[k]){
521  canCheckGlobal[jTGC1]=false;
522  canCheckSector[jTGC1]=false;
523  }
524  }// TGC Stations
525  }// WireStrip
526 
527  // Cut Segment if no stations can be checked
528  if((!canCheckGlobal[0])&&(!canCheckGlobal[1])&&(!canCheckGlobal[2])&&(!canCheckGlobal[3]))continue;
530  // Segment has passed all global cuts, fill histogram variables
531  // If no segments have already passed all cuts
532  if(nValidatedSegm==0){
533  for(int jTGC=0;jTGC<4;jTGC++){// TGC Stations
534  // Assign values of variables to fill histograms
535  // posThetaFill = segm1PosThe;
536  // posPhiFill = segm1PosPhi;
537 
538  //canCheckGlobalFill[jTGC] = canCheckGlobal[jTGC];
539  canCheckSectorFill[jTGC] = canCheckSector[jTGC];
540 
541  TGCstation_StationFEFill[jTGC] = TGCstation_StationFE[jTGC];
542  TGCstation_StationEtaFill[jTGC] = TGCstation_StationEta[jTGC];
543  TGCstation_StationPhiFill[jTGC] = TGCstation_StationPhi[jTGC];
544  }
545  for(int l=0;l<9;l++){
546  for(int k=0;k<2;k++){
547  sectorhitregisteredFill[l][k] = sectorhitregistered[l][k];
548  // hitregisteredFill[l][k] = hitregistered[l][k];
549  }
550  }
551  }// TGC Stations
552  nValidatedSegm++;
553  }// nSegm1
554 
555  // If only one Segment was validated on this side
556  if(nValidatedSegm==1){
557  for(int l=0;l<9;l++){//Layer
558  int stationIndex = TGClayer2stationindex(l);
559  for(int k=0;k<2;k++){// WireStrip
560  // If this station can be checked
561  if(canCheckSectorFill[stationIndex]){
562  if((TGCstation_StationFEFill[stationIndex]<0)||(TGCstation_StationEtaFill[stationIndex]==0)||(TGCstation_StationPhiFill[stationIndex]==0)){
563  ATH_MSG_WARNING( "MidstationOnly: canCheckSector passed for jTGC=" << stationIndex
564  << " but, FE="<<TGCstation_StationFEFill[stationIndex]
565  << " Eta="<<TGCstation_StationEtaFill[stationIndex]
566  << " Phi=" << TGCstation_StationPhiFill[stationIndex] );
567  continue;
568  }
569  // Get Sector histogram indexes
570  int stationMap_EtaIndex=getStationMapIndex(1, l, TGCstation_StationFEFill[stationIndex], TGCstation_StationEtaFill[stationIndex], TGCstation_StationPhiFill[stationIndex]);
571  int stationMap_PhiIndex=getStationMapIndex(2, l, TGCstation_StationFEFill[stationIndex], TGCstation_StationEtaFill[stationIndex], TGCstation_StationPhiFill[stationIndex]);
572 
573  // Fill Sector efficiency histograms
574  if(sectorhitregisteredFill[l][k]){// Hit in Sector matches extrapolated track
575  m_eff_stationmapmid[i][k][1]->Fill(stationMap_EtaIndex, stationMap_PhiIndex);
576  }
577  m_eff_stationmapmid[i][k][2]->Fill(stationMap_EtaIndex, stationMap_PhiIndex);
578  }
579  }// WireStrip
580  }// layer
581  }
582 
583  }// AC
584 
585 }// End of Function
TH2::Fill
int Fill(double, double)
Definition: rootspy.cxx:382
Muon::nsw::STGTPSegments::moduleIDBits::stationPhi
constexpr uint8_t stationPhi
station Phi 1 to 8
Definition: NSWSTGTPDecodeBitmaps.h:129
Muon::MuonPrepDataContainer
Template for Muon PRD containers (which are basically collections of MuonPrepDataCollections).
Definition: MuonPrepDataContainer.h:42
dumpTgcDigiDeadChambers.gasGap
list gasGap
Definition: dumpTgcDigiDeadChambers.py:33
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
Trk::SurfaceBounds::inside
virtual bool inside(const Amg::Vector2D &locpo, double tol1=0., double tol2=0.) const =0
Each Bounds has a method inside, which checks if a LocalPosition is inside the bounds.
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
MdtVsTgcRawDataValAlg::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MdtVsTgcRawDataValAlg.h:69
dumpTgcDigiDeadChambers.stationName
dictionary stationName
Definition: dumpTgcDigiDeadChambers.py:30
MuonGM::TgcReadoutElement::isForward
bool isForward() const
Returns true if the chamber is mounted on the most inner ring, i.e. a TxF chamber.
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Muon::MuonSegment::rioOnTrack
const Trk::RIO_OnTrack * rioOnTrack(unsigned int) const
returns the RIO_OnTrack (also known as ROT) objects depending on the integer
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:187
EventPrimitivesHelpers.h
MuonContainer.h
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Muon::TgcPrepData::detectorElement
virtual const MuonGM::TgcReadoutElement * detectorElement() const override final
Returns the detector element corresponding to this PRD The pointer will be zero if the det el is not ...
Definition: TgcPrepData.h:120
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
Trk::RIO_OnTrack
Definition: RIO_OnTrack.h:70
MuonGM::MuonReadoutElement::getLongSsize
double getLongSsize() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:199
MdtVsTgcRawDataValAlg::m_tgc_prdcompsag
TH1 * m_tgc_prdcompsag[2][2][4]
Definition: MdtVsTgcRawDataValAlg.h:149
Muon::MuonSegment::numberOfContainedROTs
unsigned int numberOfContainedROTs() const
number of RIO_OnTracks
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:199
MdtVsTgcRawDataValAlg::TGClayer2stationindex
int TGClayer2stationindex(int l)
Definition: MdtVsTgcRawData_TGCEffCheck.cxx:234
MuonGM::MuonReadoutElement::getSsize
double getSsize() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:196
MdtVsTgcRawDataValAlg::m_TREarray
const MuonGM::TgcReadoutElement * m_TREarray[8][2][9][49]
Definition: MdtVsTgcRawDataValAlg.h:94
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
GeoPrimitives.h
lumiFormat.i
int i
Definition: lumiFormat.py:92
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:530
MdtVsTgcRawDataValAlg::m_mvt_extrprdsag2
TH1 * m_mvt_extrprdsag2[2][4][2][2][4]
Definition: MdtVsTgcRawDataValAlg.h:148
MuonGM::TgcReadoutElement
A TgcReadoutElement corresponds to a single TGC chamber; therefore typically a TGC station contains s...
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/TgcReadoutElement.h:42
IdentifiableContainerMT::end
const_iterator end() const
return const_iterator for end of container
Definition: IdentifiableContainerMT.h:242
TrkEventPrimitivesDict.h
IdentifiableContainerMT::const_iterator
Definition: IdentifiableContainerMT.h:82
IdentifiableContainerMT::begin
const_iterator begin() const
return const_iterator for first entry
Definition: IdentifiableContainerMT.h:236
MuonGM::MuonReadoutElement::globalPosition
const Amg::Vector3D globalPosition() const
Definition: MuonDetDescr/MuonReadoutGeometry/src/MuonReadoutElement.cxx:47
TH1::Fill
int Fill(double)
Definition: rootspy.cxx:285
MdtVsTgcRawDataValAlg::MidstationOnlyCheck
void MidstationOnlyCheck(std::vector< const Muon::MuonSegment * >(&sortedSegments)[2][4], std::vector< const Muon::MuonSegment * >(&disqualifiedSegments)[2][4], const Muon::TgcPrepDataContainer *tgc_prepcontainer)
Definition: MdtVsTgcRawData_MidstationMatching.cxx:34
EventPrimitives.h
dumpTgcDigiThreshold.isStrip
list isStrip
Definition: dumpTgcDigiThreshold.py:33
RIO_OnTrack.h
MdtVsTgcRawDataValAlg::getStationMapIndex
int getStationMapIndex(int x, int l, int stationFE, int stationEta, int stationPhi)
Definition: MdtVsTgcRawData_TGCEffCheck.cxx:259
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Rtt_histogram.n1
n1
Definition: Rtt_histogram.py:21
MdtVsTgcRawDataValAlg::TGCstationname2stationindex
int TGCstationname2stationindex(int stationName)
Definition: MdtVsTgcRawData_TGCEffCheck.cxx:246
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GeoPrimitivesHelpers.h
MdtVsTgcRawDataValAlg::TGCgetlayer
int TGCgetlayer(int stationName, int g)
Definition: MdtVsTgcRawData_TGCEffCheck.cxx:196
MuonGM::MuonReadoutElement::getStationType
std::string getStationType() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:189
MuonGM::MuonReadoutElement::identify
Identifier identify() const override final
Returns the ATLAS Identifier of the MuonReadOutElement.
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:184
Muon::TgcPrepData
Class to represent TGC measurements.
Definition: TgcPrepData.h:32
MuonGM::TgcReadoutElement::length
double length() const
Definition: MuonDetDescr/MuonReadoutGeometry/src/TgcReadoutElement.cxx:53
Trk::RIO_OnTrack::identify
virtual Identifier identify() const final
return the identifier -extends MeasurementBase
Definition: RIO_OnTrack.h:155
Muon::TgcPrepData::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
Returns the global position.
Definition: TgcPrepData.h:125
Muon::MuonSegment::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
global position
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:157
MuonGM::MuonReadoutElement::sideA
bool sideA() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:207
Muon::MuonSegment
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:45
Muon::nsw::STGTPSegments::moduleIDBits::stationEta
constexpr uint8_t stationEta
1 to 3
Definition: NSWSTGTPDecodeBitmaps.h:127
MuonGM::MuonClusterReadoutElement::bounds
virtual const Trk::SurfaceBounds & bounds() const override
Return the boundaries of the element.
Definition: MuonClusterReadoutElement.h:127
MuonGM::TgcReadoutElement::globalToLocalTransf
Amg::Transform3D globalToLocalTransf(const Identifier &id) const
Returns the global -> local transformation.
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
MdtVsTgcRawDataValAlg::m_eff_stationmapmid
TH2 * m_eff_stationmapmid[2][2][4]
Definition: MdtVsTgcRawDataValAlg.h:144
MuonGM::MuonReadoutElement::getStationPhi
int getStationPhi() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:194
MuonGM::MuonReadoutElement::getStationEta
int getStationEta() const
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonReadoutElement.h:193
fitman.k
k
Definition: fitman.py:528
Generated on Tue Jul 2 2024 21:21:26 for ATLAS Offline Software by   doxygen 1.8.18