← Back to branch summary

~mice-lcr/maus/end-of-emr-run

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/SciFi/SciFiGeometryHelper.cc

  • Committer: Durga Rajaram
  • Date: 2013-08-27 04:36:50 UTC
  • mfrom: (659.1.73 rc)
  • Revision ID: durga@fnal.gov-20130827043650-me0hgsbzlzikdoik
Tags: MAUS-v0.7, MAUS-v0.7.0
MAUS-v0.7.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/SciFi/SciFiGeometryHelper.cc
 
1
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 
2
 *
 
3
 * MAUS is free software: you can redistribute it and/or modify
 
4
 * it under the terms of the GNU General Public License as published by
 
5
 * the Free Software Foundation, either version 3 of the License, or
 
6
 * (at your option) any later version.
 
7
 *
 
8
 * MAUS is distributed in the hope that it will be useful,
 
9
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
 * GNU General Public License for more details.
 
12
 *
 
13
 * You should have received a copy of the GNU General Public License
 
14
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 
15
 *
 
16
 */
 
17
 
 
18
#include "src/common_cpp/Recon/SciFi/SciFiGeometryHelper.hh"
 
19
 
 
20
namespace MAUS {
 
21
 
 
22
SciFiGeometryHelper::SciFiGeometryHelper() {}
 
23
 
 
24
SciFiGeometryHelper::SciFiGeometryHelper(const std::vector<const MiceModule*> &modules)
 
25
                                        : _modules(modules) {}
 
26
 
 
27
SciFiGeometryHelper::~SciFiGeometryHelper() {}
 
28
 
 
29
void SciFiGeometryHelper::Build() {
 
30
  // Iterate over existing modules, adding planes to the map.
 
31
  std::vector<const MiceModule*>::iterator iter;
 
32
  for ( iter = _modules.begin(); iter != _modules.end(); iter++ ) {
 
33
    const MiceModule* module = (*iter);
 
34
    if ( module->propertyExists("Tracker", "int") &&
 
35
         module->propertyExists("Station", "int") &&
 
36
         module->propertyExists("Plane", "int") ) {
 
37
      int tracker_n = module->propertyInt("Tracker");
 
38
      int station_n = module->propertyInt("Station");
 
39
      int plane_n   = module->propertyInt("Plane");
 
40
 
 
41
      double pitch        = module->propertyDouble("Pitch");
 
42
      double centralfibre = module->propertyDouble("CentralFibre");
 
43
      ThreeVector direction(0., 1., 0.);
 
44
 
 
45
      // G4RotationMatrix global_fibre_rotation = G4RotationMatrix(module->globalRotation());
 
46
      const MiceModule* plane = module->mother();
 
47
      G4RotationMatrix internal_fibre_rotation(module->relativeRotation(module->mother() // plane
 
48
                                               ->mother()));  // tracker/ station??
 
49
 
 
50
      direction     *= internal_fibre_rotation;
 
51
 
 
52
      // The plane rotation wrt to the solenoid. Identity matrix for tracker 1,
 
53
      // [ -1, 0, 0],[ 0, 1, 0],[ 0, 0, -1] for tracker 0 (180 degrees rot. around y).
 
54
      // const MiceModule* plane = module->mother();
 
55
      G4RotationMatrix plane_rotation(plane->relativeRotation(plane->mother()  // tracker
 
56
                                                              ->mother()));    // solenoid
 
57
 
 
58
      ThreeVector position  = clhep_to_root(module->globalPosition());
 
59
      ThreeVector reference = GetReferenceFramePosition(tracker_n);
 
60
 
 
61
      ThreeVector tracker_ref_frame_pos = position-reference;
 
62
      tracker_ref_frame_pos *= plane_rotation;
 
63
 
 
64
      SciFiPlaneGeometry this_plane;
 
65
      this_plane.Direction     = direction;
 
66
      this_plane.Position      = tracker_ref_frame_pos;
 
67
      this_plane.CentralFibre  = centralfibre;
 
68
      this_plane.Pitch         = pitch;
 
69
      int plane_id =  3*(station_n-1) + (plane_n+1);
 
70
      plane_id     = ( tracker_n == 0 ? -plane_id : plane_id );
 
71
      _geometry_map.insert(std::make_pair(plane_id, this_plane));
 
72
      _field_value[tracker_n] = FieldValue(reference, plane_rotation);
 
73
    }
 
74
  }
 
75
}
 
76
 
 
77
double SciFiGeometryHelper::FieldValue(ThreeVector global_position,
 
78
                                       G4RotationMatrix plane_rotation) {
 
79
  double EMfield[6]  = {0., 0., 0., 0., 0., 0.};
 
80
  double position[4] = {global_position.x(), global_position.y(), global_position.z(), 0.};
 
81
  BTFieldConstructor* field = Globals::GetMCFieldConstructor();
 
82
  field->GetElectroMagneticField()->GetFieldValue(position, EMfield);
 
83
  ThreeVector B_field(EMfield[0], EMfield[1], EMfield[2]);
 
84
  B_field *= plane_rotation;
 
85
  double Tracker_Bz = B_field.z();
 
86
  return Tracker_Bz;
 
87
}
 
88
 
 
89
const MiceModule* SciFiGeometryHelper::FindPlane(int tracker, int station, int plane) {
 
90
  const MiceModule* this_plane = NULL;
 
91
  for ( unsigned int j = 0; !this_plane && j < _modules.size(); ++j ) {
 
92
    // Find the right module
 
93
    if ( _modules[j]->propertyExists("Tracker", "int") &&
 
94
         _modules[j]->propertyExists("Station", "int") &&
 
95
         _modules[j]->propertyExists("Plane", "int")  &&
 
96
         _modules[j]->propertyInt("Tracker") ==
 
97
         tracker &&
 
98
         _modules[j]->propertyInt("Station") ==
 
99
         station &&
 
100
         _modules[j]->propertyInt("Plane") ==
 
101
         plane ) {
 
102
         // Save the module
 
103
      this_plane = _modules[j];
 
104
    }
 
105
  }
 
106
  if ( this_plane == NULL ) {
 
107
    throw(Squeal(Squeal::nonRecoverable,
 
108
    "Failed to find tracker plane.",
 
109
    "SciFiGeometryHelper::find_plane"));
 
110
  }
 
111
  return this_plane;
 
112
}
 
113
 
 
114
ThreeVector SciFiGeometryHelper::GetReferenceFramePosition(int tracker) {
 
115
  // Reference plane is plane 0, station 1 of current tracker.
 
116
  int station = 1;
 
117
  int plane   = 0;
 
118
  const MiceModule* reference_plane = NULL;
 
119
  reference_plane = FindPlane(tracker, station, plane);
 
120
 
 
121
  assert(reference_plane != NULL);
 
122
  ThreeVector reference_pos = clhep_to_root(reference_plane->globalPosition());
 
123
 
 
124
  return reference_pos;
 
125
}
 
126
 
 
127
void SciFiGeometryHelper::DumpPlanesInfo() {
 
128
  std::map<int, SciFiPlaneGeometry>::iterator plane;
 
129
  for ( plane = _geometry_map.begin(); plane != _geometry_map.end(); ++plane ) {
 
130
    Squeak::mout(Squeak::info) << "Plane ID: " << plane->first << "\n"
 
131
                               << "Direction: "<< plane->second.Direction << "\n"
 
132
                               << "Position: " << plane->second.Position << "\n"
 
133
                               << "CentralFibre: "<< plane->second.CentralFibre << "\n"
 
134
                               << "Pitch: "       << plane->second.Pitch << "\n";
 
135
  }
 
136
}
 
137
 
 
138
} // ~namespace MAUS