FairRoot/PandaRoot
testAlignMatrices.C
Go to the documentation of this file.
1 #include <TH1D.h>
2 #include <TCanvas.h>
3 #include <TGeoManager.h>
4 #include <TGeoMatrix.h>
5 #include <TGeoPhysicalNode.h>
6 #include "PndLmdGeometryHelper.h"
7 #include <TGeoNode.h>
8 #include <TFile.h>
9 #include <TChain.h>
10 #include <TClonesArray.h>
11 
12 #include <FairRunSim.h>
13 
14 #include <iostream>
15 #include <fstream>
16 #include <map>
17 #include <string>
18 #include <sstream>
19 #include <PndLmdHitPair.h>
20 
21 #include "json.hpp"
22 
23 std::map<std::string, TGeoHMatrix> *matricesMisaligned;
24 std::map<std::string, TGeoHMatrix> *matricesIdeal;
25 PndLmdGeometryHelper *helper;
26 
27 using std::cout;
28 using std::string;
29 using json = nlohmann::json;
30 
31 void initDummySimulation() {
32  cout << "creating dummy simulation to populate gGeoManager...\n";
33  TString geometryFile = "Luminosity-Detector.root";
34  TString simOutput = "./dummy.root";
35  FairRunSim *fRun = new FairRunSim();
36  fRun->SetName("TGeant4");
37  fRun->SetOutputFile(simOutput);
38  fRun->SetMaterials("media_pnd.geo");
39  FairModule *Cave = new PndCave("CAVE");
40  Cave->SetGeometryFileName("pndcave.geo");
41  fRun->AddModule(Cave);
42  FairModule *Pipe = new PndPipe("PIPE");
43  Pipe->SetGeometryFileName("beampipe_201309.root");
44  fRun->AddModule(Pipe);
45  FairModule *Magnet = new PndMagnet("MAGNET");
46  Magnet->SetGeometryFileName("FullSolenoid_V842.root");
47  fRun->AddModule(Magnet);
48  FairModule *Dipole = new PndMagnet("MAGNET");
49  Dipole->SetGeometryFileName("dipole.geo");
50  fRun->AddModule(Dipole);
51  PndLmdDetector *Lum = new PndLmdDetector("LUM", kTRUE);
52  Lum->SetExclusiveSensorType("LumActive"); //ignore MVD
53  Lum->SetGeometryFileName(geometryFile);
54  fRun->AddModule(Lum);
55  fRun->Init();
56 }
57 
58 void applyMisalignmentToGGeoManager(TString filename) {
59  bool misaligned = true;
60  if (misaligned) {
61  //load matrices
62  TFile *misalignmentMatrixRootfile = new TFile(filename, "READ");
63 
64  if (misalignmentMatrixRootfile->IsOpen()) {
65  std::map<std::string, TGeoHMatrix> *matrices;
66 
67  gDirectory->GetObject("PndLmdMisalignMatrices", matrices);
68  misalignmentMatrixRootfile->Close();
69 
70  std::cout << matrices->size()
71  << " matrices successfully read from file.\napplying misalignment.\n";
72 
73  //iterate over matrices
74  for (auto const& entry : *matrices) {
75  TString volPath = entry.first;
76 
77  gGeoManager->cd(volPath);
78 
79  TGeoNode* currentNode = gGeoManager->GetCurrentNode();
80  TGeoMatrix* matrixToNode = currentNode->GetMatrix();
81 
82  TGeoHMatrix misalignedMatrixToNode = *matrixToNode * entry.second;
83 
84  //this is just for clarity, can probably be removed
85  TGeoHMatrix* newMatrixToNode = new TGeoHMatrix(misalignedMatrixToNode); // new matrix, representing real position
86 
87  TGeoPhysicalNode* physicalNode = gGeoManager->MakePhysicalNode(volPath);
88 
89  physicalNode->Align(newMatrixToNode);
90  }
91  std::cout << "all misalignments applied.\n";
92  }
93  else {
94  std::cout << "file could not be read\n";
95  }
96  }
97 }
98 
99 void unApplyMisalignmentToGGeoManager(TString filename) {
100  bool misaligned = true;
101  if (misaligned) {
102  //load matrices
103  TFile *misalignmentMatrixRootfile = new TFile(filename, "READ");
104 
105  if (misalignmentMatrixRootfile->IsOpen()) {
106  std::map<std::string, TGeoHMatrix> *matrices;
107 
108  gDirectory->GetObject("PndLmdMisalignMatrices", matrices);
109  misalignmentMatrixRootfile->Close();
110 
111  std::cout << "resetting geometry!.\n";
112 
113  //iterate over matrices
114  for (auto const& entry : *matrices) {
115  TString volPath = entry.first;
116 
117  gGeoManager->cd(volPath);
118 
119  TGeoNode* currentNode = gGeoManager->GetCurrentNode();
120  TGeoMatrix* matrixToNode = currentNode->GetMatrix();
121 
122  // invert again!
123  TGeoHMatrix inverseMisalign = entry.second.Inverse();
124 
125  TGeoHMatrix misalignedMatrixToNode = *matrixToNode * inverseMisalign;
126 
127  //this is just for clarity, can probably be removed
128  TGeoHMatrix* newMatrixToNode = new TGeoHMatrix(misalignedMatrixToNode); // new matrix, representing real position
129 
130  TGeoPhysicalNode* physicalNode = gGeoManager->MakePhysicalNode(volPath);
131 
132  physicalNode->Align(newMatrixToNode);
133  }
134  }
135  else {
136  std::cout << "file could not be read\n";
137  }
138  }
139 }
140 
141 //this function is ugly and stiched together into a barely functioning blob of ugly.
142 //I know. I'm going to fix this (probably, some time in the future)
143 TGeoHMatrix readMatrixFromDisk(std::string filename) {
144 
145  int lines = 0;
146  stringstream *valueStream = new stringstream();
147 
148  ifstream ifs;
149  ifs.open(filename.c_str());
150  if (!ifs.is_open()) {
151  cout << "could not read ";
152  cout << filename;
153  cout << "! aborting." << "\n";
154  return NULL;
155  }
156  string linebuffer; // = new string();
157  while (std::getline(ifs, linebuffer)) {
158  *valueStream << linebuffer << "\n";
159  lines++;
160  }
161  ifs.close();
162 
163  std::stringstream *ss;
164  std::stringstream iss;
165  int i = 0;
166 
167  ss = valueStream;
168  std::string line, token;
169  std::vector<std::vector<double> > data;
170  while (getline(*ss, line)) {
171  iss << line;
172  std::vector<double> tempvec;
173  while (getline(iss, token, ',')) {
174  tempvec.push_back(std::stod(token));
175  }
176  data.push_back(tempvec);
177  ++i;
178  iss.clear();
179  }
180  ss->str("");
181  delete ss;
182  iss.str("");
183 
184  //cast to TGeoHMatrix
185 
186  if (data.size() < 1) {
187  cout << "warning! can't read matrix from file " << filename << "\n";
188  exit(1);
189  }
190  if (data[0].size() < 1) {
191  cout << "warning! can't read matrix from file " << filename << "\n";
192  exit(1);
193  }
194 
195  double rot[9];
196  double trans[3];
197 
198  int rows, columns;
199  rows = data.size();
200  columns = data[0].size();
201 
202  if (rows != 4 || columns != 4) {
203  cerr << "wrong matrix dimensions.\n";
204  exit(1);
205  }
206 
207  rot[0] = data[0][0];
208  rot[1] = data[0][1];
209  rot[2] = data[0][2];
210  rot[3] = data[1][0];
211  rot[4] = data[1][1];
212  rot[5] = data[1][2];
213  rot[6] = data[2][0];
214  rot[7] = data[2][1];
215  rot[8] = data[2][2];
216 
217  trans[0] = data[0][3];
218  trans[1] = data[1][3];
219  trans[2] = data[2][3];
220 
221  TGeoHMatrix result;
222 
223  result.SetRotation(rot);
224  result.SetDx(trans[0]);
225  result.SetDy(trans[1]);
226  result.SetDz(trans[2]);
227 
228  return result;
229 }
230 
231 std::map<std::string, TGeoHMatrix> * readRootMatrices(TString filename) {
232  cout << "reading files from file: " << filename << "\n";
233 
234  TFile *misalignmentMatrixRootfile = new TFile(filename, "READ");
235 
236  if (misalignmentMatrixRootfile->IsOpen()) {
237  std::map<std::string, TGeoHMatrix> *matrices;
238 
239  gDirectory->GetObject("PndLmdMisalignMatrices", matrices);
240  misalignmentMatrixRootfile->Close();
241  cout << matrices->size() << " matrices successfully read from file.\n";
242 
243  return matrices;
244  }
245  else {
246  cout << "file could not be read\n";
247  return NULL;
248  }
249 
250 }
251 
252 void printMatrixDiff(TGeoHMatrix mat1, TGeoHMatrix mat2) {
253 
254  const double *translationDesign = mat1.GetTranslation();
255  const double *translationICP = mat2.GetTranslation();
256  const double *rotationDesign = mat1.GetRotationMatrix();
257  const double *rotationICP = mat2.GetRotationMatrix();
258 
259  double dx = translationDesign[0] - translationICP[0];
260  double dy = translationDesign[1] - translationICP[1];
261  double da = rotationDesign[1] - rotationICP[1];
262 
263  cout << "\n=-=-=-=-=-=-=-=-=-=\n";
264  cout << "differences are:\ndx: " << dx * 1e4 << "\ndy: " << dy * 1e4 << "\n";
265  cout << "=-=-=-=-=-=-=-=-=-=\n";
266 }
267 
268 // this depends on the geometry currently loaded to gGeoManager!
269 TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo) {
270 
271  // the passive transformation that base-transforms a vector in Panda Global to the system of Sensor x
272  TGeoHMatrix toSen0 = helper->getMatrixPndGlobalToSensor(sensorOne);
273  TGeoHMatrix toSen1 = helper->getMatrixPndGlobalToSensor(sensorTwo);
274 
275  return toSen0.Inverse() * toSen1;
276 
277 }
278 
279 TGeoHMatrix baseTransformation(TGeoHMatrix &input, TGeoHMatrix &toBaseMatrix) {
280  return TGeoHMatrix(toBaseMatrix * input * toBaseMatrix.Inverse());
281 }
282 
283 TGeoHMatrix getOverlapMatrixLikeICP(PndLmdOverlapInfo info) {
284 
285  //prepare misalignment matrices
286  TGeoHMatrix misalignmentToSensor0 = (*matricesMisaligned)[info.path1];
287  TGeoHMatrix misalignmentToSensor1 = (*matricesMisaligned)[info.path2];
288  // transform second misaignment Matrix to new base
289  TGeoHMatrix Sen0ToSen1 = getMatrixSensorToSensor(info.id1, info.id2);
290  TGeoHMatrix misSen1inSen0 = baseTransformation(misalignmentToSensor1, Sen0ToSen1);
291 
292  return misalignmentToSensor0.Inverse() * misSen1inSen0;
293 }
294 
295 // WARNING. the geometry inside gGeomanager must be aligned for this to work!
296 TGeoHMatrix getMisalignedOverlapFromGeoManager(PndLmdOverlapInfo &info) {
297 
298  TGeoHMatrix sen0to1 = getMatrixSensorToSensor(info.id1, info.id2);
299  TGeoHMatrix sen0to1ICPcorr = getOverlapMatrixLikeICP(info);
300  // this order is important!
301  return sen0to1ICPcorr * sen0to1;
302 }
303 
304 // WARNING. the geometry inside gGeomanager must be aligned for this to work!
305 TGeoHMatrix getMisalignedOverlapFromICP(PndLmdOverlapInfo &info, std::string ICPmatrix) {
306 
307  TGeoHMatrix sen0to1 = getMatrixSensorToSensor(info.id1, info.id2);
308  TGeoHMatrix sen0to1ICPcorr = readMatrixFromDisk(ICPmatrix);
309  // this order is important!
310  return sen0to1ICPcorr * sen0to1;
311 }
312 
313 PndLmdOverlapInfo& getSmallOverlapInfo(std::vector<PndLmdOverlapInfo> &infos, int smallOverlap) {
314  //PndLmdOverlapInfo result;
315  for (auto &info : infos) {
316  int smallOverlapHere = info.overlapID % 10;
317  if (smallOverlap == smallOverlapHere) {
318  return info;
319  }
320  }
321  return infos[0];
322 }
323 
324 std::vector<double> getMatrixDiff(TGeoHMatrix &mat1, TGeoHMatrix &mat2) {
325  const double *translationDesign = mat1.GetTranslation();
326  const double *translationICP = mat2.GetTranslation();
327  const double *rotationDesign = mat1.GetRotationMatrix();
328  const double *rotationICP = mat2.GetRotationMatrix();
329 
330  double dx = translationDesign[0] - translationICP[0];
331  double dy = translationDesign[1] - translationICP[1];
332  double da = rotationDesign[1] - rotationICP[1];
333 
334  std::vector<double> result;
335  result.push_back(dx);
336  result.push_back(dy);
337  result.push_back(da);
338  return result;
339 }
340 
341 void buildCyclic(int alignParam, string cut) {
342 
343  TString idealMatrices = "idealMatrices.root";
344 
345  TString misalignedMatrices = "";
346  std::string pathPrefix = "/home/arbeit/RedPro3TB/simulationData/";
347  //std::string pathPost = "LMDmatrices-cut" + cut + "/";
348  std::string pathPost = "LMDmatrices-python-cut-" + cut + "-2D/";
349  std::string path;
350  std::string pdfPath;
351 
352  switch (alignParam) {
353  case 0:
354  path = "2018-08-himster2-aligned/";
355  pdfPath = "aligned/cut-" + cut + "/";
356  break;
357  case 10:
358  misalignedMatrices = "misalignMatrices-SensorsOnly-10.root";
359  path = "2018-08-himster2-misalign-10u/";
360  pdfPath = "misalign-10u/cut-" + cut + "/";
361  break;
362  case 50:
363  misalignedMatrices = "misalignMatrices-SensorsOnly-50.root";
364  path = "2018-08-himster2-misalign-50u/";
365  pdfPath = "misalign-50u/cut-" + cut + "/";
366  break;
367  case 100:
368  misalignedMatrices = "misalignMatrices-SensorsOnly-100.root";
369  path = "2018-08-himster2-misalign-100u/";
370  pdfPath = "misalign-100u/cut-" + cut + "/";
371  break;
372  case 150:
373  misalignedMatrices = "misalignMatrices-SensorsOnly-150.root";
374  path = "2018-08-himster2-misalign-150u/";
375  pdfPath = "misalign-150u/cut-" + cut + "/";
376  break;
377  case 200:
378  misalignedMatrices = "misalignMatrices-SensorsOnly-200.root";
379  path = "2018-08-himster2-misalign-200u/";
380  pdfPath = "misalign-200u/cut-" + cut + "/";
381  break;
382  case 250:
383  misalignedMatrices = "misalignMatrices-SensorsOnly-250.root";
384  path = "2018-08-himster2-misalign-250u/";
385  pdfPath = "misalign-250u/cut-" + cut + "/";
386  break;
387  }
388 
389  matricesMisaligned = readRootMatrices(misalignedMatrices);
390  matricesIdeal = readRootMatrices(idealMatrices);
391 
392  std::map<int, TGeoHMatrix> matricesByHand;
393  std::map<int, TGeoHMatrix> matricesByGeoManager;
394 
395  bool compareWithIdeal = false;
396  string ext = "cm.mat";
397 
398  if (compareWithIdeal) {
399  cout << "getting matrices by hand | ICP-like matrices from gGeoManager\n";
400  }
401  else {
402  cout << "getting matrices by hand | real ICP matrices\n";
403  }
404 
405  path = pathPrefix + path + pathPost;
406 
407  gSystem->Exec(("mkdir -p " + pdfPath).c_str());
408 
409  for (int iHalf = 0; iHalf < 2; iHalf++) {
410  for (int iPlane = 0; iPlane < 4; iPlane++) {
411  for (int iModule = 0; iModule < 5; iModule++) {
412 
413  // the order in the vector is not defined, use getSmallOverlapInfo() for that
414  std::vector<PndLmdOverlapInfo> overlaps = helper->getOverlapInfos(iHalf, iPlane, iModule);
415 
416  // load all 9 overlap matrices
417  TGeoHMatrix mat0;
418  TGeoHMatrix mat1;
419  TGeoHMatrix mat2;
420  TGeoHMatrix mat3;
421  TGeoHMatrix mat4;
422  TGeoHMatrix mat5;
423  TGeoHMatrix mat6;
424  TGeoHMatrix mat7;
425  TGeoHMatrix mat8;
426 
427  if (compareWithIdeal) {
428  mat0 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 0));
429  mat1 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 1));
430  mat2 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 2));
431  mat3 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 3));
432  mat4 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 4));
433  mat5 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 5));
434  mat6 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 6));
435  mat7 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 7));
436  mat8 = getMisalignedOverlapFromGeoManager(getSmallOverlapInfo(overlaps, 8));
437  }
438  else {
439  string icp0 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 0).overlapID) + ext;
440  mat0 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 0), icp0);
441  string icp1 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 1).overlapID) + ext;
442  mat1 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 1), icp1);
443  string icp2 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 2).overlapID) + ext;
444  mat2 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 2), icp2);
445  string icp3 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 3).overlapID) + ext;
446  mat3 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 3), icp3);
447  string icp4 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 4).overlapID) + ext;
448  mat4 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 4), icp4);
449  string icp5 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 5).overlapID) + ext;
450  mat5 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 5), icp5);
451  string icp6 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 6).overlapID) + ext;
452  mat6 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 6), icp6);
453  string icp7 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 7).overlapID) + ext;
454  mat7 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 7), icp7);
455  string icp8 = path + "/m" + std::to_string(getSmallOverlapInfo(overlaps, 8).overlapID) + ext;
456  mat8 = getMisalignedOverlapFromICP(getSmallOverlapInfo(overlaps, 8), icp8);
457  }
458 
459  // get 0->1 and 5->6
460  TGeoHMatrix sen0to1 = getMatrixSensorToSensor(0, 1);
461  TGeoHMatrix sen5to6 = getMatrixSensorToSensor(5, 6);
462 
463  TGeoHMatrix sen1to2 = mat4 * mat5.Inverse();
464  TGeoHMatrix sen1to3 = mat4 * mat1.Inverse();
465  TGeoHMatrix sen1to4a = mat4 * mat5.Inverse() * mat6 * mat2.Inverse();
466  TGeoHMatrix sen1to4b = mat4 * mat1.Inverse() * mat7 * mat8.Inverse();
467 
468  TGeoHMatrix sen1to5 = mat0; // well, kind of...
469  TGeoHMatrix sen1to6 = mat4 * mat1.Inverse() * mat3;
470  TGeoHMatrix sen1to7 = mat4 * mat1.Inverse() * mat7;
471  TGeoHMatrix sen1to8 = mat4;
472  TGeoHMatrix sen1to9a = mat4 * mat5.Inverse() * mat6;
473  TGeoHMatrix sen1to9b = mat4 * mat1.Inverse() * mat7 * mat8.Inverse() * mat2;
474 
475  int matrixId = 1000 * iHalf + 100 * iPlane + 10 * iModule;
476 
477  matricesByHand[matrixId + 2] = sen1to2;
478  matricesByHand[matrixId + 3] = sen1to3;
479  matricesByHand[matrixId + 4] = sen1to4a;
480  matricesByHand[matrixId + 5] = sen1to5; // this is not yet correct
481  matricesByHand[matrixId + 6] = sen1to6;
482  matricesByHand[matrixId + 7] = sen1to7;
483  matricesByHand[matrixId + 8] = sen1to8;
484  matricesByHand[matrixId + 9] = sen1to9a;
485  }
486  }
487  }
488 
489  if (alignParam > 0) {
490  applyMisalignmentToGGeoManager(misalignedMatrices);
491  }
492 
493  cout << "getting matrices by geoManager after misalignment\n";
494  for (int iHalf = 0; iHalf < 2; iHalf++) {
495  for (int iPlane = 0; iPlane < 4; iPlane++) {
496  for (int iModule = 0; iModule < 5; iModule++) {
497 
498  // the order in the vector is not defined, use getSmallOverlapInfo() for that
499  std::vector<PndLmdOverlapInfo> overlaps = helper->getOverlapInfos(iHalf, iPlane, iModule);
500  // get everything from gGeoManager
501  int matrixId = 1000 * iHalf + 100 * iPlane + 10 * iModule;
502 
503  PndLmdOverlapInfo info0 = getSmallOverlapInfo(overlaps, 0);
504  PndLmdOverlapInfo info1 = getSmallOverlapInfo(overlaps, 1);
505  PndLmdOverlapInfo info2 = getSmallOverlapInfo(overlaps, 2);
506  PndLmdOverlapInfo info3 = getSmallOverlapInfo(overlaps, 3);
507  PndLmdOverlapInfo info4 = getSmallOverlapInfo(overlaps, 4);
508  PndLmdOverlapInfo info5 = getSmallOverlapInfo(overlaps, 5);
509  PndLmdOverlapInfo info6 = getSmallOverlapInfo(overlaps, 6);
510  PndLmdOverlapInfo info7 = getSmallOverlapInfo(overlaps, 7);
511  PndLmdOverlapInfo info8 = getSmallOverlapInfo(overlaps, 8);
512 
513  // get sensors 1 and 2 from the overlapInfos
514  int sensor0 = info0.id1;
515  int sensor1 = info4.id1;
516  int sensor2 = info5.id1;
517  int sensor3 = info1.id1;
518  int sensor4 = info2.id1;
519  int sensor5 = info0.id2;
520  int sensor6 = info3.id2;
521  int sensor7 = info7.id2;
522  int sensor8 = info4.id2;
523  int sensor9 = info6.id2;
524 
525  TGeoHMatrix sen1to2 = getMatrixSensorToSensor(sensor1, sensor2);
526  TGeoHMatrix sen1to3 = getMatrixSensorToSensor(sensor1, sensor3);
527  TGeoHMatrix sen1to4a = getMatrixSensorToSensor(sensor1, sensor4);
528  TGeoHMatrix sen1to5 = getMatrixSensorToSensor(sensor0, sensor5);
529  TGeoHMatrix sen1to6 = getMatrixSensorToSensor(sensor1, sensor6);
530  TGeoHMatrix sen1to7 = getMatrixSensorToSensor(sensor1, sensor7);
531  TGeoHMatrix sen1to8 = getMatrixSensorToSensor(sensor1, sensor8);
532  TGeoHMatrix sen1to9a = getMatrixSensorToSensor(sensor1, sensor9);
533 
534  matricesByGeoManager[matrixId + 2] = sen1to2;
535  matricesByGeoManager[matrixId + 3] = sen1to3;
536  matricesByGeoManager[matrixId + 4] = sen1to4a;
537  matricesByGeoManager[matrixId + 5] = sen1to5; // this is not yet correct
538  matricesByGeoManager[matrixId + 6] = sen1to6;
539  matricesByGeoManager[matrixId + 7] = sen1to7;
540  matricesByGeoManager[matrixId + 8] = sen1to8;
541  matricesByGeoManager[matrixId + 9] = sen1to9a;
542  }
543  }
544  }
545 
546  //reset geometry for next run!
547  if (alignParam > 0) {
548  unApplyMisalignmentToGGeoManager(misalignedMatrices);
549  }
550 
551  cout << "making histograms\n";
552  double dalpha, dx, dy, dxTarget, dyTarget, daTarget, dxICP, dyICP, daICP;
553  TH1D histA("h1", "h1", 25, -1, -1);
554  TH1D histX("h2", "h2", 25, -1, -1);
555  TH1D histY("h3", "h3", 25, -1, -1);
556 
557  for (auto &i : matricesByHand) {
558  int key = i.first;
559  auto thisDiff = getMatrixDiff(matricesByHand[key], matricesByGeoManager[key]);
560 
561  // hist das shizzle
562  histA.Fill(thisDiff[2] * 1e3);
563  histX.Fill(thisDiff[0] * 1e4);
564  histY.Fill(thisDiff[1] * 1e4);
565  }
566 
567  //name histograms, clean up later
568  histA.GetXaxis()->SetTitle("d#alpha [#mrad]");
569  histX.GetXaxis()->SetTitle("dx [#mum]");
570  histY.GetXaxis()->SetTitle("dy [#mum]");
571  histA.GetYaxis()->SetTitle("Entries");
572  histX.GetYaxis()->SetTitle("Entries");
573  histY.GetYaxis()->SetTitle("Entries");
574  histA.SetTitle("d#alpha-combined");
575  histX.SetTitle("dx-combined");
576  histY.SetTitle("dy-combined");
577 
578  TCanvas canvas("c1", "c1", 800, 600);
579  canvas.cd();
580  histA.Draw();
581  canvas.Print((pdfPath + "dAlpha-combined.pdf").c_str());
582 
583  histX.Draw();
584  canvas.Print((pdfPath + "dx-combined.pdf").c_str());
585 
586  histY.Draw();
587  canvas.Print((pdfPath + "dy-combined.pdf").c_str());
588 
589 }
590 
591 void saveMatricesToJson() {
592  auto overlaps = helper->getOverlapInfos();
593 
594  cout << "got " << overlaps.size() << " overlaps\n";
595 
596  json j;
597 
598  for (auto &overlap : overlaps) {
599  //cout << overlap << "\n";
600 
601  double *mat1 = new double[16];
602  overlap.mat1.GetHomogenousMatrix(mat1);
603  double *mat2 = new double[16];
604  overlap.mat2.GetHomogenousMatrix(mat2);
605 
606  j[std::to_string(overlap.overlapID)]["overlapID"] = overlap.overlapID;
607  j[std::to_string(overlap.overlapID)]["id1"] = overlap.id1;
608  j[std::to_string(overlap.overlapID)]["id2"] = overlap.id2;
609  j[std::to_string(overlap.overlapID)]["path1"] = overlap.path1;
610  j[std::to_string(overlap.overlapID)]["path2"] = overlap.path2;
611  j[std::to_string(overlap.overlapID)]["matrix1"] = {
612  mat1[0], mat1[1], mat1[2], mat1[12],
613  mat1[4], mat1[5], mat1[6], mat1[13],
614  mat1[8], mat1[9], mat1[10], mat1[14],
615  mat1[3], mat1[7], mat1[11], mat1[15]
616  };
617  j[std::to_string(overlap.overlapID)]["matrix2"] = {
618  mat2[0], mat2[1], mat2[2], mat2[12],
619  mat2[4], mat2[5], mat2[6], mat2[13],
620  mat2[8], mat2[9], mat2[10], mat2[14],
621  mat2[3], mat2[7], mat2[11], mat2[15]
622  };
623  }
624 
625  cout << "got " << j.size() << " json objects\n";
626 
627  std::ofstream o("detectorMatricesIdeal.json");
628  o << std::setw(2) << j << std::endl;
629 
630 }
631 
632 void saveDetMatricesToJson() {
633 
634  // what do I want to save?
635 
636  /*
637  panda -> lmd
638  lmd -> top half
639  lmd -> bottom half
640  halfes -> planes
641 
642  */
643  cout << setprecision(16);
644  gGeoManager->cd("/cave_1/lmd_root_0/");
645  auto mat = gGeoManager->GetCurrentMatrix();
646  mat->Print();
647 
648  return;
649 }
650 
651 // compares ICP matrices found with geometry shiftes vs
652 // ICP matrices found with data shifted
653 void compareShiftDataShiftGeo() {
654 
655  std::string shiftDataPath = "/home/arbeit/RedPro3TB/simulationData/2018-10-himster2-misalignData-200u/";
656  std::string shiftGeoPath = "/home/arbeit/RedPro3TB/simulationData/2018-08-himster2-misalign-200u/";
657  std::string pathPre = "LMDmatrices-python-cut-";
658  std::string pathPost = "-2D";
659  std::string ext = "cm.mat";
660  std::string cutPath = "";
661  std::string filename;
662  std::string pdfPathPre = "ShiftGeovsShiftData/";
663  std::string pdfPath;
664 
665  // for now, this is only for 100u as the Himster2 is down
666 
667  // do for all cuts
668  //std::vector<string> cuts { "0", "0.5", "1", "3", "5" };
669  std::vector<string> cuts { "0"};
670 
671  for (auto &cut : cuts) {
672 
673  TH1D histA("h1", "h1", 25, -1, -1);
674  TH1D histX("h2", "h2", 25, -1, -1);
675  TH1D histY("h3", "h3", 25, -1, -1);
676 
677  // prepare cut path and pdf path, make pdf path
678  cutPath = pathPre + cut + pathPost;
679  pdfPath = pdfPathPre + "cut-" + cut + "/";
680  gSystem->Exec(("mkdir -p " + pdfPath).c_str());
681 
682  // for all overlaps
683  auto overlaps = helper->getOverlapInfos();
684  for (auto &overlap : overlaps) {
685 
686  // read shifted Geo
687  filename = shiftGeoPath + cutPath + "/m" + std::to_string(overlap.overlapID) + ext;
688  auto matShiftGeo = readMatrixFromDisk(filename);
689 
690  // read shifted data
691  filename = shiftDataPath + cutPath + "/m" + std::to_string(overlap.overlapID) + ext;
692  auto matShiftData = readMatrixFromDisk(filename);
693 
694  //matShiftData = TGeoHMatrix(); // just testing if we're really comparing the correct matrices
695 
696  // compare
697  auto thisDiff = getMatrixDiff(matShiftGeo, matShiftData);
698 
699  // hist das shizzle
700  histA.Fill(thisDiff[2] * 1e3);
701  histX.Fill(thisDiff[0] * 1e4);
702  histY.Fill(thisDiff[1] * 1e4);
703 
704  }
705 
706  //name histograms, clean up later
707  histA.GetXaxis()->SetTitle("d#alpha [#mrad]");
708  histX.GetXaxis()->SetTitle("dx [#mum]");
709  histY.GetXaxis()->SetTitle("dy [#mum]");
710  histA.GetYaxis()->SetTitle("Entries");
711  histX.GetYaxis()->SetTitle("Entries");
712  histY.GetYaxis()->SetTitle("Entries");
713  histA.SetTitle("ICP | shift Geo vs Shift Data, d#alpha");
714  histX.SetTitle("ICP | shift Geo vs Shift Data, px");
715  histY.SetTitle("ICP | shift Geo vs Shift Data, py");
716 
717  TCanvas canvas("c1", "c1", 800, 600);
718  canvas.cd();
719  histA.Draw();
720  canvas.Print((pdfPath + "dAlpha.pdf").c_str());
721 
722  histX.Draw();
723  canvas.Print((pdfPath + "dx.pdf").c_str());
724 
725  histY.Draw();
726  canvas.Print((pdfPath + "dy.pdf").c_str());
727 
728  }
729 }
730 
731 void compareICPmatrices(int alignParam, string cut) {
732 
733  TString idealMatrices = "idealMatrices.root";
734 
735  TString misalignedMatrices = "";
736  std::string pathMisalign;
737  std::string pathPrefix = "/home/arbeit/RedPro3TB/simulationData/";
738  //std::string pathPost = "LMDmatrices-cut" + cut + "/";
739  std::string pathPost = "LMDmatrices-python-cut-" + cut + "-2D/";
740  std::string pdfPath;
741  std::string ext = "cm.mat";
742 
743  std::string fullPath = "/home/roman/temp/himster2misaligned-move-data/Pairs/LMDmatrices/";
744 
745  switch (alignParam) {
746  case 0:
747  pathMisalign = "2018-08-himster2-aligned/";
748  pdfPath = "aligned/cut-" + cut + "/singleMatrices/";
749  break;
750  case 10:
751  misalignedMatrices = "misalignMatrices-SensorsOnly-10.root";
752  pathMisalign = "2018-08-himster2-misalign-10u/";
753  pdfPath = "misalign-10u/cut-" + cut + "/singleMatrices/";
754  break;
755  case 50:
756  misalignedMatrices = "misalignMatrices-SensorsOnly-50.root";
757  pathMisalign = "2018-08-himster2-misalign-50u/";
758  pdfPath = "misalign-50u/cut-" + cut + "/singleMatrices/";
759  break;
760  case 100:
761  misalignedMatrices = "misalignMatrices-SensorsOnly-100.root";
762  pathMisalign = "2018-08-himster2-misalign-100u/";
763  pdfPath = "misalign-100u/cut-" + cut + "/singleMatrices/";
764  break;
765  case 150:
766  misalignedMatrices = "misalignMatrices-SensorsOnly-150.root";
767  pathMisalign = "2018-08-himster2-misalign-150u/";
768  pdfPath = "misalign-150u/cut-" + cut + "/singleMatrices/";
769  break;
770  case 200:
771  misalignedMatrices = "misalignMatrices-SensorsOnly-200.root";
772  pathMisalign = "2018-08-himster2-misalign-200u/";
773  pdfPath = "misalign-200u/cut-" + cut + "/singleMatrices/";
774  break;
775  case 250:
776  misalignedMatrices = "misalignMatrices-SensorsOnly-250.root";
777  pathMisalign = "2018-08-himster2-misalign-250u/";
778  pdfPath = "misalign-250u/cut-" + cut + "/singleMatrices/";
779  break;
780  }
781 
782  // DON'T DO THIS, the function getOverlapLikeICP already considers misalignment
783  if (alignParam != 0) {
784  matricesMisaligned = readRootMatrices(misalignedMatrices);
785  }
786 
787  matricesIdeal = readRootMatrices(idealMatrices);
788 
789  std::string path = pathPrefix + pathMisalign + pathPost;
790  gSystem->Exec(("mkdir -p " + pdfPath).c_str());
791 
792  cout << "making histograms\n";
793  double dalpha, dx, dy, dxTarget, dyTarget, daTarget, dxICP, dyICP, daICP;
794  TH1D histA("h1", "h1", 25, -1, -1);
795  TH1D histX("h2", "h2", 25, -1, -1);
796  TH1D histY("h3", "h3", 25, -1, -1);
797 
798  auto overlaps = helper->getOverlapInfos();
799  for (auto &overlap : overlaps) {
800 
801  std::string filename;
802 
803  filename = path + "/m" + std::to_string(overlap.overlapID) + ext;
804 
805  TGeoHMatrix matIdeal;
806 
807  if (alignParam != 0) {
808  matIdeal = getOverlapMatrixLikeICP(overlap);
809  }
810  else {
811  //should be identity matrix already
812  }
813 
814  auto matICP = readMatrixFromDisk(filename);
815  auto thisDiff = getMatrixDiff(matIdeal, matICP);
816 
817  if (thisDiff[1] * 1e4 > 40.0) {
818  cout << "Alarm! " << overlap.overlapID << "\n";
819  }
820 
821  // hist das shizzle
822  histA.Fill(thisDiff[2] * 1e3);
823  histX.Fill(thisDiff[0] * 1e4);
824  histY.Fill(thisDiff[1] * 1e4);
825  }
826 
827  //name histograms, clean up later
828  histA.GetXaxis()->SetTitle("d#alpha [#mrad]");
829  histX.GetXaxis()->SetTitle("dx [#mum]");
830  histY.GetXaxis()->SetTitle("dy [#mum]");
831  histA.GetYaxis()->SetTitle("Entries");
832  histX.GetYaxis()->SetTitle("Entries");
833  histY.GetYaxis()->SetTitle("Entries");
834  histA.SetTitle("ICP vs ideal, d#alpha");
835  histX.SetTitle("ICP vs ideal, px");
836  histY.SetTitle("ICP vs ideal, py");
837 
838  TCanvas canvas("c1", "c1", 800, 600);
839  canvas.cd();
840  histA.Draw();
841  canvas.Print((pdfPath + "dAlpha.pdf").c_str());
842 
843  histX.Draw();
844  canvas.Print((pdfPath + "dx.pdf").c_str());
845 
846  histY.Draw();
847  canvas.Print((pdfPath + "dy.pdf").c_str());
848 
849 }
850 
851 void checkPairOrientation(){
852 
853 
854  // read pair file
855 
856  // check is sensor orientation is correct
857 
858  // check if dx-dy is weird
859 }
860 
861 int testAlignMatrices() {
862 
863  /*
864  * some example values:
865  *
866  * "/home/arbeit/RedPro3TB/simulationData/2018-04-19-slight-misalignment/himster/m0cm.mat"
867  * "misalignMatrices-SensorsOnly-100u.root"
868  * "/cave_1/lmd_root_0/half_0/plane_0/module_1/sensor_0"
869  */
870 
871  cout << "starting.\n";
872 
873  initDummySimulation();
874  helper = &(PndLmdGeometryHelper::getInstance());
875 
876  //histICPmatrices();
877  //histDairXYZdistances();
878 
879  // TODO: abstract path, cut and 2d/3d to function parameters
880 
881  //compareShiftDataShiftGeo();
882  //return 0;
883 
884  //compareICPmatrices(100, "0");
885 
886  //std::vector<string> cuts { "0", "0.5", "1", "3", "5" };
887 // std::vector<string> cuts { "0"};
888 
889 
890 // for (auto &cut : cuts) {
891 
892 // compareICPmatrices(0, cut);
893 // compareICPmatrices(10, cut);
894 // compareICPmatrices(50, cut);
895 // compareICPmatrices(100, cut);
896 // compareICPmatrices(150, cut);
897 // compareICPmatrices(200, cut);
898 //
899 // buildCyclic(0, cut);
900 // buildCyclic(10, cut);
901 // buildCyclic(50, cut);
902 // buildCyclic(100, cut);
903 // buildCyclic(150, cut);
904 // buildCyclic(200, cut);
905 
906 // }
907 
908  saveMatricesToJson();
909  //saveDetMatricesToJson();
910 
911  // compareMatrices(matrices, "/LMDMatrices/");
912  // temporary fix to avoid double frees at the destruction of te program for pandaroot/fairroot with root6
913  gGeoManager->GetListOfVolumes()->Delete();
914  gGeoManager->GetListOfShapes()->Delete();
915  delete gGeoManager;
916 
917  return 0;
918 }
initDummySimulation
void initDummySimulation()
Definition: testAlignMatrices.C:31
compareShiftDataShiftGeo
void compareShiftDataShiftGeo()
Definition: testAlignMatrices.C:653
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
getOverlapMatrixLikeICP
TGeoHMatrix getOverlapMatrixLikeICP(PndLmdOverlapInfo info)
Definition: testAlignMatrices.C:283
checkPairOrientation
void checkPairOrientation()
Definition: testAlignMatrices.C:851
getMatrixDiff
std::vector< double > getMatrixDiff(TGeoHMatrix &mat1, TGeoHMatrix &mat2)
Definition: testAlignMatrices.C:324
i
Int_t i
Definition: run_full.C:25
PndSdsDetector::SetExclusiveSensorType
void SetExclusiveSensorType(const TString sens)
Definition: PndSdsDetector.cxx:423
baseTransformation
TGeoHMatrix baseTransformation(TGeoHMatrix &input, TGeoHMatrix &toBaseMatrix)
Definition: testAlignMatrices.C:279
exit
exit(0)
PndLmdGeometryHelper::getInstance
static PndLmdGeometryHelper & getInstance()
Definition: PndLmdGeometryHelper.h:125
cuts
TString cuts[MAX]
Definition: autocutx.C:35
readMatrixFromDisk
TGeoHMatrix readMatrixFromDisk(std::string filename)
Definition: testAlignMatrices.C:143
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
matricesMisaligned
std::map< std::string, TGeoHMatrix > * matricesMisaligned
Definition: testAlignMatrices.C:23
applyMisalignmentToGGeoManager
void applyMisalignmentToGGeoManager(TString filename)
Definition: testAlignMatrices.C:58
printMatrixDiff
void printMatrixDiff(TGeoHMatrix mat1, TGeoHMatrix mat2)
Definition: testAlignMatrices.C:252
buildCyclic
void buildCyclic(int alignParam, string cut)
Definition: testAlignMatrices.C:341
getMisalignedOverlapFromGeoManager
TGeoHMatrix getMisalignedOverlapFromGeoManager(PndLmdOverlapInfo &info)
Definition: testAlignMatrices.C:296
PndLmdGeometryHelper::getOverlapInfos
std::vector< PndLmdOverlapInfo > getOverlapInfos(int iHalf=-1, int iPlane=-1, int iModule=-1)
Definition: PndLmdGeometryHelper.cxx:78
json
nlohmann::json json
Definition: testAlignMatrices.C:29
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
cut
double cut[MAX]
Definition: autocutx.C:36
trans
TGeoTranslation * trans
Definition: createROOTEnDrcGeo.C:65
helper
PndLmdGeometryHelper * helper
Definition: testAlignMatrices.C:25
Dipole
FairModule * Dipole
Definition: sim_emc_apd.C:40
saveDetMatricesToJson
void saveDetMatricesToJson()
Definition: testAlignMatrices.C:632
Cave
FairModule * Cave
Definition: sim_emc_apd.C:32
readRootMatrices
std::map< std::string, TGeoHMatrix > * readRootMatrices(TString filename)
Definition: testAlignMatrices.C:231
PndLmdGeometryHelper::getMatrixPndGlobalToSensor
const TGeoHMatrix getMatrixPndGlobalToSensor(const int sensorId)
Definition: PndLmdGeometryHelper.cxx:222
matricesIdeal
std::map< std::string, TGeoHMatrix > * matricesIdeal
Definition: testAlignMatrices.C:24
compareICPmatrices
void compareICPmatrices(int alignParam, string cut)
Definition: testAlignMatrices.C:731
saveMatricesToJson
void saveMatricesToJson()
Definition: testAlignMatrices.C:591
sensor2
TGeoVolume * sensor2
Definition: createTestSensorsTGM.C:61
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
getSmallOverlapInfo
PndLmdOverlapInfo & getSmallOverlapInfo(std::vector< PndLmdOverlapInfo > &infos, int smallOverlap)
Definition: testAlignMatrices.C:313
simOutput
TString simOutput
Definition: Lars/runMvdSim.C:10
ext
TNtuple * ext
Definition: reco_muo.C:24
testAlignMatrices
int testAlignMatrices()
Definition: testAlignMatrices.C:861
rot
TGeoRotation rot
Definition: createRootFscGeometryFile.C:124
unApplyMisalignmentToGGeoManager
void unApplyMisalignmentToGGeoManager(TString filename)
Definition: testAlignMatrices.C:99
json
nlohmann::json json
Definition: convertRootMatricesToJSON.C:5
Pipe
FairModule * Pipe
Definition: sim_emc_apd.C:44
getMatrixSensorToSensor
TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo)
Definition: testAlignMatrices.C:269
Magnet
FairModule * Magnet
Definition: sim_emc_apd.C:36
getMisalignedOverlapFromICP
TGeoHMatrix getMisalignedOverlapFromICP(PndLmdOverlapInfo &info, std::string ICPmatrix)
Definition: testAlignMatrices.C:305
filename
const string filename
Definition: Modify_root_geometry.C:23
PndCave
Definition: PndCave.h:8
Generated on Sun Aug 18 2019 03:02:04 for FairRoot/PandaRoot by   doxygen 1.8.8