FairRoot/PandaRoot
Typedefs | Functions | Variables
testAlignMatrices.C File Reference
#include <TH1D.h>
#include <TCanvas.h>
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoPhysicalNode.h>
#include "PndLmdGeometryHelper.h"
#include <TGeoNode.h>
#include <TFile.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <FairRunSim.h>
#include <iostream>
#include <fstream>
#include <map>
#include <string>
#include <sstream>
#include <PndLmdHitPair.h>
#include "json.hpp"

Go to the source code of this file.

Typedefs

using json = nlohmann::json
 

Functions

void initDummySimulation ()
 
void applyMisalignmentToGGeoManager (TString filename)
 
void unApplyMisalignmentToGGeoManager (TString filename)
 
TGeoHMatrix readMatrixFromDisk (std::string filename)
 
std::map< std::string,
TGeoHMatrix > * 
readRootMatrices (TString filename)
 
void printMatrixDiff (TGeoHMatrix mat1, TGeoHMatrix mat2)
 
TGeoHMatrix getMatrixSensorToSensor (int sensorOne, int sensorTwo)
 
TGeoHMatrix baseTransformation (TGeoHMatrix &input, TGeoHMatrix &toBaseMatrix)
 
TGeoHMatrix getOverlapMatrixLikeICP (PndLmdOverlapInfo info)
 
TGeoHMatrix getMisalignedOverlapFromGeoManager (PndLmdOverlapInfo &info)
 
TGeoHMatrix getMisalignedOverlapFromICP (PndLmdOverlapInfo &info, std::string ICPmatrix)
 
PndLmdOverlapInfogetSmallOverlapInfo (std::vector< PndLmdOverlapInfo > &infos, int smallOverlap)
 
std::vector< double > getMatrixDiff (TGeoHMatrix &mat1, TGeoHMatrix &mat2)
 
void buildCyclic (int alignParam, string cut)
 
void saveMatricesToJson ()
 
void saveDetMatricesToJson ()
 
void compareShiftDataShiftGeo ()
 
void compareICPmatrices (int alignParam, string cut)
 
void checkPairOrientation ()
 
int testAlignMatrices ()
 

Variables

std::map< std::string,
TGeoHMatrix > * 
matricesMisaligned
 
std::map< std::string,
TGeoHMatrix > * 
matricesIdeal
 
PndLmdGeometryHelperhelper
 

Typedef Documentation

using json = nlohmann::json

Definition at line 29 of file testAlignMatrices.C.

Function Documentation

void applyMisalignmentToGGeoManager ( TString  filename)

Definition at line 58 of file testAlignMatrices.C.

References gGeoManager, and TString.

Referenced by buildCyclic().

58  {
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 }
TGeoManager * gGeoManager
const string filename
TGeoHMatrix baseTransformation ( TGeoHMatrix &  input,
TGeoHMatrix &  toBaseMatrix 
)

Definition at line 279 of file testAlignMatrices.C.

Referenced by getOverlapMatrixLikeICP().

279  {
280  return TGeoHMatrix(toBaseMatrix * input * toBaseMatrix.Inverse());
281 }
void buildCyclic ( int  alignParam,
string  cut 
)

Definition at line 341 of file testAlignMatrices.C.

References applyMisalignmentToGGeoManager(), dx, dy, ext, getMatrixDiff(), getMatrixSensorToSensor(), getMisalignedOverlapFromGeoManager(), getMisalignedOverlapFromICP(), PndLmdGeometryHelper::getOverlapInfos(), getSmallOverlapInfo(), i, PndLmdOverlapInfo::id1, PndLmdOverlapInfo::id2, matricesIdeal, matricesMisaligned, readRootMatrices(), sensor2, TString, and unApplyMisalignmentToGGeoManager().

341  {
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) {
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 }
double dy
std::vector< double > getMatrixDiff(TGeoHMatrix &mat1, TGeoHMatrix &mat2)
Int_t i
Definition: run_full.C:25
std::map< std::string, TGeoHMatrix > * matricesMisaligned
void applyMisalignmentToGGeoManager(TString filename)
TGeoHMatrix getMisalignedOverlapFromGeoManager(PndLmdOverlapInfo &info)
std::vector< PndLmdOverlapInfo > getOverlapInfos(int iHalf=-1, int iPlane=-1, int iModule=-1)
double cut[MAX]
Definition: autocutx.C:36
PndLmdGeometryHelper * helper
std::map< std::string, TGeoHMatrix > * readRootMatrices(TString filename)
std::map< std::string, TGeoHMatrix > * matricesIdeal
TGeoVolume * sensor2
double dx
PndLmdOverlapInfo & getSmallOverlapInfo(std::vector< PndLmdOverlapInfo > &infos, int smallOverlap)
TNtuple * ext
Definition: reco_muo.C:24
void unApplyMisalignmentToGGeoManager(TString filename)
TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo)
TGeoHMatrix getMisalignedOverlapFromICP(PndLmdOverlapInfo &info, std::string ICPmatrix)
void checkPairOrientation ( )

