ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4RPSubsystem.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4RPSubsystem.cc
1 //____________________________________________________________________________..
2 //
3 #include "EICG4RPSubsystem.h"
4 
5 #include "EICG4RPDetector.h"
7 
8 #include <phparameter/PHParameters.h>
9 
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHIODataNode.h>
15 #include <phool/PHNode.h>
16 #include <phool/PHNodeIterator.h>
17 #include <phool/PHObject.h>
18 #include <phool/getClass.h>
19 
20 #include <TSystem.h>
21 
22 #include <cmath>
23 #include <cstdlib>
24 #include <fstream>
25 #include <iostream>
26 #include <sstream>
27 #include <utility>
28 
29 //_______________________________________________________________________
30 EICG4RPSubsystem::EICG4RPSubsystem(const std::string &name, const int lyr)
31  : PHG4DetectorSubsystem(name, lyr)
32  , m_Detector(nullptr)
33  , m_SteppingAction(nullptr)
34 {
35  // call base class method which will set up parameter infrastructure
36  // and call our SetDefaultParameters() method
38 }
39 //_______________________________________________________________________
41 {
42  // create detector
43  m_Detector = new EICG4RPDetector(this, topNode, GetParams(), Name(), GetLayer());
46 
47  if (GetParams()->get_int_param("active"))
48  {
49  PHNodeIterator iter(topNode);
50  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
51 
52  std::string nodename;
53  std::string nodenameVirtSheet;
54 
55  if (SuperDetector() != "NONE")
56  {
57  // create super detector subnodes
58  PHNodeIterator iter_dst(dstNode);
59  PHCompositeNode *superSubNode = dynamic_cast<PHCompositeNode *>(iter_dst.findFirst("PHCompositeNode", SuperDetector()));
60  if (!superSubNode)
61  {
62  superSubNode = new PHCompositeNode(SuperDetector());
63  dstNode->addNode(superSubNode);
64  }
65  dstNode = superSubNode;
66 
67  nodename = "G4HIT_" + SuperDetector();
68  }
69  else
70  {
71  nodename = "G4HIT_" + Name();
72  }
73 
74  nodenameVirtSheet = nodename + "_VirtSheet";
75 
76  PHG4HitContainer *rp_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
77  PHG4HitContainer *rp_hitsVirtSheet = findNode::getClass<PHG4HitContainer>(topNode, nodenameVirtSheet);
78 
79  if( ! rp_hits )
80  {
81  dstNode->addNode(new PHIODataNode<PHObject>(rp_hits = new PHG4HitContainer(nodename), nodename, "PHObject"));
82  }
83  if( ! rp_hitsVirtSheet )
84  {
85  dstNode->addNode(new PHIODataNode<PHObject>(rp_hits = new PHG4HitContainer(nodenameVirtSheet), nodenameVirtSheet, "PHObject"));
86  }
87 
88  rp_hits->AddLayer(GetLayer());
89 
90  auto *tmp = new EICG4RPSteppingAction(this, m_Detector, GetParams());
91  tmp->HitNodeName(nodename);
92  tmp->HitNodeNameVirt(nodenameVirtSheet);
94  }
95  else if (GetParams()->get_int_param("blackhole"))
96  {
98  }
99  if (m_SteppingAction)
100  {
101  (dynamic_cast<EICG4RPSteppingAction *>(m_SteppingAction))->SaveAllHits(m_SaveAllHitsFlag);
102  }
103  return 0;
104 }
105 //_______________________________________________________________________
107 {
108  // pass top node to stepping action so that it gets
109  // relevant nodes needed internally
110  if (m_SteppingAction)
111  {
113  }
114  return 0;
115 }
116 //_______________________________________________________________________
117 void EICG4RPSubsystem::Print(const std::string &what) const
118 {
119  if (m_Detector)
120  {
121  m_Detector->Print(what);
122  }
123  return;
124 }
125 
126 //_______________________________________________________________________
128 {
129  return m_Detector;
130 }
131 
132 //_______________________________________________________________________
134 {
135 
136  std::string filename;
137  std::string calibroot = std::string(getenv("CALIBRATIONROOT"));
138 
139  if ( calibroot.empty() )
140  {
141  std::cout << PHWHERE << "no CALIBRATIONROOT environment variable" << std::endl;
142  gSystem->Exit(1);
143  }
144 
145  filename = calibroot + "/RomanPots/RP_parameters_IP6.dat";
146 
147  set_default_string_param("parameter_file", filename);
148  set_default_double_param("FFenclosure_center", 2500);
149  set_default_int_param("layerNumber", 1);
150  set_default_double_param("detid", 0.); //detector id
151  set_default_int_param("lightyield", 0);
152  set_default_int_param("use_g4steps", 0);
153  set_default_double_param("tmin", NAN);
154  set_default_double_param("tmax", NAN);
155 
156  set_default_double_param("Number_layers" , 0);
157  set_default_double_param("Sensor_size" , 0);
158 
159  // maximum of 4 layers for the IP8 config, only layers 1 and 2 used for IP6
160  for( int l=1; l <= 4; l++ ) {
161  set_default_double_param(Form("Layer%i_pos_x", l) , 0);
162  set_default_double_param(Form("Layer%i_pos_y", l) , 0);
163  set_default_double_param(Form("Layer%i_pos_z", l) , 0);
164  set_default_double_param(Form("Layer%i_size_x", l) , 0);
165  set_default_double_param(Form("Layer%i_size_y", l) , 0);
166  set_default_double_param(Form("Layer%i_beamHoleHW_x", l) , 0);
167  set_default_double_param(Form("Layer%i_beamHoleHW_y", l) , 0);
168  set_default_double_param(Form("Layer%i_Si_size_z", l) , 0);
169  set_default_double_param(Form("Layer%i_Cu_size_z", l) , 0);
170  set_default_double_param(Form("Layer%i_rot_y", l) , 0);
171  }
172 
173 }
174 
175 //_______________________________________________________________________
177 {
178  set_string_param("parameter_file", filename );
179 
180  std::ifstream infile;
181  std::string line;
182 
183  std::string paramFile = filename;
184  infile.open( paramFile );
185 
186  if(!infile.is_open()) {
187  std::cout << PHWHERE << "ERROR in EICG4RPSubsystem: Failed to open parameter file " << paramFile << std::endl;
188  gSystem->Exit(1);
189  }
190 
191  while( std::getline(infile, line) ) {
192 
193  std::string name;
194  double value;
195 
196  std::istringstream iss( line );
197 
198  // skip comment lines
199  if( line.find("#") != std::string::npos ) { continue; }
200  if( line.empty() == true ) { continue; }
201 
202  // grab the line
203  if( !(iss >> name >> value) ) {
204  std::cout << PHWHERE << " RP parameters: could not decode " << line << std::endl;
205  gSystem->Exit(1);
206  }
207 
208  // simply set all parameters NOT related to the beam hole
209  if( name.find("beamHoleHW") == std::string::npos ) {
210  set_double_param( name, value );
211  }
212  else { // parse beam hole size based on chosen beam energy and config
213 
214  std::string prefix = name.substr( 0, name.find("__") ); // LayerN_beamHoleHW_x(or)y
215  std::string expected_string = prefix;
216 
217  if( m_beamProfile.find("eA") != std::string::npos ) { // Heavy ions
218  expected_string += "__eA";
219 
220  double ionE_diff_1 = fabs( m_ionE - 110 );
221  double ionE_diff_2 = fabs( m_ionE - 41 );
222  double elecE_diff_1 = fabs( m_elecE - 18 );
223  double elecE_diff_2 = fabs( m_elecE - 10 );
224  double elecE_diff_3 = fabs( m_elecE - 5 );
225 
226  if( ionE_diff_2 < ionE_diff_1 ) {
227  expected_string += "_41x5";
228  }
229  else {
230  if( elecE_diff_3 < elecE_diff_2 ) { expected_string += "_110x5"; }
231  else if( elecE_diff_2 < elecE_diff_1 ) { expected_string += "_110x10"; }
232  else { expected_string += "_110x18"; }
233  }
234  }
235  else { // protons
236  if( m_beamProfile.find("high-acceptance") != std::string::npos ) {
237  expected_string += "__epHA";
238  }
239  else {
240  expected_string += "__epHD";
241  }
242 
243  double ionE_diff_1 = fabs( m_ionE - 275 );
244  double ionE_diff_2 = fabs( m_ionE - 100 );
245  double ionE_diff_3 = fabs( m_ionE - 41 );
246  double elecE_diff_1 = fabs( m_elecE - 18 );
247  double elecE_diff_2 = fabs( m_elecE - 10 );
248  double elecE_diff_3 = fabs( m_elecE - 5 );
249 
250  if( ionE_diff_3 < ionE_diff_2 ) { // 41 GeV protons
251  expected_string += "_41x5";
252  }
253  else if( ionE_diff_2 < ionE_diff_1 ) { // 100 GeV protons
254  if( elecE_diff_3 < elecE_diff_2 ) { expected_string += "_100x5"; }
255  else { expected_string += "_100x10"; }
256  }
257  else {
258  if( elecE_diff_2 < elecE_diff_1 ) { expected_string += "_275x10"; }
259  else { expected_string += "_275x18"; }
260  }
261  }
262 
263  // skip undesired beam hole configuration
264  if( name.compare( expected_string ) != 0 ) { continue; }
265 
266  set_double_param( prefix, value);
267  }
268 
269  }
270 
271  infile.close();
272 
273 }