FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
run_tutorial4_createParameterFile.C
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 
9 // forward declaration
10 void writeVectorToFile(std::vector<Double_t>& vec, ofstream* myfile);
11 
12 void run_tutorial4_createParameterFile()
13 {
14 
15  // Create shift in x, y, z when true
16  Bool_t ShiftX = kTRUE;
17  Bool_t ShiftY = kTRUE;
18  Bool_t ShiftZ = kFALSE;
19 
20  // Create a rotation around x, y, z when true
21  Bool_t RotX = kFALSE;
22  Bool_t RotY = kFALSE;
23  Bool_t RotZ = kFALSE;
24 
25  // main container to store all relevant alignment informatio
26  std::map<UInt_t, std::pair<TString, TGeoCombiTrans*>> misalignInfo;
27 
28  // Open file and get the stired TGeoManager
29  TString mcEngine = "TGeant3";
30  TString dir = getenv("VMCWORKDIR");
31 
32  TString outDir = "./data/";
33 
34  TString geoFile = "makeGeoFile_";
35  geoFile = outDir + geoFile + mcEngine + "_full.root";
36 
37  TFile* f = new TFile(geoFile, "READ");
38  if (!f) {
39  std::cout << "Could not open the file" << std::endl;
40  exit(1);
41  }
42 
43  TList* l = f->GetListOfKeys();
44 
45  TKey* key;
46  TIter next(l);
47 
48  TGeoManager* geoMan{nullptr};
49 
50  while ((key = (TKey*)next())) {
51  if (key->ReadObj()->InheritsFrom("TGeoManager")) {
52  geoMan = static_cast<TGeoManager*>(key->ReadObj());
53  break;
54  }
55  }
56 
57  if (!geoMan) {
58  f->Close();
59  std::cout << "Could not find TGeoManager inside the file" << std::endl;
60  exit(1);
61  }
62 
63  Float_t shiftx{0.};
64  Float_t shifty{0.};
65  Float_t shiftz{0.};
66  Float_t rotx{0.};
67  Float_t roty{0.};
68  Float_t rotz{0.};
69 
70  TGeoNode* topNode = gGeoManager->GetTopNode();
71  TObjArray* nodes = topNode->GetNodes();
72  for (Int_t iNode = 0; iNode < nodes->GetEntriesFast(); iNode++) {
73  TGeoNode* node = static_cast<TGeoNode*>(nodes->At(iNode));
74  if (!TString(node->GetName()).Contains("tutorial4"))
75  continue; // trd_vXXy top node, e.g. trd_v13a, trd_v14b
76  TGeoNode* station = node;
77  TObjArray* layers = station->GetNodes();
78  for (Int_t iLayer = 0; iLayer < layers->GetEntriesFast(); iLayer++) {
79  TGeoNode* layer = static_cast<TGeoNode*>(layers->At(iLayer));
80  TString path = TString("/") + topNode->GetName() + "/" + station->GetName() + "/" + layer->GetName();
81 
82  UInt_t hash = path.Hash();
83  if (ShiftX) {
84  shiftx = gRandom->Uniform(-1., 1.);
85  }
86  if (ShiftY) {
87  shifty = gRandom->Uniform(-1., 1.);
88  }
89  if (ShiftZ) {
90  shiftz = gRandom->Uniform(-1., 1.);
91  }
92 
93  TGeoTranslation* trans{nullptr};
94  TGeoRotation* rot{nullptr};
95  if (iLayer > 0 & iLayer < layers->GetEntriesFast() - 1) {
96  trans = new TGeoTranslation(shiftx, shifty, shiftz);
97  rot = new TGeoRotation("temp", rotx, roty, rotz);
98  } else {
99  trans = new TGeoTranslation(0., 0., 0.);
100  rot = new TGeoRotation("temp", 0., 0., 0.);
101  }
102  TGeoCombiTrans* combi = new TGeoCombiTrans(*trans, *rot);
103 
104  if (misalignInfo.find(hash) != misalignInfo.end()) {
105  cerr << "[ERROR ] Hash collision" << endl;
106  cerr << "[ERROR ] Hash value " << hash << " already exist." << endl;
107  auto it = misalignInfo.find(hash);
108  cerr << "[ERROR ] TGeoNode: " << it->second.first << " and " << path << " both have hash " << hash
109  << endl;
110  cerr << "[FATAL ] This should never happen. Stop execution." << endl;
111  exit(1);
112  }
113 
114  misalignInfo[hash] = make_pair(path, combi);
115  }
116  }
117 
118  // Close the file with the TGeoManager after all information was extracted
119  f->Close();
120 
121  cout << "There are in total " << misalignInfo.size() << " misalligned nodes" << endl;
122 
123  for (auto& elem : misalignInfo) {
124  cout << "Module " << elem.first << " " << elem.second.first << " " << elem.second.second << endl;
125  }
126 
127  // Write the paramter file
128  ofstream myfile;
129  myfile.open("example.par");
130 
131  myfile << "##############################################################################" << endl;
132  myfile << "# Class: FairTutorialDetMissallignPar" << endl;
133  myfile << "# Context: TestDefaultContext" << endl;
134  myfile << "##############################################################################" << endl;
135  myfile << "[FairTutorialDetMissallignPar]" << endl;
136  myfile << "//----------------------------------------------------------------------------" << endl;
137  myfile << "NrOfDetectors: Int_t " << misalignInfo.size() << endl;
138 
139  // vectors for the different components of the rotation matrix and
140  // the translation vector.
141  std::vector<Double_t> r11;
142  std::vector<Double_t> r12;
143  std::vector<Double_t> r13;
144  std::vector<Double_t> r21;
145  std::vector<Double_t> r22;
146  std::vector<Double_t> r23;
147  std::vector<Double_t> r31;
148  std::vector<Double_t> r32;
149  std::vector<Double_t> r33;
150  std::vector<Double_t> tx;
151  std::vector<Double_t> ty;
152  std::vector<Double_t> tz;
153  std::vector<TString> path;
154  std::vector<UInt_t> hash;
155 
156  Int_t counter = 0;
157 
158  myfile << "DetectorId: UInt_t \\ " << endl;
159  myfile << " ";
160  for (auto& elem : misalignInfo) {
161  const Double_t* rotation = elem.second.second->GetRotationMatrix();
162  const Double_t* translation = elem.second.second->GetTranslation();
163  r11.push_back(rotation[0]);
164  r12.push_back(rotation[1]);
165  r13.push_back(rotation[2]);
166  r21.push_back(rotation[3]);
167  r22.push_back(rotation[4]);
168  r23.push_back(rotation[5]);
169  r31.push_back(rotation[6]);
170  r32.push_back(rotation[7]);
171  r33.push_back(rotation[8]);
172  tx.push_back(translation[0]);
173  ty.push_back(translation[1]);
174  tz.push_back(translation[2]);
175  hash.push_back(elem.first);
176  path.push_back(elem.second.first);
177 
178  if ((8 == counter)) {
179  myfile << elem.first << " \\" << endl;
180  myfile << " ";
181  counter = 0;
182  } else {
183  myfile << elem.first << " ";
184  ++counter;
185  }
186  }
187  myfile << endl;
188 
189  // Due to performance issues when reading the parameters it is not possible
190  // to store the information about rotation and translation for each
191  // misaligned node. Instead store the values for one component of rotation
192  // or translation for all nodes in one array. This reduces the number of
193  // parameters from a (possible) very large number to 12 (9 components for
194  // the rotation matrix and 3 for the translation matrix. The proper
195  // matrices will be created when reading the parameter container.
196 
197  myfile << "R11: Double_t \\ " << endl;
198  myfile << " ";
199  writeVectorToFile(r11, &myfile);
200  myfile << "R12: Double_t \\ " << endl;
201  myfile << " ";
202  writeVectorToFile(r12, &myfile);
203  myfile << "R13: Double_t \\ " << endl;
204  myfile << " ";
205  writeVectorToFile(r13, &myfile);
206 
207  myfile << "R21: Double_t \\ " << endl;
208  myfile << " ";
209  writeVectorToFile(r21, &myfile);
210  myfile << "R22: Double_t \\ " << endl;
211  myfile << " ";
212  writeVectorToFile(r22, &myfile);
213  myfile << "23: Double_t \\ " << endl;
214  myfile << " ";
215  writeVectorToFile(r23, &myfile);
216 
217  myfile << "R31: Double_t \\ " << endl;
218  myfile << " ";
219  writeVectorToFile(r31, &myfile);
220  myfile << "R32: Double_t \\ " << endl;
221  myfile << " ";
222  writeVectorToFile(r32, &myfile);
223  myfile << "33: Double_t \\ " << endl;
224  myfile << " ";
225  writeVectorToFile(r33, &myfile);
226 
227  myfile << "TX: Double_t \\ " << endl;
228  myfile << " ";
229  writeVectorToFile(tx, &myfile);
230  myfile << "TY: Double_t \\ " << endl;
231  myfile << " ";
232  writeVectorToFile(ty, &myfile);
233  myfile << "TZ: Double_t \\ " << endl;
234  myfile << " ";
235  writeVectorToFile(tz, &myfile);
236 
237  myfile << "##############################################################################" << endl;
238 
239  myfile.close();
240 
241  /*
242  // To check if the vector sorting works and one gets in the end the same
243  // results as when using the map
244  myfile.open ("test1.txt");
245  for (Int_t iCount = 0; iCount < r11.size(); iCount++) {
246  myfile << "Module "<< hash[iCount] << " " << path[iCount] << " "
247  << r11[iCount] << " "
248  << r12[iCount] << " "
249  << r13[iCount] << " "
250  << r21[iCount] << " "
251  << r22[iCount] << " "
252  << r23[iCount] << " "
253  << r31[iCount] << " "
254  << r32[iCount] << " "
255  << r33[iCount] << " "
256  << tx[iCount] << " "
257  << ty[iCount] << " "
258  << tz[iCount]
259  << endl;
260  }
261  myfile.close();
262 */
263 
264  myfile.open("misalignment.txt");
265  for (auto& elem : misalignInfo) {
266  const Double_t* rotation = elem.second.second->GetRotationMatrix();
267  const Double_t* translation = elem.second.second->GetTranslation();
268  myfile << "Module " << elem.first << " " << elem.second.first << " " << rotation[0] << " " << rotation[1] << " "
269  << rotation[2] << " " << rotation[3] << " " << rotation[4] << " " << rotation[5] << " " << rotation[6]
270  << " " << rotation[7] << " " << rotation[8] << " " << translation[0] << " " << translation[1] << " "
271  << translation[2] << endl;
272  }
273  myfile.close();
274 }
275 
276 void writeVectorToFile(std::vector<Double_t>& vec, ofstream* myfile)
277 {
278  Int_t counter = 0;
279  for (auto& elem : vec) {
280  if ((8 == counter)) {
281  *myfile << elem << " \\" << endl;
282  *myfile << " ";
283  counter = 0;
284  } else {
285  *myfile << elem << " ";
286  ++counter;
287  }
288  }
289  *myfile << endl;
290 }
run_tutorial4_createParameterFile
void run_tutorial4_createParameterFile()
Definition: run_tutorial4_createParameterFile.C:12
writeVectorToFile
void writeVectorToFile(std::vector< Double_t > &vec, ofstream *myfile)
Definition: run_tutorial4_createParameterFile.C:276