Definition at line 851 of file testAlignMatrices.C.

851  {
852 
853 
854  // read pair file
855 
856  // check is sensor orientation is correct
857 
858  // check if dx-dy is weird
859 }
void compareICPmatrices ( int  alignParam,
string  cut 
)

Definition at line 731 of file testAlignMatrices.C.

References dx, dy, ext, filename, getMatrixDiff(), PndLmdGeometryHelper::getOverlapInfos(), getOverlapMatrixLikeICP(), matricesIdeal, matricesMisaligned, readMatrixFromDisk(), readRootMatrices(), and TString.

731  {
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 }
double dy
TGeoHMatrix getOverlapMatrixLikeICP(PndLmdOverlapInfo info)
std::vector< double > getMatrixDiff(TGeoHMatrix &mat1, TGeoHMatrix &mat2)
TGeoHMatrix readMatrixFromDisk(std::string filename)
std::map< std::string, TGeoHMatrix > * matricesMisaligned
std::vector< PndLmdOverlapInfo > getOverlapInfos(int iHalf=-1, int iPlane=-1, int iModule=-1)
double cut[MAX]
Definition: autocutx.C:36
PndLmdGeometryHelper * helper
std::map< std::string, TGeoHMatrix > * readRootMatrices(TString filename)
std::map< std::string, TGeoHMatrix > * matricesIdeal
double dx
TNtuple * ext
Definition: reco_muo.C:24
const string filename
void compareShiftDataShiftGeo ( )

Definition at line 653 of file testAlignMatrices.C.

References cut, cuts, ext, filename, getMatrixDiff(), PndLmdGeometryHelper::getOverlapInfos(), and readMatrixFromDisk().

653  {
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 }
std::vector< double > getMatrixDiff(TGeoHMatrix &mat1, TGeoHMatrix &mat2)
TString cuts[MAX]
Definition: autocutx.C:35
TGeoHMatrix readMatrixFromDisk(std::string filename)
std::vector< PndLmdOverlapInfo > getOverlapInfos(int iHalf=-1, int iPlane=-1, int iModule=-1)
double cut[MAX]
Definition: autocutx.C:36
PndLmdGeometryHelper * helper
TNtuple * ext
Definition: reco_muo.C:24
const string filename
std::vector<double> getMatrixDiff ( TGeoHMatrix &  mat1,
TGeoHMatrix &  mat2 
)

Definition at line 324 of file testAlignMatrices.C.

References dx, and dy.

Referenced by buildCyclic(), compareICPmatrices(), and compareShiftDataShiftGeo().

324  {
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 }
double dy
double dx
TGeoHMatrix getMatrixSensorToSensor ( int  sensorOne,
int  sensorTwo 
)

Definition at line 269 of file testAlignMatrices.C.

References PndLmdGeometryHelper::getMatrixPndGlobalToSensor().

Referenced by buildCyclic(), getMisalignedOverlapFromGeoManager(), getMisalignedOverlapFromICP(), and getOverlapMatrixLikeICP().

269  {
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 }
PndLmdGeometryHelper * helper
const TGeoHMatrix getMatrixPndGlobalToSensor(const int sensorId)
TGeoHMatrix getMisalignedOverlapFromGeoManager ( PndLmdOverlapInfo info)

Definition at line 296 of file testAlignMatrices.C.

References getMatrixSensorToSensor(), getOverlapMatrixLikeICP(), PndLmdOverlapInfo::id1, and PndLmdOverlapInfo::id2.

Referenced by buildCyclic().

296  {
297 
298  TGeoHMatrix sen0to1 = getMatrixSensorToSensor(info.id1, info.id2);
299  TGeoHMatrix sen0to1ICPcorr = getOverlapMatrixLikeICP(info);
300  // this order is important!
301  return sen0to1ICPcorr * sen0to1;
302 }
TGeoHMatrix getOverlapMatrixLikeICP(PndLmdOverlapInfo info)
TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo)
TGeoHMatrix getMisalignedOverlapFromICP ( PndLmdOverlapInfo info,
std::string  ICPmatrix 
)

Definition at line 305 of file testAlignMatrices.C.

References getMatrixSensorToSensor(), PndLmdOverlapInfo::id1, PndLmdOverlapInfo::id2, and readMatrixFromDisk().

Referenced by buildCyclic().

