ATLAS Offline Software
TileMuonFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /********************************************************************
6  *
7  * NAME: TileMuonFitter
8  * PACKAGE: offline/TileCalorimeter/TileCosmicAlgs
9  *
10  * AUTHOR : J. Maneira
11  * CREATED: 10-Jul-2006
12  * REVISED: 25-Jan-2008 Added Hough Transform Method by Luciano M. de Andrade
13  *
14  *
15  * 08-Apr-2009 Remove masking of bad cells (now done by
16  * TileCellBuilder)
17  * Optimized loop over CaloCellContainer to read only
18  Tile Cells
19  *
20  *
21  *
22  * PURPOSE: Fit a straight muon track through a selected subset
23  * of Tile Cells, calculate horizontal plane(or vertical for
24  * single beam) crossing time. Use track parameters to
25  * estimate energy and path length.
26  *
27  * Input: TileCellContainer
28  * Output: ComissioningEvent/ComTime object
29  * TileEvent/TileCosmicMuon object
30  * Parameters:
31  *
32  *
33  ********************************************************************/
34 // Tile includes
35 #include "TileCosmicAlgs/TileMuonFitter.h"
36 #include "TileIdentifier/TileHWID.h"
37 #include "TileDetDescr/TileDetDescrManager.h"
38 #include "TileEvent/TileCosmicMuon.h"
39 
40 // Calo includes
41 #include "CaloDetDescr/CaloDetDescrElement.h"
42 #include "CaloIdentifier/CaloID.h"
43 
44 // Athena includes
45 #include "EventContainers/SelectAllObject.h"
46 #include "StoreGate/ReadHandle.h"
47 #include "StoreGate/WriteHandle.h"
48 #include "Identifier/HWIdentifier.h"
49 
50 // Gaudi includes
51 #include "GaudiKernel/ISvcLocator.h"
52 
53 // External includes
54 #include "Minuit2/MnMigrad.h"
55 #include "Minuit2/MnUserParameterState.h"
56 #include "Minuit2/FunctionMinimum.h"
57 
58 // ROOT includes
59 #include "CLHEP/Vector/ThreeVector.h"
60 #include "CLHEP/Units/SystemOfUnits.h"
61 #include "CLHEP/Units/PhysicalConstants.h"
62 
63 #include <cmath>
64 #ifndef __APPLE__
65 #define pi M_PIl
66 #else
67 #define pi M_PI
68 #endif
69 
70 using CLHEP::Hep3Vector;
71 using CLHEP::c_light;
72 
73 // Hough Transform Constants --------------------------------------------------
74 
75 #define BIN_RES_RXY 100.0 // bin resolution for Hough parameter rxy (in milimeters)
76 #define BIN_RES_RZY 100.0 // bin resolution for Hough parameter rzy (in milimeters)
77 #define BIN_RES_AXY (2.0*pi/128.0) // bin resolution for Hough parameter axy (in radian)
78 #define BIN_RES_AZY (2.0*pi/128.0) // bin resolution for Hough parameter azy (in radian)
79 #define MAX_NCELLS 500U // max number of cells to compute HT
80 #define MIN_NCELLS 3U // min number of cells in track
81 #define MIN_ENERGY 1.0F // minimal track energy (in GeV)
82 #define MIN_ENERGY_MEV 1000.0F // minimal track energy (in MeV)
83 #define SHIFT_X 5000.0 // displacement in X
84 #define SHIFT_Z 8000.0 // displacement in Z
85 
86 // -----------------------------------------------------------------------------
87 
88 using namespace ROOT::Minuit2;
89 
90 const CaloCell_ID::SUBCALO TileMuonFitter::m_caloIndex = CaloCell_ID::TILE;
91 
92 //
93 // Constructor
94 //
95 TileMuonFitter::TileMuonFitter(const std::string& name, ISvcLocator* pSvcLocator)
96  : AthAlgorithm(name, pSvcLocator)
97  , m_tileID(nullptr)
98  , m_tileHWID(nullptr)
99  , m_tileMgr(nullptr)
100  , m_caloCells(nullptr)
101  , m_theTrack(nullptr)
102  , m_nCells(0)
103  , m_meanX(0.0)
104  , m_meanY(0.0)
105  , m_meanZ(0.0)
106  , m_weightedMeanX(0.0)
107  , m_weightedMeanY(0.0)
108  , m_weightedMeanZ(0.0)
109  , m_maxBottomIndex(0)
110  , m_maxTopIndex(0)
111  , m_reg1to2(false)
112 {
113  declareProperty("DoHoughTransform", m_doHoughTransform = true);
114  declareProperty("EThreshold", m_eThreshold = 250.);
115  declareProperty("DeltaTimeCut", m_deltaTimeCut = 6.);
116  declareProperty("MinimumCells", m_minimumCells = 2);
117  declareProperty("DoWeighted", m_doWeighted = true);
118  declareProperty("DoDensityWeighting", m_doDensity = true);
119  declareProperty("BeamType", m_beamType = "cosmics");
120 }
121 
122 // ********************************************************************
123 // ********************************************************************
124 // ********************************************************************
125 TileMuonFitter::~TileMuonFitter() {
126 }
127 
128 // ********************************************************************
129 // ********************************************************************
130 // ********************************************************************
131 StatusCode TileMuonFitter::initialize() {
132 
133  CHECK( detStore()->retrieve(m_tileMgr) );
134  CHECK( detStore()->retrieve(m_tileID) );
135  CHECK( detStore()->retrieve(m_tileHWID) );
136 
137  if (m_doHoughTransform) {
138  m_cosmicMuonContainerKey = "TileCosmicMuonHT";
139  }
140 
141  ATH_CHECK( m_cellContainerKey.initialize() );
142  ATH_CHECK( m_cosmicMuonContainerKey.initialize() );
143  ATH_CHECK( m_comTimeKey.initialize() );
144 
145  m_cellEnergy.resize(m_tileID->cell_hash_max());
146  m_cellWeight.resize(m_tileID->cell_hash_max());
147  m_cellTime.resize(m_tileID->cell_hash_max());
148  m_cellDeltaTime.resize(m_tileID->cell_hash_max());
149  m_cellHash.resize(m_tileID->cell_hash_max());
150  m_cellPosition.resize(m_tileID->cell_hash_max());
151 
152  if (m_minimumCells < 2) m_minimumCells = 2;
153 
154  ATH_MSG_INFO( "Loading: " << m_cellContainerKey.key() );
155  ATH_MSG_INFO( "Recording: " << m_cosmicMuonContainerKey.key() << " and " << m_comTimeKey.key() );
156  ATH_MSG_INFO( "Cell Energy Threshold: " << m_eThreshold << " (MeV)" );
157  ATH_MSG_INFO( "Delta Time Cell Cut: " << m_deltaTimeCut << " (ns)" );
158  ATH_MSG_INFO( "Minimum Number of cells for fit: " << m_minimumCells );
159 
160  if (m_doWeighted && m_doDensity)
161  ATH_MSG_INFO( "Weighting cell position and time with energy density." );
162 
163  if (m_doWeighted && !m_doDensity)
164  ATH_MSG_INFO( "Weighting cell position and time with energy." );
165 
166  if (!m_doWeighted)
167  ATH_MSG_INFO( "Not weighting cell position and time." );
168 
169  if (!m_beamType.compare("collisions")
170  || !m_beamType.compare("singlebeam")
171  || !m_beamType.compare("cosmics")) {
172 
173  ATH_MSG_INFO( "Setting up TMF for beam type = " << m_beamType );
174  } else {
175  ATH_MSG_WARNING( "BeamType " << m_beamType
176  << " from jobproperties not recognized. Setting to cosmics for TMF only." );
177 
178  m_beamType = "cosmics";
179  }
180 
181  m_tileDD_radiusLB.resize(4);
182  m_tileDD_radiusEB.resize(4);
183  m_tileDD_zEBA.resize(2);
184  m_tileDD_zEBC.resize(3);
185  m_tileDD_zLB.resize(2);
186  bool hack = false;
187 
188  std::vector<TileDetDescriptor*>::const_iterator first = m_tileMgr->tile_descriptors_begin();
189  std::vector<TileDetDescriptor*>::const_iterator last = m_tileMgr->tile_descriptors_end();
190  for (; first != last; ++first) {
191  const TileDetDescriptor * tileDD = (*first);
192  //tileDD->print();
193  if (tileDD->n_samp() != 3) break;
194  if (tileDD->n_eta(0) == 10 && tileDD->sign_eta() == 1) { //LBA
195  hack = (tileDD->dr(0) > 300);
196  m_tileDD_zLB[1] = tileDD->zcenter(0) + tileDD->dz(0) / 2.;
197  for (int is = 0; is < 3; is++)
198  m_tileDD_radiusLB[is] = tileDD->rcenter(is) - tileDD->dr(is) / 2.;
199  m_tileDD_radiusLB[3] = tileDD->rcenter(2) + tileDD->dr(2) / 2.;
200  } else if (tileDD->n_eta(0) == 10 && tileDD->sign_eta() == -1) { //LBC
201  m_tileDD_zLB[0] = tileDD->zcenter(0) - tileDD->dz(0) / 2.;
202  } else if (tileDD->n_eta(0) == 5 && tileDD->sign_eta() == 1) { //EBA
203  m_tileDD_zEBA[0] = tileDD->zcenter(0) - tileDD->dz(0) / 2.;
204  m_tileDD_zEBA[1] = tileDD->zcenter(0) + tileDD->dz(0) / 2.;
205  for (int is = 0; is < 3; is++)
206  m_tileDD_radiusEB[is] = tileDD->rcenter(is) - tileDD->dr(is) / 2.;
207  m_tileDD_radiusEB[3] = tileDD->rcenter(2) + tileDD->dr(2) / 2.;
208  } else if (tileDD->n_eta(0) == 5 && tileDD->sign_eta() == -1) { //EBC
209  m_tileDD_zEBC[0] = tileDD->zcenter(0) - tileDD->dz(0) / 2.;
210  m_tileDD_zEBC[1] = tileDD->zcenter(0) + tileDD->dz(0) / 2.;
211  }
212 
213  }
214 // WARNING: ugly hack! Remove when Geo-Model is updated
215 // this is needed to remove from the cell lengths the spacer
216 // plate after tilerow 11 and the front plate before tilerow 1
217 // see http://indico.cern.ch/conferenceDisplay.py?confId=59722
218  if (hack) {
219  ATH_MSG_INFO( "Old geometry detected: moving inner/outer radius by +10/-40 mm " );
220  m_tileDD_radiusLB[3] -= 40; // reduce outer radius by 40 mm
221  m_tileDD_radiusEB[3] -= 40; // reduce outer radius by 40 mm
222  m_tileDD_radiusLB[0] += 10; // increase inner radius by 10 mm
223  m_tileDD_radiusEB[0] += 10; // increase inner radius by 10 mm
224  }
225 // end of hack
226  ATH_MSG_INFO( "Geometry read from TileDetDescr for track length calculation:" );
227  for (int is = 0; is < 4; is++)
228  ATH_MSG_INFO( "Radius sampling " << is
229  << " LB: " << m_tileDD_radiusLB[is]
230  << " EB: " << m_tileDD_radiusEB[is] );
231 
232  ATH_MSG_INFO( "LB z between " << m_tileDD_zLB[0] << " and " << m_tileDD_zLB[1] );
233  ATH_MSG_INFO( "EBA z between " << m_tileDD_zEBA[0] << " and " << m_tileDD_zEBA[1] );
234  ATH_MSG_INFO( "EBC z between " << m_tileDD_zEBC[0] << " and " << m_tileDD_zEBC[1] );
235  ATH_MSG_INFO( "TileMuonFitter initialization completed" );
236 
237  return StatusCode::SUCCESS;
238 }
239 
240 // ********************************************************************
241 // ********************************************************************
242 // ********************************************************************
243 StatusCode TileMuonFitter::finalize() {
244 
245  ATH_MSG_INFO( "Finalizing" );
246 
247  return StatusCode::SUCCESS;
248 }
249 
250 // ********************************************************************
251 // ********************************************************************
252 // ********************************************************************
253 StatusCode TileMuonFitter::execute() {
254 
255  ATH_MSG_DEBUG( " start execute " );
256 
257  int fitStatus = 2; //not even try
258  setEventDefaults();
259 
260  SG::ReadHandle<CaloCellContainer> cellContainer(m_cellContainerKey);
261  if (!cellContainer.isValid()) {
262  ATH_MSG_WARNING( " Could not find container " << m_cellContainerKey.key() );
263  } else {
264  ATH_MSG_DEBUG( "Container " << m_cellContainerKey.key() << " with CaloCells found " );
265 
266  m_caloCells = cellContainer.cptr();
267 
268  // check on number of tile cells
269  if (m_caloCells->nCellsCalo(m_caloIndex) == 0) {
270  ATH_MSG_DEBUG( "no TileCells in CellContainer " << m_cellContainerKey.key() << "> for this event!" );
271  } else {
272  buildCells();
273  }
274 
275  if (eventSelection()) {
276  if (!m_doHoughTransform)
277  fitStatus = fitTrack();
278  else
279  fitStatus = houghTrack();
280  if (fitStatus) calculateTime();
281  }
282 
283  buildTileCosmicMuon(fitStatus);
284  buildComTime(fitStatus);
285  }
286 
287  // Execution completed.
288  ATH_MSG_DEBUG( "TileMuonFitter execution completed." );
289 
290  return StatusCode::SUCCESS;
291 }
292 
293 // ********************************************************************
294 // ********************************************************************
295 // ********************************************************************
296 
297 void TileMuonFitter::setEventDefaults() {
298  m_nCells = 0;
299  for (unsigned int i = 0; i < m_tileID->cell_hash_max(); i++) {
300  m_cellEnergy[i] = 0.0;
301  m_cellWeight[i] = 0.0;
302  m_cellTime[i] = 0.0;
303  m_cellDeltaTime[i] = 999.;
304  m_cellHash[i] = 0;
305  m_cellPosition[i].setX(0.0);
306  m_cellPosition[i].setY(0.0);
307  m_cellPosition[i].setZ(0.0);
308  }
309 
310  m_linePar.clear();
311  m_zeroCrossingTime.clear();
312  m_fitMinimum.clear();
313 }
314 
315 // ********************************************************************
316 // ********************************************************************
317 // ********************************************************************
318 
319 void TileMuonFitter::buildCells() {
320 
321  // CaloCellContainer::const_iterator collItr = celcoll->begin();
322  // CaloCellContainer::const_iterator lastColl = celcoll->end();
323 
324  size_t collItr = m_caloCells->indexFirstCellCalo(m_caloIndex);
325  size_t lastColl = m_caloCells->indexLastCellCalo(m_caloIndex);
326 
327  for (; collItr != lastColl; ++collItr) {
328 
329  const CaloCell* cell = (*m_caloCells)[collItr];
330 
331  const TileCell* tilecell = dynamic_cast<const TileCell*>(cell);
332  if (tilecell == 0) continue;
333 
334  Identifier cell_id = cell->ID();
335  IdentifierHash cell_idhash = m_tileID->cell_hash(cell_id);
336  CaloDetDescrElement *caloDDE = m_tileMgr->get_cell_element(cell_id);
337 
338  double ener = cell->energy();
339  double time = cell->time();
340  double deltatime = tilecell->timeDiff();
341 
342  if (ener < m_eThreshold) continue; // below threshold removal
343  if (fabs(caloDDE->eta()) < 0.05 && caloDDE->r() > 3500) {
344  ATH_MSG_DEBUG( "Cell eta: " << caloDDE->eta()
345  << " Cell r: " << caloDDE->r()
346  << " Tile diff: " << tilecell->timeDiff() );
347 
348  deltatime = 0;
349  }
350  if (fabs(deltatime) > m_deltaTimeCut) continue; // cut high time imbalance cells
351 
352  double volume = caloDDE->volume();
353  if (volume == 0.0) {
354  ATH_MSG_DEBUG( "Warning: Skipping cell with zero volume!" );
355  continue;
356  }
357 
358  if (m_doDensity)
359  m_cellWeight[m_nCells] = ener / volume;
360  else
361  m_cellWeight[m_nCells] = ener;
362  m_cellEnergy[m_nCells] = ener;
363  m_cellTime[m_nCells] = time;
364  m_cellDeltaTime[m_nCells] = deltatime;
365  m_cellHash[m_nCells] = cell_idhash;
366  m_cellPosition[m_nCells].setX(caloDDE->x());
367  m_cellPosition[m_nCells].setY(caloDDE->y());
368  m_cellPosition[m_nCells].setZ(caloDDE->z());
369  m_nCells++;
370 
371  } //end of CaloCellContainer cycle
372 }
373 
374 // ********************************************************************
375 // ********************************************************************
376 // ********************************************************************
377 
378 bool TileMuonFitter::eventSelection() {
379  // simple check if there are cells on top and bottom halfs
380  // returns 1 if both top and bottom have energy,
381  // JM, Aug08: Change to allow for single beam halo muon
382  // looks for A/C side instead of top/bottom
383 
384  if (m_nCells < m_minimumCells) return false;
385 
386  double maxReg1Energy;
387  double maxReg2Energy;
388  double cellPosition;
389  int maxReg1Index;
390  int maxReg2Index;
391  bool selection = false;
392  maxReg1Energy = maxReg2Energy = 0.;
393  maxReg1Index = maxReg2Index = -1;
394 
395  for (int i = 0; i < m_nCells; i++) {
396  if (!m_beamType.compare("singlebeam")) {
397  cellPosition = m_cellPosition[i].getZ();
398  } else if (!m_beamType.compare("collisions")) {
399  cellPosition = 1.;
400  } else {
401  cellPosition = m_cellPosition[i].getY();
402  }
403 
404  if (cellPosition > 0) {
405  if (m_cellEnergy[i] > maxReg1Energy) {
406  maxReg1Energy = m_cellEnergy[i];
407  maxReg1Index = i;
408  }
409  } else {
410  if (m_cellEnergy[i] > maxReg2Energy) {
411  maxReg2Energy = m_cellEnergy[i];
412  maxReg2Index = i;
413  }
414  }
415  }
416 
417  m_reg1to2 = 1;
418  if (maxReg1Index >= 0 && maxReg2Index >= 0) {
419  if (m_cellTime[maxReg2Index] < m_cellTime[maxReg1Index]) m_reg1to2 = 0;
420  }
421 
422  if (!m_beamType.compare("collisions")) {
423  selection = (maxReg1Energy > m_eThreshold);
424  } else {
425  selection = (maxReg2Energy > m_eThreshold && maxReg1Energy > m_eThreshold);
426  }
427 
428  return selection;
429 }
430 
431 // ********************************************************************
432 // ********************************************************************
433 // ********************************************************************
434 
435 int TileMuonFitter::fitTrack() {
436 
437  ATH_MSG_DEBUG( "Fitting with Std method" << m_nCells << " cells." );
438 
439  std::vector<double> X;
440  std::vector<double> Y;
441  std::vector<double> Z;
442  X.resize(m_nCells);
443  Y.resize(m_nCells);
444  Z.resize(m_nCells);
445  for (int i = 0; i < m_nCells; i++) {
446 
447  X[i] = m_cellPosition[i].getX();
448  Y[i] = m_cellPosition[i].getY();
449  Z[i] = m_cellPosition[i].getZ();
450  double rr = sqrt(X[i] * X[i] + Y[i] * Y[i]);
451 
452  ATH_MSG_DEBUG( "Cell " << i
453  << " Energy " << m_cellEnergy[i]
454  << " Weight " << m_cellWeight[i]
455  << " Time " << m_cellTime[i]
456  << " r " << rr
457  << " X " << X[i]
458  << " Y " << Y[i]
459  << " Z " << Z[i] );
460 
461  }
462 
463  // Minuit Method ------------------------------------------------------------
464 
465  m_theTrack = new TileMuonTrackDistance(X, Y, Z, m_cellWeight);
466  m_theTrack->SetWeighted(m_doWeighted);
467  m_theTrack->Means();
468 
469  double aa, bb, cc, dd;
470  double erraa, errbb, errcc, errdd;
471  aa = bb = cc = dd = erraa = errbb = errcc = errdd = 0.0;
472  if (m_cellPosition[m_maxTopIndex].getX() == m_cellPosition[m_maxBottomIndex].getX()) {
473  bb = 0.;
474  dd = 0.;
475  } else {
476  bb = (m_cellPosition[m_maxTopIndex].getY() - m_cellPosition[m_maxBottomIndex].getY())
477  / (m_cellPosition[m_maxTopIndex].getX() - m_cellPosition[m_maxBottomIndex].getX());
478  dd = (m_cellPosition[m_maxTopIndex].getZ() - m_cellPosition[m_maxBottomIndex].getZ())
479  / (m_cellPosition[m_maxTopIndex].getX() - m_cellPosition[m_maxBottomIndex].getX());
480  }
481 
482  MnUserParameters upar;
483  upar.Add("bb", bb, 0.1);
484  upar.Add("dd", dd, 0.1);
485 
486  MnMigrad migrad(*m_theTrack, upar);
487  FunctionMinimum min = migrad();
488 
489  bool myIsValid = min.IsValid();
490 
491  m_fitMinimum.push_back(0);
492 
493  if (myIsValid) {
494  bb = min.UserState().Value("bb");
495  dd = min.UserState().Value("dd");
496  aa = m_theTrack->GetMeanY() - bb * m_theTrack->GetMeanX();
497  cc = m_theTrack->GetMeanZ() - dd * m_theTrack->GetMeanX();
498  errbb = min.UserState().Error("bb");
499  errdd = min.UserState().Error("dd");
500  erraa = fabs(errbb * m_theTrack->GetMeanX());
501  errcc = fabs(errdd * m_theTrack->GetMeanX());
502 
503  addTrack(aa, bb, cc, dd);
504 
505  m_fitMinimum[m_fitMinimum.size() - 1] = min.UserState().Fval();
506  }
507 
508  if (myIsValid) {
509  if (msgLvl(MSG::DEBUG)) {
510  msg(MSG::DEBUG) << "Minimum Value aa: " << aa
511  << " bb: " << bb
512  << " cc: " << cc
513  << " dd " << dd << endmsg;
514  msg(MSG::DEBUG) << "Minimum Error aa: " << erraa
515  << " bb: " << errbb
516  << " cc: " << errcc
517  << " dd " << errdd << endmsg;
518  }
519  return 1;
520  } else {
521  ATH_MSG_DEBUG( "Could not minimize" );
522  return 0;
523  }
524 
525 }
526 // ********************************************************************
527 // ********************************************************************
528 // ********************************************************************
529 
530 void TileMuonFitter::calculateTime() {
531 
532  if (!m_beamType.compare("singlebeam") || !m_beamType.compare("collisions"))
533  calculateTimeAtZequal0();
534  else
535  calculateTimeAtYequal0();
536 
537 }
538 // ********************************************************************
539 // ********************************************************************
540 // ********************************************************************
541 
542 void TileMuonFitter::calculateTimeAtYequal0() {
543 
544  // Uses the track parameters to calculate:
545  // time, position at horizontal plane
546  // Since y = a + b*x and z = c + d*x
547  // y == 0 => x = -a/b && z = c - d*a/b
548  // Time is average from each corrected with tof along track
549 
550 
551  for (unsigned int n = 0; n < m_linePar.size(); n++) {
552  double p0 = m_linePar[n][0];
553  double p1 = m_linePar[n][1];
554  double p2 = m_linePar[n][2];
555  double p3 = m_linePar[n][3];
556 
557  ATH_MSG_DEBUG( "Starting CalculateTime at Y=0 m_linePar: "
558  << p0 << " " << p1 << " " << p2 << " " << p3 );
559 
560  double correction = 0, weight = 1.0, weightsum = 0, time = 0;
561 
562  if (p1 == 0) {
563  m_zeroCrossingTime.push_back(0.0);
564  } else {
565  Hep3Vector theZeroCrossingPosition(-1.0 * p0 / p1, 0.0, p2 - 1.0 * p3 * p0 / p1);
566 
567  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Using cells to calculate time at Y=0: ";
568 
569  for (int i = 0; i < m_nCells; i++) {
570  if (m_doHoughTransform && m_cellInfo[i].track_index != (int) n) continue; // use only cells that belong to the current track
571 
572  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << i << "/ ";
573  Hep3Vector dataP(m_cellPosition[i].getX(), m_cellPosition[i].getY(),
574  m_cellPosition[i].getZ());
575  Hep3Vector lineP(m_theTrack->ClosestPoint(&dataP, m_linePar[n]));
576  correction = (theZeroCrossingPosition - lineP).mag() / c_light;
577 
578  if (m_doWeighted) weight = m_cellEnergy[i];
579  if ((m_reg1to2 && (m_cellPosition[i].getY() > 0))
580  || (!m_reg1to2 && (m_cellPosition[i].getY() < 0))) {
581  correction = 1.0 * correction;
582  } else {
583  correction = -1.0 * correction;
584  }
585 
586  time += weight * (m_cellTime[i] + correction);
587  weightsum += weight;
588  }
589 
590  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << endmsg;
591 
592  if (weightsum > 0)
593  m_zeroCrossingTime.push_back(time / weightsum);
594  else
595  m_zeroCrossingTime.push_back(0.0);
596 
597  ATH_MSG_DEBUG( "Time at Y=0: " << m_zeroCrossingTime[n] );
598  }
599  }
600 }
601 // ********************************************************************
602 // ********************************************************************
603 // ********************************************************************
604 
605 void TileMuonFitter::calculateTimeAtZequal0() {
606 
607  // Uses the track parameters to calculate:
608  // time, position at vertical plane
609  // Since y = a + b*x and z = c + d*x
610  // z == 0 => x = -c/d && y = a - b*c/d
611  // Time is average from each corrected with tof along track
612 
613 
614  for (unsigned int n = 0; n < m_linePar.size(); n++) {
615  double p0 = m_linePar[n][0];
616  double p1 = m_linePar[n][1];
617  double p2 = m_linePar[n][2];
618  double p3 = m_linePar[n][3];
619 
620  ATH_MSG_DEBUG( "Starting CalculateTime at Z=0 m_linePar: "
621  << p0 << " " << p1 << " " << p2 << " " << p3 );
622 
623  double correction = 0, weight = 1.0, weightsum = 0, time = 0;
624 
625  if (p3 == 0) {
626  m_zeroCrossingTime.push_back(0.0);
627  } else {
628  Hep3Vector theZeroCrossingPosition(-1.0 * p2 / p3, p0 - 1.0 * p1 * p2 / p3, 0.0);
629 
630  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Using cells to calculate time at Z=0: ";
631 
632  for (int i = 0; i < m_nCells; i++) {
633  if (m_doHoughTransform && m_cellInfo[i].track_index != (int) n) continue; // use only cells that belong to the current track
634 
635  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << i << "/ ";
636  Hep3Vector dataP(m_cellPosition[i].getX(), m_cellPosition[i].getY(),
637  m_cellPosition[i].getZ());
638  Hep3Vector lineP(m_theTrack->ClosestPoint(&dataP, m_linePar[n]));
639  correction = (theZeroCrossingPosition - lineP).mag() / c_light;
640 
641  if (m_doWeighted) weight = m_cellEnergy[i];
642  if ((m_reg1to2 && (m_cellPosition[i].getZ() > 0))
643  || (!m_reg1to2 && (m_cellPosition[i].getZ() < 0)))
644  correction = 1.0 * correction;
645  else
646  correction = -1.0 * correction;
647 
648  time += weight * (m_cellTime[i] + correction);
649  weightsum += weight;
650  }
651 
652  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << endmsg;
653 
654  if (weightsum > 0)
655  m_zeroCrossingTime.push_back(time / weightsum);
656  else
657  m_zeroCrossingTime.push_back(0.0);
658 
659  ATH_MSG_DEBUG( "Time: " << m_zeroCrossingTime[n] );
660  }
661  }
662 }
663 
664 // ********************************************************************
665 // ********************************************************************
666 // ********************************************************************
667 
668 void TileMuonFitter::buildTileCosmicMuon(int fitok) {
669 
670  if (!m_beamType.compare("singlebeam") || !m_beamType.compare("collisions"))
671  buildTileCosmicMuonAtZequal0(fitok);
672  else
673  buildTileCosmicMuonAtYequal0(fitok);
674 
675 }
676 
677 // ********************************************************************
678 // ********************************************************************
679 // ********************************************************************
680 
681 void TileMuonFitter::buildTileCosmicMuonAtYequal0(int fitok) {
682 
683  // Uses the track parameters to calculate:
684  // position and direction at horizontal at horizontal plane
685  // Since y = a + b*x and z = c + d*x
686  // y == 0 => x = -a/b && z = c - d*a/b
687  // Direction is (1,b,d)/sqrt(1 + b*b + d*d)
688  //
689  // Then put the result in a ComTime object
690 
691  SG::WriteHandle<TileCosmicMuonContainer> cosmicMuonContainer (m_cosmicMuonContainerKey);
692  StatusCode sc = cosmicMuonContainer.record(std::make_unique<TileCosmicMuonContainer>());
693  if (sc.isFailure()) {
694  ATH_MSG_FATAL( "Cannot record TileCosmicMuonContainer" << cosmicMuonContainer.key() );
695  }
696 
697  for (unsigned int n = 0; n < m_linePar.size(); n++) {
698  double p0 = m_linePar[n][0];
699  double p1 = m_linePar[n][1];
700  double p2 = m_linePar[n][2];
701  double p3 = m_linePar[n][3];
702 
703  ATH_MSG_DEBUG( "Starting BuildTileCosmicMuon at y=0 m_linePar: "
704  << p0 << " " << p1 << " " << p2 << " " << p3 );
705 
706  TileCosmicMuon *theTileCosmicMuon = new TileCosmicMuon();
707 
708  if (!(fitok == 1) || p1 == 0.0) {
709  cosmicMuonContainer->push_back(theTileCosmicMuon);
710  } else {
711  Hep3Vector theDirection(1.0, p1, p3);
712  theDirection = theDirection.unit();
713  if ((m_reg1to2 && p1 > 0.0) || (!m_reg1to2 && p1 < 0.0)) theDirection *= -1.0;
714 
715  std::vector<double> ltop, lbot, etop, ebot;
716  std::vector<IdentifierHash> cells;
717  std::vector<double> segPath;
718  std::vector<int> segPartition, segModule, segSampling;
719 
720  if (m_linePar.size() > 0) trackIntersection(ltop, lbot, n);
721  if (m_linePar.size() > 0)
722  trackSegmentIntersection(segPath, segPartition, segModule, segSampling, n);
723  if (m_linePar.size() > 0) energyInTrack(etop, ebot, cells, n);
724 
725  theTileCosmicMuon->SetTime(m_zeroCrossingTime[n]);
726  theTileCosmicMuon->SetPositionX(-1.0 * p0 / p1);
727  theTileCosmicMuon->SetPositionY(0.);
728  theTileCosmicMuon->SetPositionZ(p2 - 1.0 * p3 * p0 / p1);
729  theTileCosmicMuon->SetDirectionPhi(theDirection.phi());
730  theTileCosmicMuon->SetDirectionTheta(theDirection.theta());
731  theTileCosmicMuon->SetFitQuality(m_fitMinimum[n]);
732  theTileCosmicMuon->SetFitNCells(m_nCells);
733  theTileCosmicMuon->SetPathTop(ltop);
734  theTileCosmicMuon->SetPathBottom(lbot);
735  theTileCosmicMuon->SetSegmentPath(segPath);
736  theTileCosmicMuon->SetSegmentPartition(segPartition);
737  theTileCosmicMuon->SetSegmentModule(segModule);
738  theTileCosmicMuon->SetSegmentSampling(segSampling);
739  theTileCosmicMuon->SetEnergyTop(etop);
740  theTileCosmicMuon->SetEnergyBottom(ebot);
741  theTileCosmicMuon->SetTrackCellHash(cells);
742 
743  cosmicMuonContainer->push_back(theTileCosmicMuon);
744  }
745 
746  }
747 
748  if (fitok != 2) delete m_theTrack;
749 }
750 // ********************************************************************
751 // ********************************************************************
752 // ********************************************************************
753 
754 void TileMuonFitter::buildTileCosmicMuonAtZequal0(int fitok) {
755 
756  // Uses the track parameters to calculate:
757  // position and direction at horizontal at horizontal plane
758  // Since y = a + b*x and z = c + d*x
759  // z == 0 => x = -c/d && y = a - b*c/d
760  // Direction is (1,b,d)/sqrt(1 + b*b + d*d)
761  //
762  // Then put the result in a ComTime object
763 
764  SG::WriteHandle<TileCosmicMuonContainer> cosmicMuonContainer (m_cosmicMuonContainerKey);
765  StatusCode sc = cosmicMuonContainer.record(std::make_unique<TileCosmicMuonContainer>());
766  if (sc.isFailure()) {
767  ATH_MSG_FATAL( "Cannot record TileCosmicMuonContainer" << cosmicMuonContainer.key() );
768  }
769 
770  for (unsigned int n = 0; n < m_linePar.size(); n++) {
771  double p0 = m_linePar[n][0];
772  double p1 = m_linePar[n][1];
773  double p2 = m_linePar[n][2];
774  double p3 = m_linePar[n][3];
775 
776  ATH_MSG_DEBUG( "Starting BuildTileCosmicMuon at z=0 m_linePar: "
777  << p0 << " " << p1 << " " << p2 << " " << p3 );
778 
779  TileCosmicMuon *theTileCosmicMuon = new TileCosmicMuon();
780 
781  if (!(fitok == 1) || p3 == 0.0) {
782  cosmicMuonContainer->push_back(theTileCosmicMuon);
783  } else {
784  Hep3Vector theDirection(1.0, p1, p3);
785  theDirection = theDirection.unit();
786  if ((m_reg1to2 && p3 > 0.0) || (!m_reg1to2 && p3 < 0.0)) theDirection *= -1.0;
787 
788  std::vector<double> ltop, lbot, etop, ebot;
789  std::vector<IdentifierHash> cells;
790  std::vector<double> segPath;
791  std::vector<int> segPartition, segModule, segSampling;
792 
793  if (m_linePar.size() > 0) trackIntersection(ltop, lbot, n);
794  if (m_linePar.size() > 0)
795  trackSegmentIntersection(segPath, segPartition, segModule, segSampling, n);
796  if (m_linePar.size() > 0) energyInTrack(etop, ebot, cells, n);
797 
798  theTileCosmicMuon->SetTime(m_zeroCrossingTime[n]);
799  theTileCosmicMuon->SetPositionX(-1.0 * p2 / p3);
800  theTileCosmicMuon->SetPositionY(p0 - 1.0 * p1 * p2 / p3);
801  theTileCosmicMuon->SetPositionZ(0.);
802  theTileCosmicMuon->SetDirectionPhi(theDirection.phi());
803  theTileCosmicMuon->SetDirectionTheta(theDirection.theta());
804  theTileCosmicMuon->SetFitQuality(m_fitMinimum[n]);
805  theTileCosmicMuon->SetFitNCells(m_nCells);
806  theTileCosmicMuon->SetPathTop(ltop);
807  theTileCosmicMuon->SetPathBottom(lbot);
808  theTileCosmicMuon->SetSegmentPath(segPath);
809  theTileCosmicMuon->SetSegmentPartition(segPartition);
810  theTileCosmicMuon->SetSegmentModule(segModule);
811  theTileCosmicMuon->SetSegmentSampling(segSampling);
812  theTileCosmicMuon->SetEnergyTop(etop);
813  theTileCosmicMuon->SetEnergyBottom(ebot);
814  theTileCosmicMuon->SetTrackCellHash(cells);
815 
816  cosmicMuonContainer->push_back(theTileCosmicMuon);
817  }
818 
819  }
820 
821  if (fitok != 2) delete m_theTrack;
822 }
823 
824 // ********************************************************************
825 // ********************************************************************
826 // ********************************************************************
827 void TileMuonFitter::trackIntersection(std::vector<double> & ltop, std::vector<double> & lbot,
828  int index) {
829 
830  const double p0 = m_linePar[index][0];
831  const double p1 = m_linePar[index][1];
832  const double p2 = m_linePar[index][2];
833  const double p3 = m_linePar[index][3];
834 
835  ATH_MSG_DEBUG( "Starting TrackIntersection: " );
836  int is, ip;
837  std::vector<Hep3Vector> intPoint; //intersection points
838  std::vector<Hep3Vector> topIPO; //intersection points on top hemisphere, order according to y
839  std::vector<Hep3Vector> bottomIPO; //intersection points on bottom hemisphere, order according to y
840  std::vector<Hep3Vector> temp3v;
841  std::vector<Hep3Vector> plane;
842  Hep3Vector horizontalPlane;
843 
844  plane.resize(6);
845  lbot.resize(3);
846  ltop.resize(3);
847  temp3v.resize(2);
848  for (is = 0; is < 3; is++) {
849  ltop[is] = 0.0;
850  lbot[is] = 0.0;
851  }
853  for (is = 0; is < 4; is++) {
854 
856  double aux_squared = m_tileDD_radiusLB[is] * m_tileDD_radiusLB[is] * (1 + p1 * p1) - p0 * p0;
857  if (aux_squared > 0) {
858  double x[2] = { (-1.0 * p0 * p1 + sqrt(aux_squared)) / (1 + p1 * p1), (-1.0 * p0 * p1
859  - sqrt(aux_squared)) / (1 + p1 * p1) };
860  for (uint8_t i = 0; i < 2; i++) {
861  temp3v[i].set(x[i], p0 + p1 * x[i], p2 + p3 * x[i]);
862  if (checkLBz(temp3v[i].z())) intPoint.push_back(temp3v[i]);
863  }
864  }
865 
867  aux_squared = m_tileDD_radiusEB[is] * m_tileDD_radiusEB[is] * (1 + p1 * p1) - p0 * p0;
868  if (aux_squared > 0) {
869  double x[2] = { (-1.0 * p0 * p1 + sqrt(aux_squared)) / (1 + p1 * p1), (-1.0 * p0 * p1
870  - sqrt(aux_squared)) / (1 + p1 * p1) };
871  for (uint8_t i = 0; i < 2; i++) {
872  temp3v[i].set(x[i], p0 + p1 * x[i], p2 + p3 * x[i]);
873  if (checkEBz(temp3v[i].z())) intPoint.push_back(temp3v[i]);
874  }
875  }
876 
877  }
878  uint8_t ncylint = intPoint.size();
879 
881  for (ip = 0; ip < 6; ip++)
882  plane[ip].set(0, 0, 0);
883  if (p3 != 0) {
884  const double inv_p3 = 1. / p3;
885  for (ip = 0; ip < 2; ip++) {
886  plane[ip].set((m_tileDD_zLB[ip] - p2) * inv_p3, p0 + p1 * (m_tileDD_zLB[ip] - p2) * inv_p3,
887  m_tileDD_zLB[ip]);
888 
889  plane[ip + 2].set((m_tileDD_zEBA[ip] - p2) * inv_p3, p0 + p1 * (m_tileDD_zEBA[ip] - p2) * inv_p3,
890  m_tileDD_zEBA[ip]);
891 
892  plane[ip + 4].set((m_tileDD_zEBC[ip] - p2) * inv_p3, p0 + p1 * (m_tileDD_zEBC[ip] - p2) * inv_p3,
893  m_tileDD_zEBC[ip]);
894 
895  for (is = 0; is < 3; is++) {
896  if (checkLBr(plane[ip].perp(), is)) intPoint.push_back(plane[ip]);
897  if (checkEBr(plane[ip + 2].perp(), is)) intPoint.push_back(plane[ip + 2]);
898  if (checkEBr(plane[ip + 4].perp(), is)) intPoint.push_back(plane[ip + 4]);
899  }
900  }
901  }
902 
903  ATH_MSG_DEBUG( "Intersection with cylinders: " << ncylint
904  << " with planes: " << intPoint.size() - ncylint );
905 
909  horizontalPlane.set(0, 0, 0);
910  uint16_t idx_hor = 99;
911  uint16_t i, j, k, jmax;
912  if (p1 != 0) {
913  const double inv_p1 = 1. / p1;
915  horizontalPlane.set(-1.0 * p0 * inv_p1, 0, p2 - p3 * p0 * inv_p1);
916  if (checkLBz(horizontalPlane.z()) && checkLBr(horizontalPlane.x())) {
917  for (i = 0; i < 3; i++)
918  if (checkLBr(horizontalPlane.x(), i)) idx_hor = i;
919  }
920  if (checkEBz(horizontalPlane.z()) && checkEBr(horizontalPlane.x())) {
921  for (i = 0; i < 3; i++)
922  if (checkEBr(horizontalPlane.x(), i)) idx_hor = i;
923  }
924  }
925 
926  ATH_MSG_DEBUG( "Horizontal Plane idx: " << idx_hor
927  << " X: " << horizontalPlane.x()
928  << " Z: " << horizontalPlane.z() );
929 
932  uint16_t npoints = intPoint.size();
933  float ymax;
934  bool filled;
935  std::vector<uint8_t> neworder;
936  neworder.resize(0);
937  for (i = 0; i < npoints; i++) {
938  ymax = -99999999.0;
939  jmax = 9999;
940  for (j = 0; j < npoints; j++) {
941  filled = false;
942  for (k = 0; k < neworder.size(); k++) {
943  if (j == neworder[k]) filled = true;
944  }
945  if (filled) continue;
946  if (intPoint[j].y() > ymax) {
947  ymax = intPoint[j].y();
948  jmax = j;
949  }
950  }
951  if (jmax < 9999) neworder.push_back(jmax);
952  }
953 
954  if (neworder.size() != npoints)
955  ATH_MSG_ERROR( " Npoints " << npoints
956  << " New order: " << neworder.size() );
957 
962  if (idx_hor < 99) bottomIPO.push_back(horizontalPlane);
963 
964  for (i = 0; i < npoints; i++) {
965  if (intPoint[neworder[i]].y() > 0)
966  topIPO.push_back(intPoint[neworder[i]]);
967  else
968  bottomIPO.push_back(intPoint[neworder[i]]);
969  }
970  if (idx_hor < 99) topIPO.push_back(horizontalPlane);
971 
972  ATH_MSG_DEBUG( "Number of points on top: " << topIPO.size()
973  << " and on bottom: " << bottomIPO.size() );
974 
975  if (topIPO.size() > 1 && bottomIPO.size()) {
976 
983  Hep3Vector temp_midpoint;
984  int temp_idx;
985  for (i = 0; i < topIPO.size() - 1; i++) {
986  ATH_MSG_DEBUG( "Top points : " << i
987  << " X: " << topIPO[i].x()
988  << " Y: " << topIPO[i].y()
989  << " Z: " << topIPO[i].z() );
990 
991  temp_midpoint = (topIPO[i] + topIPO[i + 1]) / 2.;
992  ATH_MSG_DEBUG( "Temp_midz top: " << temp_midpoint.z() );
993 
994  if (checkLBz(temp_midpoint.z()))
995  temp_idx = whichLBr(temp_midpoint.perp());
996  else if (checkEBz(temp_midpoint.z()))
997  temp_idx = whichEBr(temp_midpoint.perp());
998  else
999  continue;
1000 
1001  ATH_MSG_DEBUG( "Temp_idx : " << temp_idx );
1002 
1003  if (temp_idx < 0 || temp_idx > 2) continue;
1004  ltop[temp_idx] += (topIPO[i] - topIPO[i + 1]).mag();
1005 
1006  ATH_MSG_DEBUG( "ltop : " << ltop[temp_idx] );
1007 
1008  }
1009 
1010  ATH_MSG_DEBUG( "Top points : " << topIPO.size() - 1
1011  << " X: " << topIPO[topIPO.size() - 1].x()
1012  << " Y: " << topIPO[topIPO.size() - 1].y()
1013  << " Z: " << topIPO[topIPO.size() - 1].z() );
1014 
1015 
1016  for (i = 0; i < bottomIPO.size() - 1; i++) {
1017  ATH_MSG_DEBUG( "Bot points : " << i
1018  << " X: " << bottomIPO[i].x()
1019  << " Y: " << bottomIPO[i].y()
1020  << " Z: " << bottomIPO[i].z() );
1021 
1022 
1023  temp_midpoint = (bottomIPO[i] + bottomIPO[i + 1]) / 2.;
1024  ATH_MSG_DEBUG( "Temp_midz bottom: " << temp_midpoint.z() );
1025 
1026  if (checkLBz(temp_midpoint.z()))
1027  temp_idx = whichLBr(temp_midpoint.perp());
1028  else if (checkEBz(temp_midpoint.z()))
1029  temp_idx = whichEBr(temp_midpoint.perp());
1030  else
1031  continue;
1032 
1033  ATH_MSG_DEBUG( "Temp_idx : " << temp_idx );
1034 
1035  if (temp_idx < 0 || temp_idx > 2) continue;
1036  lbot[temp_idx] += (bottomIPO[i] - bottomIPO[i + 1]).mag();
1037 
1038  ATH_MSG_DEBUG( "lbot : " << lbot[temp_idx] );
1039  }
1040 
1041  ATH_MSG_DEBUG( "Bot points : " << bottomIPO.size() - 1
1042  << " X: " << bottomIPO[bottomIPO.size() - 1].x()
1043  << " Y: " << bottomIPO[bottomIPO.size() - 1].y()
1044  << " Z: " << bottomIPO[bottomIPO.size() - 1].z() );
1045 
1046  }
1047  ATH_MSG_DEBUG( "Leaving TrackIntersection" );
1048 }
1049 
1050 // ********************************************************************
1051 // ********************************************************************
1052 // ********************************************************************
1053 void TileMuonFitter::trackSegmentIntersection(std::vector<double> & segPath,
1054  std::vector<int> & segPartition, std::vector<int> & segModule, std::vector<int> & segSampling,
1055  int index) {
1056 
1057  double p0 = m_linePar[index][0];
1058  double p1 = m_linePar[index][1];
1059  double p2 = m_linePar[index][2];
1060  double p3 = m_linePar[index][3];
1061 
1062  ATH_MSG_DEBUG( "Starting TrackIntersection: " );
1063 
1064  std::vector<Hep3Vector> intPoint; //intersection points
1065  std::vector<Hep3Vector> ordPoint; //ordered intersection points (according to y)
1066 
1067  std::vector<Hep3Vector> temp3v;
1068  std::vector<Hep3Vector> plane;
1069  std::vector<Hep3Vector> modPlane;
1070 
1071  plane.resize(6);
1072  modPlane.clear();
1073  modPlane.resize(32);
1074  temp3v.resize(2);
1075  segPath.clear();
1076  segPartition.clear();
1077  segModule.clear();
1078  segSampling.clear();
1079 
1081  int is, ip;
1082  uint16_t i, j, k, jmax;
1083  for (is = 0; is < 4; is++) {
1084 
1086  double aux_squared = m_tileDD_radiusLB[is] * m_tileDD_radiusLB[is] * (1 + p1 * p1) - p0 * p0;
1087  if (aux_squared > 0) {
1088  double x[2] = { (-1.0 * p0 * p1 + sqrt(aux_squared)) / (1 + p1 * p1), (-1.0 * p0 * p1
1089  - sqrt(aux_squared)) / (1 + p1 * p1) };
1090  for (i = 0; i < 2; i++) {
1091  temp3v[i].set(x[i], p0 + p1 * x[i], p2 + p3 * x[i]);
1092  if (checkLBz(temp3v[i].z())) intPoint.push_back(temp3v[i]);
1093  }
1094  }
1095 
1097  aux_squared = m_tileDD_radiusEB[is] * m_tileDD_radiusEB[is] * (1 + p1 * p1) - p0 * p0;
1098  if (aux_squared > 0) {
1099  double x[2] = { (-1.0 * p0 * p1 + sqrt(aux_squared)) / (1 + p1 * p1), (-1.0 * p0 * p1
1100  - sqrt(aux_squared)) / (1 + p1 * p1) };
1101  for (i = 0; i < 2; i++) {
1102  temp3v[i].set(x[i], p0 + p1 * x[i], p2 + p3 * x[i]);
1103  if (checkEBz(temp3v[i].z())) intPoint.push_back(temp3v[i]);
1104  }
1105  }
1106 
1107  }
1108  uint16_t ncylint = intPoint.size();
1109 
1111  for (ip = 0; ip < 6; ip++)
1112  plane[ip].set(0, 0, 0);
1113  if (p3 != 0) {
1114  for (ip = 0; ip < 2; ip++) {
1115  plane[ip].set((m_tileDD_zLB[ip] - p2) / p3, p0 + p1 * (m_tileDD_zLB[ip] - p2) / p3,
1116  m_tileDD_zLB[ip]);
1117 
1118  plane[ip + 2].set((m_tileDD_zEBA[ip] - p2) / p3, p0 + p1 * (m_tileDD_zEBA[ip] - p2) / p3,
1119  m_tileDD_zEBA[ip]);
1120 
1121  plane[ip + 4].set((m_tileDD_zEBC[ip] - p2) / p3, p0 + p1 * (m_tileDD_zEBC[ip] - p2) / p3,
1122  m_tileDD_zEBC[ip]);
1123 
1124  for (is = 0; is < 3; is++) {
1125  if (checkLBr(plane[ip].perp(), is)) intPoint.push_back(plane[ip]);
1126  if (checkEBr(plane[ip + 2].perp(), is)) intPoint.push_back(plane[ip + 2]);
1127  if (checkEBr(plane[ip + 4].perp(), is)) intPoint.push_back(plane[ip + 4]);
1128  }
1129  }
1130  }
1131  uint16_t nvpint = intPoint.size() - ncylint;
1132 
1134  double tanTheta, tempX;
1135  for (ip = 0; ip < 32; ip++) {
1136  tanTheta = tan(pi / 32. * double(ip));
1137  if ((tanTheta - p1) != 0) {
1138  tempX = p0 / (tanTheta - p1);
1139  modPlane[ip].set(tempX, p0 + p1 * tempX, p2 + p3 * tempX);
1140 
1141  if (checkLBz(modPlane[ip].z()) && checkLBr(modPlane[ip].perp())) {
1142  intPoint.push_back(modPlane[ip]);
1143  }
1144 
1145  if (checkEBz(modPlane[ip].z()) && checkEBr(modPlane[ip].perp())) {
1146  intPoint.push_back(modPlane[ip]);
1147  }
1148  }
1149  }
1150  uint16_t nmpint = intPoint.size() - ncylint - nvpint;
1151  uint16_t npoints = intPoint.size();
1152 
1153  ATH_MSG_DEBUG( "Intersections: " << npoints
1154  << " cylinders: " << ncylint
1155  << " vertical planes: " << nvpint
1156  << " module planes: " << nmpint );
1157 
1159  double ymax;
1160  bool filled;
1161  std::vector<int> neworder;
1162  neworder.resize(0);
1163  for (i = 0; i < npoints; i++) {
1164  ymax = -999999.0;
1165  jmax = 9999;
1166  for (j = 0; j < npoints; j++) {
1167  filled = false;
1168  for (k = 0; k < neworder.size(); k++) {
1169  if (j == neworder[k]) filled = true;
1170  }
1171  if (filled) continue;
1172  if (intPoint[j].y() > ymax) {
1173  ymax = intPoint[j].y();
1174  jmax = j;
1175  }
1176  }
1177  if (jmax < 9999) neworder.push_back(jmax);
1178  }
1179 
1180  if (neworder.size() != npoints)
1181  ATH_MSG_DEBUG( " Npoints " << npoints << " New order: " << neworder.size() );
1182 
1183  for (i = 0; i < npoints; i++)
1184  ordPoint.push_back(intPoint[neworder[i]]);
1185 
1186  ATH_MSG_DEBUG( "Number of points: " << ordPoint.size() );
1187 
1194  if (ordPoint.size() > 1) {
1195  Hep3Vector temp_midpoint;
1196  int temp_idr, temp_idp, temp_idm;
1197  bool extra_debug = false;
1198  for (i = 0; i < ordPoint.size(); i++) {
1199  ATH_MSG_DEBUG( "Ordered points : " << i
1200  << " X: " << ordPoint[i].x()
1201  << " Y: " << ordPoint[i].y()
1202  << " Z: " << ordPoint[i].z() );
1203 
1204  if (i == ordPoint.size() - 1) break;
1205 
1206  temp_midpoint = (ordPoint[i] + ordPoint[i + 1]) / 2.;
1207 
1208  if (checkLBz(temp_midpoint.z())) {
1209  temp_idr = whichLBr(temp_midpoint.perp());
1210  temp_idp = (temp_midpoint.z() > 0) ? 1 : 2;
1211  } else if (checkEBz(temp_midpoint.z())) {
1212  temp_idr = whichEBr(temp_midpoint.perp());
1213  temp_idp = (temp_midpoint.z() > 0) ? 3 : 4;
1214  } else
1215  continue;
1216 
1217  temp_idm = whichModule(temp_midpoint);
1218  if (temp_idr < 0 || temp_idr > 2) continue;
1219  if (temp_idm < 0 || temp_idm > 63) continue;
1220 
1221  segPath.push_back((ordPoint[i] - ordPoint[i + 1]).mag());
1222  segPartition.push_back(temp_idp);
1223  segModule.push_back(temp_idm);
1224  segSampling.push_back(temp_idr);
1225 
1226  }
1227 
1228  ATH_MSG_DEBUG( "Number of segments : " << segPath.size() );
1229 
1230  for (i = 0; i < segPath.size(); i++) {
1231  if (segPath[i] > 1100) extra_debug = true;
1232  ATH_MSG_DEBUG( "Segment " << i
1233  << " Path : " << segPath[i]
1234  << " Partition : " << segPartition[i]
1235  << " Module : " << segModule[i]
1236  << " Sampling : " << segSampling[i] );
1237 
1238  }
1239 
1240  if (msgLvl(MSG::DEBUG) && extra_debug) {
1241  msg(MSG::DEBUG) << "Large segment p0: " << p0
1242  << " p1 " << p1
1243  << " p2 " << p2
1244  << " p3 " << p3 << endmsg;
1245 
1246  msg(MSG::DEBUG) << "Intersections: " << npoints
1247  << " cylinders: " << ncylint
1248  << " vertical planes: " << nvpint
1249  << " module planes: " << nmpint << endmsg;
1250 
1251  for (i = 0; i < ordPoint.size(); i++) {
1252  msg(MSG::DEBUG) << "Ordered points : " << i
1253  << " X: " << ordPoint[i].x()
1254  << " Y: " << ordPoint[i].y()
1255  << " Z: " << ordPoint[i].z() << endmsg;
1256  }
1257 
1258  for (i = 0; i < segPath.size(); i++) {
1259  msg(MSG::DEBUG) << "Segment " << i
1260  << " Path : " << segPath[i]
1261  << " Partition : " << segPartition[i]
1262  << " Module : " << segModule[i]
1263  << " Sampling : " << segSampling[i] << endmsg;
1264  }
1265  }
1266 
1267  }
1268 
1269  ATH_MSG_DEBUG( "Leaving TrackSegmentIntersection" );
1270 }
1271 // ********************************************************************
1272 // ********************************************************************
1273 // ********************************************************************
1274 bool TileMuonFitter::checkLBz(double x1) {
1275  if (x1 > m_tileDD_zLB[0] && x1 < m_tileDD_zLB[1])
1276  return true;
1277  else
1278  return false;
1279 }
1280 // ********************************************************************
1281 // ********************************************************************
1282 // ********************************************************************
1283 bool TileMuonFitter::checkEBz(double x1) {
1284  return (checkEBAz(x1) || checkEBCz(x1));
1285 }
1286 // ********************************************************************
1287 // ********************************************************************
1288 // ********************************************************************
1289 bool TileMuonFitter::checkEBAz(double x1) {
1290  if (x1 > m_tileDD_zEBA[0] && x1 < m_tileDD_zEBA[1])
1291  return true;
1292  else
1293  return false;
1294 }
1295 // ********************************************************************
1296 // ********************************************************************
1297 // ********************************************************************
1298 bool TileMuonFitter::checkEBCz(double x1) {
1299  if (x1 > m_tileDD_zEBC[0] && x1 < m_tileDD_zEBC[1])
1300  return true;
1301  else
1302  return false;
1303 }
1304 // ********************************************************************
1305 // ********************************************************************
1306 // ********************************************************************
1307 bool TileMuonFitter::checkLBr(double x1) {
1308  if (x1 < 0) x1 = -1.0 * x1;
1309  if (x1 > m_tileDD_radiusLB[0] && x1 < m_tileDD_radiusLB[3])
1310  return true;
1311  else
1312  return false;
1313 }
1314 // ********************************************************************
1315 // ********************************************************************
1316 // ********************************************************************
1317 bool TileMuonFitter::checkEBr(double x1) {
1318  if (x1 < 0.0) x1 = -1.0 * x1;
1319  if (x1 > m_tileDD_radiusEB[0] && x1 < m_tileDD_radiusEB[3])
1320  return true;
1321  else
1322  return false;
1323 }
1324 // ********************************************************************
1325 // ********************************************************************
1326 // ********************************************************************
1327 bool TileMuonFitter::checkLBr(double x1, uint8_t s1) {
1328  if (s1 > 3U) return false;
1329  if (x1 < 0.0) x1 = -1.0 * x1;
1330  if (x1 > m_tileDD_radiusLB[s1] && x1 < m_tileDD_radiusLB[s1 + 1])
1331  return true;
1332  else
1333  return false;
1334 }
1335 // ********************************************************************
1336 // ********************************************************************
1337 // ********************************************************************
1338 bool TileMuonFitter::checkEBr(double x1, uint8_t s1) {
1339  if (s1 > 3U) return false;
1340  if (x1 < 0.0) x1 = -1.0 * x1;
1341  if (x1 > m_tileDD_radiusEB[s1] && x1 < m_tileDD_radiusEB[s1 + 1])
1342  return true;
1343  else
1344  return false;
1345 }
1346 // ********************************************************************
1347 // ********************************************************************
1348 // ********************************************************************
1349 int TileMuonFitter::whichEBr(double x1) {
1350  int radidx = -1;
1351  for (int i = 0; i < 3; i++) {
1352  if (checkEBr(x1, i)) radidx = i;
1353  }
1354  return radidx;
1355 }
1356 // ********************************************************************
1357 // ********************************************************************
1358 // ********************************************************************
1359 int TileMuonFitter::whichLBr(double x1) {
1360  int radidx = -1;
1361  for (int i = 0; i < 3; i++) {
1362  if (checkLBr(x1, i)) radidx = i;
1363  }
1364  return radidx;
1365 }
1366 // ********************************************************************
1367 // ********************************************************************
1368 // ********************************************************************
1369 int TileMuonFitter::whichModule(Hep3Vector tempvec) {
1370  double phi = tempvec.phi();
1371  int mod;
1372  if (phi < 0.0) phi += 2.0 * pi;
1373  if (phi > 2.0 * pi)
1374  mod = -1;
1375  else
1376  mod = int(phi * (64.0 / (2.0 * pi)));
1377  return mod;
1378 }
1379 // ********************************************************************
1380 // ********************************************************************
1381 // ********************************************************************
1382 void TileMuonFitter::energyInTrack(std::vector<double> & etop, std::vector<double> & ebot
1383  , std::vector<IdentifierHash> & cells, int index) {
1384 
1385  int is;
1386  double distance_cut[2][3] = { { 300., 375., 860. }, { 750., 750., 1700. } };
1387  etop.resize(3);
1388  ebot.resize(3);
1389  for (is = 0; is < 3; is++) {
1390  etop[is] = 0.;
1391  ebot[is] = 0.;
1392  }
1393  Hep3Vector dataP;
1394  Hep3Vector lineP;
1395 
1396  size_t collItr = m_caloCells->indexFirstCellCalo(m_caloIndex);
1397  size_t lastColl = m_caloCells->indexLastCellCalo(m_caloIndex);
1398 
1399  for (; collItr != lastColl; ++collItr) {
1400 
1401  const CaloCell* cell = (*m_caloCells)[collItr];
1402  const TileCell* tilecell = dynamic_cast<const TileCell*>(cell);
1403  if (tilecell == 0) continue;
1404 
1405  Identifier cell_id = cell->ID();
1406  IdentifierHash cell_idhash = m_tileID->cell_hash(cell_id);
1407  CaloDetDescrElement *caloDDE = m_tileMgr->get_cell_element(cell_id);
1408 
1409  double volume = caloDDE->volume();
1410  if (volume == 0.0) {
1411  ATH_MSG_DEBUG( "Warning: Skipping cell with zero volume!" );
1412  continue;
1413  }
1414 
1415  int sample = m_tileID->sample(cell_id);
1416  int section = m_tileID->section(cell_id);
1417  if (section < 1 || section > 2 || sample < 0 || sample >= 3) break;
1418  dataP.set(caloDDE->x(), caloDDE->y(), caloDDE->z());
1419  lineP = m_theTrack->ClosestPoint(&dataP, m_linePar[index]);
1420  double distance = (dataP - lineP).mag();
1421  if (distance < distance_cut[section - 1][sample]) {
1422  cells.push_back(cell_idhash);
1423  if (caloDDE->y() > 0.0) etop[sample] += cell->energy();
1424  else ebot[sample] += cell->energy();
1425  }
1426 
1427  } //end of CaloCellContainer cycle
1428 
1429 }
1430 
1431 // ********************************************************************
1432 // ********************************************************************
1433 // ********************************************************************
1434 
1435 void TileMuonFitter::buildComTime(int fitok) {
1436 
1437  if (!m_beamType.compare("singlebeam") || !m_beamType.compare("collisions"))
1438  buildComTimeAtZequal0(fitok);
1439  else
1440  buildComTimeAtYequal0(fitok);
1441 
1442 }
1443 // ********************************************************************
1444 // ********************************************************************
1445 // ********************************************************************
1446 
1447 void TileMuonFitter::buildComTimeAtYequal0(int fitok) {
1448 
1449  // Uses the track parameters to calculate:
1450  // position and direction at horizontal at horizontal plane
1451  // Since y = a + b*x and z = c + d*x
1452  // y == 0 => x = -a/b && z = c - d*a/b
1453  // Direction is (1,b,d)/sqrt(1 + b*b + d*d)
1454  //
1455  // Then put the result in a ComTime object
1456 
1457 
1458  double p0, p1, p2, p3;
1459  double ztime;
1460 
1461  if (m_linePar.size() == 0) {
1462  p0 = 0.0;
1463  p1 = 0.0;
1464  p2 = 0.0;
1465  p3 = 0.0;
1466  ztime = 0.0;
1467  } else {
1468  p0 = m_linePar[0][0];
1469  p1 = m_linePar[0][1];
1470  p2 = m_linePar[0][2];
1471  p3 = m_linePar[0][3];
1472  ztime = m_zeroCrossingTime[0];
1473  }
1474 
1475  ATH_MSG_DEBUG( "Starting BuildComTime at y=0 m_linePar: "
1476  << p0 << " " << p1 << " " << p2 << " " << p3 );
1477 
1478  SG::WriteHandle<ComTime> theComTime (m_comTimeKey);
1479  StatusCode sc;
1480  if (!(fitok == 1) || p1 == 0.0) {
1481  sc = theComTime.record(std::make_unique<ComTime>());
1482  } else {
1483  Hep3Vector theZeroCrossingPosition(-1.0 * p0 / p1, 0.0, p2 - 1.0 * p3 * p0 / p1);
1484  Hep3Vector theDirection(1.0, p1, p3);
1485  theDirection = theDirection.unit();
1486  if ((m_reg1to2 && p1 > 0.0) || (!m_reg1to2 && p1 < 0.0)) theDirection *= -1.0;
1487 
1488  sc = theComTime.record(std::make_unique<ComTime>(0.0, ztime, theZeroCrossingPosition, theDirection));
1489 
1490  if (msgLvl(MSG::DEBUG)) {
1491 
1492  msg(MSG::DEBUG) << "Time: " << ztime
1493  << " Position X: " << theZeroCrossingPosition.getX()
1494  << " Position Y: " << theZeroCrossingPosition.getY()
1495  << " Position Z: " << theZeroCrossingPosition.getZ() << endmsg;
1496 
1497  msg(MSG::DEBUG) << "Direction u: " << theDirection.getX()
1498  << " Direction v: " << theDirection.getY()
1499  << " Direction w: " << theDirection.getZ() << endmsg;
1500  }
1501  }
1502 
1503  if (sc.isFailure()) {
1504  ATH_MSG_FATAL( "Cannot record ComTime" );
1505  }
1506 }
1507 
1508 // ********************************************************************
1509 // ********************************************************************
1510 // ********************************************************************
1511 
1512 void TileMuonFitter::buildComTimeAtZequal0(int fitok) {
1513 
1514  // Uses the track parameters to calculate:
1515  // position and direction at vertical plane
1516  // Since y = a + b*x and z = c + d*x
1517  // z == 0 => x = -c/d && y = a - b*c/d
1518  // Direction is (1,b,d)/sqrt(1 + b*b + d*d)
1519  //
1520  // Then put the result in a ComTime object
1521 
1522  double p0, p1, p2, p3;
1523  double ztime;
1524 
1525  if (m_linePar.size() == 0) {
1526  p0 = 0.0;
1527  p1 = 0.0;
1528  p2 = 0.0;
1529  p3 = 0.0;
1530  ztime = 0.0;
1531  } else {
1532  p0 = m_linePar[0][0];
1533  p1 = m_linePar[0][1];
1534  p2 = m_linePar[0][2];
1535  p3 = m_linePar[0][3];
1536  ztime = m_zeroCrossingTime[0];
1537  }
1538 
1539  ATH_MSG_DEBUG( "Starting BuildComTime at z=0 m_linePar: "
1540  << p0 << " " << p1 << " " << p2 << " " << p3 );
1541 
1542  SG::WriteHandle<ComTime> theComTime (m_comTimeKey);
1543  StatusCode sc;
1544  if (!(fitok == 1) || p3 == 0.0) {
1545  sc = theComTime.record(std::make_unique<ComTime>());
1546  } else {
1547  Hep3Vector theZeroCrossingPosition(-1.0 * p2 / p3, p0 - 1.0 * p1 * p2 / p3, 0.0);
1548  Hep3Vector theDirection(1.0, p1, p3);
1549  theDirection = theDirection.unit();
1550  if ((m_reg1to2 && p3 > 0.0) || (!m_reg1to2 && p3 < 0.0)) theDirection *= -1.0;
1551 
1552  sc = theComTime.record(std::make_unique<ComTime>(0.0, ztime, theZeroCrossingPosition, theDirection));
1553 
1554  if (msgLvl(MSG::DEBUG)) {
1555 
1556  msg(MSG::DEBUG) << "Time: " << ztime
1557  << " Position X: " << theZeroCrossingPosition.getX()
1558  << " Position Y: " << theZeroCrossingPosition.getY()
1559  << " Position Z: " << theZeroCrossingPosition.getZ() << endmsg;
1560 
1561  msg(MSG::DEBUG) << "Direction u: " << theDirection.getX()
1562  << " Direction v: " << theDirection.getY()
1563  << " Direction w: " << theDirection.getZ() << endmsg;
1564  }
1565  }
1566 
1567  if (sc.isFailure()) {
1568  ATH_MSG_FATAL( "Cannot record ComTime" );
1569  }
1570 }
1571 
1572 // *****************************************************************************
1573 // ************************ Hough Transform Routines ***************************
1574 // *****************************************************************************
1575 
1576 // Cartesian to Hough space
1577 void TileMuonFitter::cart2hough(float x1, float y1, float x2, float y2, double &raio,
1578  double &angu) {
1579 
1580  // track should have at least a very small angle
1581  if (x1 == x2) x1 += 0.001F;
1582  if (y1 == y2) y1 += 0.001F;
1583 
1584  double a = (y1 - y2) / (x1 - x2);
1585  double b = y1 - a * x1;
1586 
1587  raio = fabs(b) / sqrt(a * a + 1.0);
1588 
1589  if (a < 0.0 && b > 0.0)
1590  angu = atan(a) + pi / 2.0;
1591  else if (a > 0.0 && b < 0.0)
1592  angu = atan(a) - pi / 2.0;
1593  else if (a < 0.0 && b < 0.0)
1594  angu = atan(a) - pi / 2.0;
1595  else
1596  angu = atan(a) + pi / 2.0;
1597 
1598  angu = angu - pi;
1599  if (angu < 0.0) angu = 2.0 * pi + angu;
1600 }
1601 
1602 // Hough to Cartesian space
1603 void TileMuonFitter::hough2cart(double r, double a, double offset, double &aa, double &bb) {
1604  double x1, y1;
1605  double x2, y2;
1606 
1607  x1 = (r - sin(a + pi)) / cos(a + pi) - offset; // to translate back
1608  y1 = 1.0;
1609 
1610  x2 = (r + sin(a + pi)) / cos(a + pi) - offset;
1611  y2 = -1.0;
1612 
1613  aa = (y2 - y1) / (x2 - x1);
1614  bb = y1 - aa * x1;
1615 }
1616 
1617 float TileMuonFitter::dist2line(CellInfo &ci, float *pos, float *w) {
1618  float v0, v1, v2, d0, d1, d2;
1619 
1620  v0 = ci.x - pos[0];
1621  v1 = ci.y - pos[1];
1622  v2 = ci.z - pos[2];
1623 
1624  d0 = v1 * w[2] - w[1] * v2;
1625  d1 = w[0] * v2 - v0 * w[2];
1626  d2 = v0 * w[1] - w[0] * v1;
1627 
1628  return sqrt(d0 * d0 + d1 * d1 + d2 * d2);
1629 }
1630 
1631 void TileMuonFitter::points2dir(CellInfo &ci1, CellInfo &ci2, float *w) {
1632  w[0] = ci1.x - ci2.x;
1633  w[1] = ci1.y - ci2.y;
1634  w[2] = ci1.z - ci2.z;
1635 
1636  float invmod = 1.0F / sqrt(w[0] * w[0] + w[1] * w[1] + w[2] * w[2]);
1637  w[0] *= invmod;
1638  w[1] *= invmod;
1639  w[2] *= invmod;
1640 }
1641 
1642 // count number of cells inside RoI (line between two cells)
1643 unsigned int TileMuonFitter::CntCells(unsigned int index1, unsigned int index2, double &skew) {
1644 
1645  unsigned int size = m_cellInfo.size();
1646  double dist;
1647 
1648  // compute direction -------------------------------------------------------
1649 
1650  float p[3] = { m_cellInfo[index1].x, m_cellInfo[index1].y, m_cellInfo[index1].z };
1651  float w[3];
1652  points2dir(m_cellInfo[index1], m_cellInfo[index2], w);
1653 
1654  // find number of cells inside RoI -----------------------------------------
1655 
1656  unsigned int cnt = 0U;
1657  skew = 0.0;
1658  for (unsigned int i = 0U; i < size; i++) {
1659  if (!m_cellInfo[i].use || i == index1 || i == index2) continue;
1660 
1661  dist = dist2line(m_cellInfo[i], p, w);
1662  if (dist < m_cellInfo[i].dist * 0.5F) {
1663  cnt++;
1664  skew += dist;
1665  }
1666  }
1667  if (cnt > 0U) skew /= (int) cnt;
1668 
1669  return cnt;
1670 }
1671 
1672 // get the line between two cells with the highest number of cells inside RoI
1673 bool TileMuonFitter::guessTrack(unsigned int &index1, unsigned int &index2) {
1674  unsigned int i, j;
1675  unsigned int NPoints = m_cellInfo.size();
1676 
1677  double skew, skew_min = 999999999.9;
1678  unsigned int cnt, cnt_max = 0U;
1679 
1680  for (i = 1; i < NPoints; i++) {
1681  for (j = 0; j < i; j++) {
1682  if (m_cellInfo[i].use && m_cellInfo[j].use) {
1683  cnt = CntCells(i, j, skew);
1684  if (cnt > cnt_max || (cnt == cnt_max && skew < skew_min)) {
1685  cnt_max = cnt;
1686  skew_min = skew;
1687  index1 = i;
1688  index2 = j;
1689  }
1690  }
1691  }
1692  }
1693 
1694  return (cnt_max >= MIN_NCELLS);
1695 }
1696 
1697 // build vector of CellInfo with all activated cells
1698 unsigned int TileMuonFitter::buildCellInfoVector() {
1699  const float distance_cut[2][3] = { { 300.0, 820.0, 470.0 }, { 350.0, 540.0, 740.0 } };
1700 
1701  m_cellInfo.clear();
1702  m_linePar.clear();
1703 
1704  for (int i = 0; i < m_nCells; i++) {
1705  CellInfo cell;
1706 
1707  cell.x = m_cellPosition[i].getX() + SHIFT_X;
1708  cell.y = m_cellPosition[i].getY();
1709  cell.z = m_cellPosition[i].getZ() + SHIFT_Z;
1710 
1711  cell.e = m_cellEnergy[i];
1712  cell.ev = m_cellWeight[i];
1713 
1714  cell.use = true;
1715  cell.is_out = true;
1716  cell.track_index = -1;
1717 
1718  int sample = m_tileID->sample(m_tileID->cell_id(m_cellHash[i]));
1719  int section = m_tileID->section(m_tileID->cell_id(m_cellHash[i]));
1720  if (section < 1 || section > 2 || sample < 0 || sample >= 3)
1721  cell.dist = 500.0;
1722  else
1723  cell.dist = distance_cut[section - 1][sample];
1724 
1725  m_cellInfo.push_back(cell);
1726  }
1727 
1728  return m_cellInfo.size();
1729 }
1730 
1731 // select cells that are inside RoI
1732 float TileMuonFitter::selectCells(float *p, float *w) {
1733  unsigned int NPoints = m_cellInfo.size();
1734  float toten = 0.0;
1735  unsigned int ntrack = m_linePar.size();
1736 
1737  for (unsigned int i = 0; i < NPoints; i++) {
1738  if (m_cellInfo[i].use == false || dist2line(m_cellInfo[i], p, w) > m_cellInfo[i].dist) {
1739  m_cellInfo[i].is_out = true;
1740  } else {
1741  m_cellInfo[i].is_out = false;
1742  m_cellInfo[i].use = false;
1743  m_cellInfo[i].track_index = ntrack;
1744 
1745  toten += m_cellInfo[i].e;
1746  }
1747  }
1748 
1749  return toten;
1750 }
1751 
1752 // if line is close to horizontal in plane zy, it should be halo event
1753 bool TileMuonFitter::isHaloMuon(double azy) {
1754  if ((azy > 1.0 * pi / 2.0 - 6.0 * BIN_RES_AZY && azy < 1.0 * pi / 2.0 + 6.0 * BIN_RES_AZY)||
1755  (azy > 3.0*pi/2.0 - 6.0*BIN_RES_AZY && azy < 3.0*pi/2.0 + 6.0*BIN_RES_AZY))return true;
1756  else return false;
1757 }
1758 
1759 void TileMuonFitter::doHough(double &rxy, double &axy, double &rzy, double &azy) {
1760  unsigned int i, j;
1761 
1762  // dont need to compute other iteration for Halo events
1763  if (isHaloMuon(azy)) return;
1764 
1765  float nminrxy = rxy - 4.0 * BIN_RES_RXY;
1766  float nmaxrxy = rxy + 4.0 * BIN_RES_RXY;
1767  float nminaxy = axy - 4.0 * BIN_RES_AXY;
1768  float nmaxaxy = axy + 4.0 * BIN_RES_AXY;
1769  float nminrzy = rzy - 5.0 * BIN_RES_RZY;
1770  float nmaxrzy = rzy + 5.0 * BIN_RES_RZY;
1771  float nminazy = azy - 5.0 * BIN_RES_AZY;
1772  float nmaxazy = azy + 5.0 * BIN_RES_AZY;
1773 
1774  float weight;
1775  float aw = 0.0F;
1776  float arxy = 0.0F;
1777  float aaxy = 0.0F;
1778  float arzy = 0.0F;
1779  float aazy = 0.0F;
1780 
1781  unsigned int NPoints = m_cellInfo.size();
1782  for (i = 1; i < NPoints; i++) {
1783  for (j = 0; j < i; j++) {
1784  cart2hough(m_cellInfo[i].x, m_cellInfo[i].y, m_cellInfo[j].x, m_cellInfo[j].y, rxy, axy);
1785  cart2hough(m_cellInfo[i].z, m_cellInfo[i].y, m_cellInfo[j].z, m_cellInfo[j].y, rzy, azy);
1786 
1787  if (rxy < nminrxy || rxy > nmaxrxy || axy < nminaxy || axy > nmaxaxy || rzy < nminrzy
1788  || rzy > nmaxrzy || azy < nminazy || azy > nmaxazy) continue;
1789 
1790  weight = m_cellInfo[i].ev * m_cellInfo[j].ev;
1791 
1792  arxy += rxy * weight;
1793  aaxy += axy * weight;
1794  arzy += rzy * weight;
1795  aazy += azy * weight;
1796  aw += weight;
1797  }
1798  }
1799 
1800  const double inv_aw = 1. / aw;
1801  rxy = arxy * inv_aw;
1802  axy = aaxy * inv_aw;
1803  rzy = arzy * inv_aw;
1804  azy = aazy * inv_aw;
1805 }
1806 
1807 int TileMuonFitter::houghTrack() {
1808 
1809  ATH_MSG_DEBUG( "Fitting with Hough method" << m_nCells << " cells." );
1810 
1811  // Build vector with cell Infos --------------------------------------------
1812  bool compute = true;
1813 
1814  unsigned int NPoints = buildCellInfoVector();
1815  if (NPoints > MAX_NCELLS) compute = false;
1816 
1817  // loop to find tracks -----------------------------------------------------
1818 
1819  while (compute) {
1820 
1821  // Get direction with maximum number of cells --------------------------
1822 
1823  unsigned int index1, index2;
1824  if (!guessTrack(index1, index2)) break;
1825 
1826  // select cells which belong to the track ------------------------------
1827 
1828  float p[3] = { m_cellInfo[index1].x, m_cellInfo[index1].y, m_cellInfo[index1].z };
1829  float w[3];
1830  points2dir(m_cellInfo[index1], m_cellInfo[index2], w);
1831 
1832  float en = selectCells(p, w);
1833  if (en < MIN_ENERGY_MEV) continue;
1834 
1835  // compute reference direction -----------------------------------------
1836 
1837  double rxy, axy, rzy, azy;
1838  cart2hough(m_cellInfo[index1].x, m_cellInfo[index1].y, m_cellInfo[index2].x, m_cellInfo[index2].y, rxy, axy);
1839  cart2hough(m_cellInfo[index1].z, m_cellInfo[index1].y, m_cellInfo[index2].z, m_cellInfo[index2].y, rzy, azy);
1840 
1841  // compute Hough Transform ---------------------------------------------
1842 
1843  doHough(rxy, axy, rzy, azy);
1844 
1845  // compute track parameters --------------------------------------------
1846 
1847  double aa, bb, cc, dd;
1848 
1849  hough2cart(rxy, axy, SHIFT_X, bb, aa);
1850  hough2cart(rzy, azy, SHIFT_Z, dd, cc);
1851 
1852  addTrack(aa, bb, (aa - cc) / dd, bb / dd);
1853 
1854  ATH_MSG_DEBUG( "Minimum Value aa: " << aa
1855  << " bb: " << bb
1856  << " cc: " << (aa - cc) / dd
1857  << " dd " << bb / dd );
1858 
1859  m_fitMinimum.push_back(1.0);
1860  }
1861 
1862  // prepare comtime ---------------------------------------------------------
1863 
1864  std::vector<double> X;
1865  std::vector<double> Y;
1866  std::vector<double> Z;
1867 
1868  X.resize(m_nCells);
1869  Y.resize(m_nCells);
1870  Z.resize(m_nCells);
1871 
1872  for (int i = 0; i < m_nCells; i++) {
1873  X[i] = m_cellPosition[i].getX();
1874  Y[i] = m_cellPosition[i].getY();
1875  Z[i] = m_cellPosition[i].getZ();
1876  }
1877 
1878  m_theTrack = new TileMuonTrackDistance(X, Y, Z, m_cellWeight);
1879  m_theTrack->SetWeighted(m_doWeighted);
1880 
1881  if (m_linePar.size() > 0) return 1;
1882  else return 0;
1883 }
1884 
1885 void TileMuonFitter::addTrack(double aa, double bb, double cc, double dd) {
1886  std::vector<double> par;
1887  par.resize(4);
1888 
1889  par[0] = aa;
1890  par[1] = bb;
1891  par[2] = cc;
1892  par[3] = dd;
1893 
1894  m_linePar.push_back(par);
1895 }
CellInfo::z
float z
Definition: TileMuonFitter.h:66
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
TileMuonFitter::m_caloIndex
static const CaloCell_ID::SUBCALO m_caloIndex
Definition: TileMuonFitter.h:262
TileDetDescrManager::tile_descriptors_begin
tile_descr_const_iterator tile_descriptors_begin() const
Definition: TileDetDescrManager.h:101
BIN_RES_AZY
#define BIN_RES_AZY
Definition: TileMuonFitter.cxx:78
TileCell
Definition: TileCell.h:57
ROOT::Minuit2::TileMuonTrackDistance::GetMeanX
double GetMeanX()
X mean getter.
Definition: TileMuonTrackDistance.h:97
TileMuonFitter::m_tileDD_zLB
std::vector< double > m_tileDD_zLB
Z bounds of LB, loaded from Detector Description.
Definition: TileMuonFitter.h:225
TileMuonFitter::setEventDefaults
void setEventDefaults()
Reset variables.
Definition: TileMuonFitter.cxx:297
beamspotman.r
def r
Definition: beamspotman.py:676
RunTileCalibRec.cells
cells
Definition: RunTileCalibRec.py:280
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
CaloCellContainer::indexLastCellCalo
int indexLastCellCalo(const CaloCell_ID::SUBCALO caloNum) const
index of last cell of given calorimeter (-2 if none) Note that it is normally more efficient to use i...
Definition: CaloCellContainer.cxx:141
ReadCellNoiseFromCoolCompare.s1
s1
Definition: ReadCellNoiseFromCoolCompare.py:378
Tile_Base_ID::cell_hash_max
size_type cell_hash_max(void) const
Definition: Tile_Base_ID.cxx:1313
TileMuonFitter::whichModule
int whichModule(CLHEP::Hep3Vector tempvec)
Returns module index for TVector3 input.
Definition: TileMuonFitter.cxx:1369
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
TileMuonFitter::buildCells
void buildCells()
Creates an internal cell container (just vectors) from the input CaloCellContainer.
Definition: TileMuonFitter.cxx:319
TileCosmicMuon::SetSegmentPath
void SetSegmentPath(const std::vector< double > &path)
Definition: TileCosmicMuon.h:69
ReadCellNoiseFromCool.cell
cell
Definition: ReadCellNoiseFromCool.py:53
Merge_LAr_blobs.filled
list filled
Definition: Merge_LAr_blobs.py:43
TileMuonFitter::m_cellTime
std::vector< double > m_cellTime
Selected cell's time.
Definition: TileMuonFitter.h:238
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
TileCosmicMuon
Class containing detailed results from TileMuonFitter.
Definition: TileCosmicMuon.h:35
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TileMuonFitter::m_zeroCrossingTime
std::vector< double > m_zeroCrossingTime
Time at y=0.
Definition: TileMuonFitter.h:249
MIN_ENERGY_MEV
#define MIN_ENERGY_MEV
Definition: TileMuonFitter.cxx:82
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:557
CaloDetDescrElement::y
float y() const
cell y
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:365
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
TileMuonFitter::trackIntersection
void trackIntersection(std::vector< double > &ltop, std::vector< double > &lbot, int index)
Calculates length of track intersection with TileCal (by sampling).
Definition: TileMuonFitter.cxx:827
TileMuonFitter::m_maxBottomIndex
int m_maxBottomIndex
Definition: TileMuonFitter.h:258
perp
Scalar perp() const
perp method - perpenticular length
Definition: AmgMatrixBasePlugin.h:44
TileMuonFitter::houghTrack
int houghTrack()
Fits a straight track to the cells centers, using a Hough Transform algorithm.
Definition: TileMuonFitter.cxx:1807
TileMuonFitter::buildTileCosmicMuon
void buildTileCosmicMuon(int fitok)
Selects between the two next methods.
Definition: TileMuonFitter.cxx:668
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
index
Definition: index.py:1
TileMuonFitter::eventSelection
bool eventSelection()
Checks if there are good cells on the top and bottom modules.
Definition: TileMuonFitter.cxx:378
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
TileMuonFitter::m_cellContainerKey
SG::ReadHandleKey< CaloCellContainer > m_cellContainerKey
Definition: TileMuonFitter.h:264
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
TileMuonFitter::buildComTimeAtYequal0
void buildComTimeAtYequal0(int fitok)
Creates output ComTime object in StoreGate.
Definition: TileMuonFitter.cxx:1447
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
TileCosmicMuon::SetPositionX
void SetPositionX(double posx)
Definition: TileCosmicMuon.h:51
TileMuonFitter::checkLBr
bool checkLBr(double x1, uint8_t s1)
Checks if x1 is within LB r coordinate bounds for sampling s1.
Definition: TileMuonFitter.cxx:1327
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell. Provides a generic interface for al...
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
TileMuonFitter::cart2hough
void cart2hough(float x1, float y1, float x2, float y2, double &raio, double &angu)
Definition: TileMuonFitter.cxx:1577
TileMuonFitter::buildComTimeAtZequal0
void buildComTimeAtZequal0(int fitok)
Creates output ComTime object in StoreGate.
Definition: TileMuonFitter.cxx:1512
TileMuonFitter::m_tileDD_radiusEB
std::vector< double > m_tileDD_radiusEB
Radial bounds of the 3 samplings in EB, loaded from Detector Description.
Definition: TileMuonFitter.h:219
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
TileMuonFitter::calculateTimeAtZequal0
void calculateTimeAtZequal0()
Extrapolates cell time to z=0.
Definition: TileMuonFitter.cxx:605
CaloID.h
TileMuonFitter::calculateTime
void calculateTime()
Calculates time in reference plane.
Definition: TileMuonFitter.cxx:530
Tile_Base_ID::sample
int sample(const Identifier &id) const
Definition: Tile_Base_ID.cxx:171
TileMuonFitter::m_cellPosition
std::vector< CLHEP::Hep3Vector > m_cellPosition
Position of selected cell's center.
Definition: TileMuonFitter.h:232
yodamerge_tmp.npoints
npoints
Definition: yodamerge_tmp.py:250
AthCommonMsg< Algorithm >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
TileMuonFitter::execute
virtual StatusCode execute()
Definition: TileMuonFitter.cxx:253
dq_defect_virtual_defect_validation.d1
d1
Definition: dq_defect_virtual_defect_validation.py:79
TileCosmicMuon::SetSegmentSampling
void SetSegmentSampling(const std::vector< int > &sampling)
Definition: TileCosmicMuon.h:75
ROOT::Minuit2
Definition: TileMuonTrackDistance.cxx:26
TileMuonFitter::dist2line
float dist2line(CellInfo &ci, float *pos, float *w)
Definition: TileMuonFitter.cxx:1617
TileMuonFitter::m_linePar
std::vector< std::vector< double > > m_linePar
Vector with the fitted four track parameters.
Definition: TileMuonFitter.h:245
TileMuonFitter::buildCellInfoVector
unsigned int buildCellInfoVector()
Definition: TileMuonFitter.cxx:1698
TileMuonFitter::m_maxTopIndex
int m_maxTopIndex
Definition: TileMuonFitter.h:259
TileCosmicMuon::SetEnergyBottom
void SetEnergyBottom(const std::vector< double > &energy)
Definition: TileCosmicMuon.h:65
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
BIN_RES_AXY
#define BIN_RES_AXY
Definition: TileMuonFitter.cxx:77
SG::VarHandleKey::key
const std::string & key() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:141
TileMuonFitter::calculateTimeAtYequal0
void calculateTimeAtYequal0()
Extrapolates cell time to y=0.
Definition: TileMuonFitter.cxx:542
BIN_RES_RXY
#define BIN_RES_RXY
Definition: TileMuonFitter.cxx:75
TileCosmicMuon::SetEnergyTop
void SetEnergyTop(const std::vector< double > &energy)
Definition: TileCosmicMuon.h:63
TileMuonFitter::m_fitMinimum
std::vector< double > m_fitMinimum
Chi-square minumum.
Definition: TileMuonFitter.h:247
x
#define x
TileDetDescrManager::tile_descriptors_end
tile_descr_const_iterator tile_descriptors_end() const
Definition: TileDetDescrManager.h:104
TileDetDescriptor::dz
float dz(unsigned int samp) const
Definition: TileDetDescriptor.h:212
TileDetDescriptor::sign_eta
int sign_eta() const
Definition: TileDetDescriptor.h:224
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
TileMuonFitter::m_cellInfo
std::vector< CellInfo > m_cellInfo
Definition: TileMuonFitter.h:177
pi
#define pi
Definition: TileMuonFitter.cxx:65
TileMuonFitter::m_cellEnergy
std::vector< double > m_cellEnergy
Selected cell's energy.
Definition: TileMuonFitter.h:234
MAX_NCELLS
#define MAX_NCELLS
Definition: TileMuonFitter.cxx:79
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
TileDetDescrManager.h
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
TileMuonFitter::whichLBr
int whichLBr(double x1)
Returns sampling index if x1 is within LB r coordinate bounds.
Definition: TileMuonFitter.cxx:1359
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:95
TileDetDescriptor::n_eta
int n_eta(unsigned int samp) const
Definition: TileDetDescriptor.h:159
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
TileCosmicMuon::SetPositionY
void SetPositionY(double posy)
Definition: TileCosmicMuon.h:52
TileMuonFitter::m_cellWeight
std::vector< double > m_cellWeight
Selected cell's weight for fit.
Definition: TileMuonFitter.h:236
TileDetDescriptor::dr
float dr(unsigned int samp) const
Definition: TileDetDescriptor.h:198
CellInfo::x
float x
Definition: TileMuonFitter.h:66
tools.zlumi_mc_cf.correction
def correction(mu, runmode, campaign, run=None)
Definition: zlumi_mc_cf.py:4
TileCosmicMuon.h
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
TileCell::timeDiff
float timeDiff(void) const
get time diff for two PMTs (data member)
Definition: TileCell.h:190
TileMuonFitter::checkEBr
bool checkEBr(double x1, uint8_t s1)
Checks if x1 is within EB r coordinate bounds for sampling s1.
Definition: TileMuonFitter.cxx:1338
WriteHandle.h
Handle class for recording to StoreGate.
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
ROOT::Minuit2::TileMuonTrackDistance::SetWeighted
void SetWeighted(bool weight)
Set weight/no-weight flag.
Definition: TileMuonTrackDistance.h:110
TileMuonFitter::m_doHoughTransform
bool m_doHoughTransform
Flag to use Hough Transform instead of Fit.
Definition: TileMuonFitter.h:215
SHIFT_Z
#define SHIFT_Z
Definition: TileMuonFitter.cxx:84
CaloCellContainer::indexFirstCellCalo
int indexFirstCellCalo(const CaloCell_ID::SUBCALO caloNum) const
index of first cell of given calorimeter (-1 if none).
Definition: CaloCellContainer.cxx:137
TileDetDescriptor::rcenter
float rcenter(unsigned int samp) const
Definition: TileDetDescriptor.h:191
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
MIN_NCELLS
#define MIN_NCELLS
Definition: TileMuonFitter.cxx:80
TileMuonFitter::checkEBAz
bool checkEBAz(double x1)
Checks if x1 is within EBA z coordinate bounds.
Definition: TileMuonFitter.cxx:1289
TileHWID.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
Trk::index1
@ index1
Definition: BoundarySurfaceFace.h:48
maskDeadModules.mod
mod
Definition: maskDeadModules.py:36
TileMuonFitter::~TileMuonFitter
virtual ~TileMuonFitter()
Definition: TileMuonFitter.cxx:125
TileMuonFitter::TileMuonFitter
TileMuonFitter(const std::string &name, ISvcLocator *pSvcLocator)
Constructor.
Definition: TileMuonFitter.cxx:95
xAOD::uint16_t
setWord1 uint16_t
Definition: eFexEMRoI_v1.cxx:93
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:116
parseMapping.v0
def v0
Definition: parseMapping.py:149
lumiFormat.i
int i
Definition: lumiFormat.py:85
TileMuonFitter::checkLBz
bool checkLBz(double x1)
Checks if x1 is within LB z coordinate bounds.
Definition: TileMuonFitter.cxx:1274
TileCosmicMuon::SetDirectionPhi
void SetDirectionPhi(double phi)
Definition: TileCosmicMuon.h:54
TileMuonFitter::guessTrack
bool guessTrack(unsigned int &index1, unsigned int &index2)
Definition: TileMuonFitter.cxx:1673
z
#define z
TileDetDescriptor::zcenter
float zcenter(unsigned int samp) const
Definition: TileDetDescriptor.h:205
TileDetDescriptor::n_samp
int n_samp() const
Definition: TileDetDescriptor.h:219
TileMuonFitter::CntCells
unsigned int CntCells(unsigned int index1, unsigned int index2, double &skew)
Definition: TileMuonFitter.cxx:1643
beamspotman.n
n
Definition: beamspotman.py:731
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TileDetDescriptor
Definition: TileDetDescriptor.h:46
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
CaloCell_Base_ID::SUBCALO
SUBCALO
enumeration of sub calorimeters
Definition: CaloCell_Base_ID.h:46
SHIFT_X
#define SHIFT_X
Definition: TileMuonFitter.cxx:83
TileCosmicMuon::SetPositionZ
void SetPositionZ(double posz)
Definition: TileCosmicMuon.h:53
CaloCellContainer::nCellsCalo
int nCellsCalo(const CaloCell_ID::SUBCALO caloNum) const
get number of cels of given calorimeter
Definition: CaloCellContainer.cxx:145
TileCosmicMuon::SetSegmentPartition
void SetSegmentPartition(const std::vector< int > &partition)
Definition: TileCosmicMuon.h:71
TRT::Track::d0
@ d0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:62
TileMuonFitter::m_cellHash
std::vector< IdentifierHash > m_cellHash
Selected cell's identifier hash.
Definition: TileMuonFitter.h:242
TileMuonFitter::m_theTrack
ROOT::Minuit2::TileMuonTrackDistance * m_theTrack
Auxiliary class representing the function to be minimized - weighted sum of squares of orthogonal dis...
Definition: TileMuonFitter.h:201
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
TileMuonFitter::m_cellDeltaTime
std::vector< double > m_cellDeltaTime
Selected cell's time difference between two PMTs.
Definition: TileMuonFitter.h:240
CellInfo
Definition: TileMuonFitter.h:63
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
CaloCell_Base_ID::TILE
@ TILE
Definition: CaloCell_Base_ID.h:46
TileCosmicMuon::SetSegmentModule
void SetSegmentModule(const std::vector< int > &module)
Definition: TileCosmicMuon.h:73
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
ROOT::Minuit2::TileMuonTrackDistance::GetMeanZ
double GetMeanZ()
Z mean getter.
Definition: TileMuonTrackDistance.h:105
TileMuonFitter::buildTileCosmicMuonAtZequal0
void buildTileCosmicMuonAtZequal0(int fitok)
Creates output TileCosmicMuon object in StoreGate.
Definition: TileMuonFitter.cxx:754
TileMuonFitter::checkEBz
bool checkEBz(double x1)
Checks if x1 is within EB z coordinate bounds.
Definition: TileMuonFitter.cxx:1283
AthAlgorithm
Definition: AthAlgorithm.h:47
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
ROOT::Minuit2::TileMuonTrackDistance::Means
void Means()
Calculates means (weighted or not) of cell positions.
Definition: TileMuonTrackDistance.cxx:28
TileMuonFitter::trackSegmentIntersection
void trackSegmentIntersection(std::vector< double > &segPath, std::vector< int > &segPartition, std::vector< int > &segModule, std::vector< int > &segSampling, int index)
Calculates length of track intersection with TileCal (by sampling and module).
Definition: TileMuonFitter.cxx:1053
TileDetDescrManager::get_cell_element
CaloDetDescrElement * get_cell_element(unsigned int cell_hash) const
Definition: TileDetDescrManager.h:156
CxxUtils::set
constexpr std::enable_if_t< is_bitmask_v< E >, E & > set(E &lhs, E rhs)
Convenience function to set bits in a class enum bitmask.
Definition: bitmask.h:232
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
ROOT::Minuit2::TileMuonTrackDistance::GetMeanY
double GetMeanY()
Y mean getter.
Definition: TileMuonTrackDistance.h:101
TileMuonFitter::whichEBr
int whichEBr(double x1)
Returns sampling index if x1 is within EB r coordinate bounds.
Definition: TileMuonFitter.cxx:1349
selection
const std::string selection
Definition: fbtTestBasics.cxx:74
Trk::index2
@ index2
Definition: BoundarySurfaceFace.h:49
HWIdentifier.h
CaloDetDescrElement::volume
float volume() const
cell volume
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:381
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
TileMuonFitter::checkEBCz
bool checkEBCz(double x1)
Checks if x1 is within EBC z coordinate bounds.
Definition: TileMuonFitter.cxx:1298
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
TileMuonFitter::buildTileCosmicMuonAtYequal0
void buildTileCosmicMuonAtYequal0(int fitok)
Creates output TileCosmicMuon object in StoreGate.
Definition: TileMuonFitter.cxx:681
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
PlotCalibFromCool.en
en
Definition: PlotCalibFromCool.py:399
TileCosmicMuon::SetPathBottom
void SetPathBottom(const std::vector< double > &path)
Definition: TileCosmicMuon.h:61
TileMuonFitter::doHough
void doHough(double &rxy, double &axy, double &rzy, double &azy)
Definition: TileMuonFitter.cxx:1759
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
SG::VarHandleBase::key
virtual const std::string & key() const override final
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleBase.cxx:64
TileCosmicMuon::SetFitNCells
void SetFitNCells(int ncells)
Definition: TileCosmicMuon.h:57
library_scraper.dd
list dd
Definition: library_scraper.py:46
TileMuonFitter::energyInTrack
void energyInTrack(std::vector< double > &etop, std::vector< double > &ebot, std::vector< IdentifierHash > &cells, int index)
Sums up energy in TileCal cells close to the track (by sampling).
Definition: TileMuonFitter.cxx:1382
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
CaloDetDescrElement::x
float x() const
cell x
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:363
TileMuonFitter::m_tileID
const TileID * m_tileID
Definition: TileMuonFitter.h:194
TileMuonFitter::m_tileDD_zEBC
std::vector< double > m_tileDD_zEBC
Z bounds of EBC, loaded from Detector Description.
Definition: TileMuonFitter.h:223
TileMuonFitter::m_doWeighted
bool m_doWeighted
Flag to weigh or not the chi-square with an energy parameter.
Definition: TileMuonFitter.h:208
ReadCellNoiseFromCoolCompare.v2
v2
Definition: ReadCellNoiseFromCoolCompare.py:364
TileMuonFitter::m_nCells
int m_nCells
Number of cells selected for fit.
Definition: TileMuonFitter.h:228
ROOT::Minuit2::TileMuonTrackDistance::ClosestPoint
CLHEP::Hep3Vector ClosestPoint(CLHEP::Hep3Vector *dataPoint, const std::vector< double > &par) const
Returns point in track defined by par closest to dataPoint.
Definition: TileMuonTrackDistance.cxx:93
TileMuonFitter::buildComTime
void buildComTime(int fitok)
Selects between the two next methods.
Definition: TileMuonFitter.cxx:1435
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:73
DeMoScan.index
string index
Definition: DeMoScan.py:364
TileMuonFitter::m_beamType
std::string m_beamType
Flag to indicate: cosmics, singlebeam or collisions.
Definition: TileMuonFitter.h:213
a
TList * a
Definition: liststreamerinfos.cxx:10
TileMuonFitter::isHaloMuon
bool isHaloMuon(double azy)
Definition: TileMuonFitter.cxx:1753
y
#define y
BIN_RES_RZY
#define BIN_RES_RZY
Definition: TileMuonFitter.cxx:76
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
TileMuonFitter::addTrack
void addTrack(double aa, double bb, double cc, double dd)
Definition: TileMuonFitter.cxx:1885
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
trigbs_pickEvents.cnt
cnt
Definition: trigbs_pickEvents.py:71
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TileMuonFitter::points2dir
void points2dir(CellInfo &ci1, CellInfo &ci2, float *w)
Definition: TileMuonFitter.cxx:1631
CaloDetDescrElement::z
float z() const
cell z
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:367
dq_defect_virtual_defect_validation.d2
d2
Definition: dq_defect_virtual_defect_validation.py:81
DeMoScan.first
bool first
Definition: DeMoScan.py:536
TileMuonFitter::initialize
virtual StatusCode initialize()
Definition: TileMuonFitter.cxx:131
DEBUG
#define DEBUG
Definition: page_access.h:11
F
#define F(x, y, z)
Definition: MD5.cxx:112
CellInfo::y
float y
Definition: TileMuonFitter.h:66
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
TileMuonFitter::m_doDensity
bool m_doDensity
Flag defining the energy weighting parameter: energy density or plain energy.
Definition: TileMuonFitter.h:211
TileMuonFitter::m_tileDD_zEBA
std::vector< double > m_tileDD_zEBA
Z bounds of EBA, loaded from Detector Description.
Definition: TileMuonFitter.h:221
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
TileCosmicMuon::SetDirectionTheta
void SetDirectionTheta(double theta)
Definition: TileCosmicMuon.h:55
TileMuonFitter::m_tileDD_radiusLB
std::vector< double > m_tileDD_radiusLB
Radial bounds of the 3 samplings in LB, loaded from Detector Description.
Definition: TileMuonFitter.h:217
TileCosmicMuon::SetTrackCellHash
void SetTrackCellHash(const std::vector< IdentifierHash > &cells)
Definition: TileCosmicMuon.h:67
TileMuonFitter::finalize
virtual StatusCode finalize()
Definition: TileMuonFitter.cxx:243
TileCosmicMuon::SetFitQuality
void SetFitQuality(double quality)
Definition: TileCosmicMuon.h:56
CaloDetDescrElement::eta
float eta() const
cell eta
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:344
TileCosmicMuon::SetPathTop
void SetPathTop(const std::vector< double > &path)
Definition: TileCosmicMuon.h:59
TileMuonFitter::m_tileHWID
const TileHWID * m_tileHWID
Definition: TileMuonFitter.h:195
Tile_Base_ID::cell_hash
IdentifierHash cell_hash(const Identifier &cell_id) const
fast conversion from ID to hash for cells
Definition: Tile_Base_ID.cxx:1030
rr
const boost::regex rr(r_r)
SelectAllObject.h
Tile_Base_ID::section
int section(const Identifier &id) const
Definition: Tile_Base_ID.cxx:147
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
TileMuonFitter::selectCells
float selectCells(float *p, float *w)
Definition: TileMuonFitter.cxx:1732
ReadHandle.h
Handle class for reading from StoreGate.
CaloDetDescrElement::r
float r() const
cell r
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:348
TileMuonFitter::m_eThreshold
double m_eThreshold
Cell energy threshold.
Definition: TileMuonFitter.h:204
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
section
void section(const std::string &sec)
Definition: TestTriggerMenuAccess.cxx:22
TileMuonFitter::hough2cart
void hough2cart(double r, double a, double offset, double &aa, double &bb)
Definition: TileMuonFitter.cxx:1603
TRTCalib_cfilter.p3
p3
Definition: TRTCalib_cfilter.py:132
Tile_Base_ID::cell_id
Identifier cell_id(const Identifier &any_id) const
Definition: Tile_Base_ID.cxx:581
TileMuonFitter::m_cosmicMuonContainerKey
SG::WriteHandleKey< TileCosmicMuonContainer > m_cosmicMuonContainerKey
Definition: TileMuonFitter.h:267
ROOT::Minuit2::TileMuonTrackDistance
Auxiliary to TileMuonFitter. Chi-square function for Minuit.
Definition: TileMuonTrackDistance.h:70
test_AnalysisBaseEventLoopJob.aa
aa
Definition: test_AnalysisBaseEventLoopJob.py:37
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
TileMuonFitter::m_reg1to2
bool m_reg1to2
Definition: TileMuonFitter.h:260
TileMuonFitter.h
TileMuonFitter::m_minimumCells
int m_minimumCells
Minimum number of cells needed for fit.
Definition: TileMuonFitter.h:230
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
TileMuonFitter::m_tileMgr
const TileDetDescrManager * m_tileMgr
Definition: TileMuonFitter.h:196
TileMuonFitter::m_caloCells
const CaloCellContainer * m_caloCells
Definition: TileMuonFitter.h:197
TileMuonFitter::m_deltaTimeCut
double m_deltaTimeCut
Cell Delta Time cut.
Definition: TileMuonFitter.h:206
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129
fitman.k
k
Definition: fitman.py:528
python.handimod.cc
int cc
Definition: handimod.py:523
ymax
double ymax
Definition: listroot.cxx:64
TileMuonFitter::fitTrack
int fitTrack()
Fits a straight track to the cell centers, using the auxiliary class TileMuonTrackDistance.
Definition: TileMuonFitter.cxx:435
TileMuonFitter::m_comTimeKey
SG::WriteHandleKey< ComTime > m_comTimeKey
Definition: TileMuonFitter.h:272
TileCosmicMuon::SetTime
void SetTime(double time)
Setters.
Definition: TileCosmicMuon.h:50
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Mon Apr 28 2025 21:20:31 for ATLAS Offline Software by   doxygen 1.8.18