305  {
306 
307  TGeoHMatrix sen0to1 = getMatrixSensorToSensor(info.id1, info.id2);
308  TGeoHMatrix sen0to1ICPcorr = readMatrixFromDisk(ICPmatrix);
309  // this order is important!
310  return sen0to1ICPcorr * sen0to1;
311 }
TGeoHMatrix readMatrixFromDisk(std::string filename)
TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo)
TGeoHMatrix getOverlapMatrixLikeICP ( PndLmdOverlapInfo  info)

Definition at line 283 of file testAlignMatrices.C.

References baseTransformation(), getMatrixSensorToSensor(), PndLmdOverlapInfo::id1, PndLmdOverlapInfo::id2, PndLmdOverlapInfo::path1, and PndLmdOverlapInfo::path2.

Referenced by compareICPmatrices(), and getMisalignedOverlapFromGeoManager().

283  {
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 }
TGeoHMatrix baseTransformation(TGeoHMatrix &input, TGeoHMatrix &toBaseMatrix)
TGeoHMatrix getMatrixSensorToSensor(int sensorOne, int sensorTwo)
PndLmdOverlapInfo& getSmallOverlapInfo ( std::vector< PndLmdOverlapInfo > &  infos,
int  smallOverlap 
)

Definition at line 313 of file testAlignMatrices.C.

References PndLmdOverlapInfo::overlapID.

Referenced by buildCyclic().

313  {
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 }
void initDummySimulation ( )

Definition at line 31 of file testAlignMatrices.C.

References Cave, Dipole, fRun, Magnet, Pipe, PndSdsDetector::SetExclusiveSensorType(), simOutput, and TString.

Referenced by testAlignMatrices().

31  {
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 }
void SetExclusiveSensorType(const TString sens)
FairRunAna * fRun
Definition: hit_dirc.C:58
FairModule * Dipole
Definition: sim_emc_apd.C:40
FairModule * Cave
Definition: sim_emc_apd.C:32
TString simOutput
FairModule * Pipe
Definition: sim_emc_apd.C:44
FairModule * Magnet
Definition: sim_emc_apd.C:36
Definition: PndCave.h:8
void printMatrixDiff ( TGeoHMatrix  mat1,
TGeoHMatrix  mat2 
)

Definition at line 252 of file testAlignMatrices.C.

References dx, and dy.

252  {
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 }
double dy
double dx
TGeoHMatrix readMatrixFromDisk ( std::string  filename)

Definition at line 143 of file testAlignMatrices.C.

References exit(), filename, i, rot, and trans.

Referenced by compareICPmatrices(), compareShiftDataShiftGeo(), and getMisalignedOverlapFromICP().

143  {
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 }
Int_t i
Definition: run_full.C:25
exit(0)
TGeoTranslation * trans
TGeoRotation rot
const string filename
std::map<std::string, TGeoHMatrix>* readRootMatrices ( TString  filename)

Definition at line 231 of file testAlignMatrices.C.

Referenced by buildCyclic(), and compareICPmatrices().

231  {
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 }
const string filename
void saveDetMatricesToJson ( )

Definition at line 632 of file testAlignMatrices.C.

References gGeoManager.

632  {
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 }
TGeoManager * gGeoManager
void saveMatricesToJson ( )

Definition at line 591 of file testAlignMatrices.C.

References PndLmdGeometryHelper::getOverlapInfos().

Referenced by testAlignMatrices().

591  {
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 }
std::vector< PndLmdOverlapInfo > getOverlapInfos(int iHalf=-1, int iPlane=-1, int iModule=-1)
PndLmdGeometryHelper * helper
nlohmann::json json
int testAlignMatrices ( )

Definition at line 861 of file testAlignMatrices.C.

References PndLmdGeometryHelper::getInstance(), gGeoManager, initDummySimulation(), and saveMatricesToJson().

861  {
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 
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 
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 }
void initDummySimulation()
static PndLmdGeometryHelper & getInstance()
TGeoManager * gGeoManager
PndLmdGeometryHelper * helper
void saveMatricesToJson()
void unApplyMisalignmentToGGeoManager ( TString  filename)

Definition at line 99 of file testAlignMatrices.C.

References gGeoManager, and TString.

Referenced by buildCyclic().

99  {
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 }
TGeoManager * gGeoManager
const string filename

Variable Documentation

std::map<std::string, TGeoHMatrix>* matricesIdeal

Definition at line 24 of file testAlignMatrices.C.

Referenced by buildCyclic(), and compareICPmatrices().

std::map<std::string, TGeoHMatrix>* matricesMisaligned

Definition at line 23 of file testAlignMatrices.C.

Referenced by buildCyclic(), and compareICPmatrices().