PndLmdDim.h

Go to the documentation of this file.

/*
 * PndLmdDim.h
 *
 * this class gives you methods to retrieve parameters of the geometry
 * dimensions are in cm and radian. this is a singleton and must be called
 * with PndLmdDim::Instance()
 *
 * Global comments:
 *
 * TODO
 *
 *
 *  Created on: Oct 5, 2012
 *      Author: promme
 *
 *      To use transformation functions you have to call Read_Transformation_matrices
 *
 *              as of 2015 maintained by Roman Klasen, klasen@kph.uni-mainz.de 

 */

#ifndef PNDLMDDIM_H_
#define PNDLMDDIM_H_

#include <cmath>
#include <TRandom.h>
#include <map>
#include <string>
#include <sstream>
#include <vector>
#include <iostream>
#include <TGeoManager.h>
#include <TGeoVolume.h>
#include <TGeoMatrix.h>
#include <TVector3.h>
#include <TMatrixT.h>
#include <TH2Poly.h>
#include <TMultiGraph.h>
#include <TGraph.h>
#include <stdlib.h>
#include <cctype>

//work with DB
/* #include<PndLmdContFact.h> */
/* #include<TList.h> */
#include <LmdTools/PndLmdAlignPar.h>
#include <TPolyLine.h>
/* #include "FairRuntimeDb.h" */
/* #include "FairRunAna.h" */
/* #include "FairRun.h" */

typedef TMatrixT<double> TMatrixD; // hmm funny, should be declared in TMatrixT.h

using namespace std;

// simple class used for a fast particle track propagation
// based on a matrix fit to simulated data
// the representation of the basis is
// x, y, z, w, px/pz, py/pz, pw
// w und pw are weights in order to
// provide homogeneous coordinates
class PndLmdDimPropMat {
public:
        double mom;
        double M[7][7];
        PndLmdDimPropMat(){
                mom = 0;
                for (int i = 0; i < 7; i++){
                        for (int j = 0; j < 7; j++){
                                M[i][j] = 0;
                        }
                }
        }
        // set a matrix which is given as an array
        void Set(const double* matrix){
                for (int i = 0; i < 7*7; i++){
                        M[i/7][i%7]=matrix[i];
                }
        }
        // set a matrix with fitted hardcoded values
        // up to now only 1.5 GeV/c and 15 GeV/c are available
        // to produce more, please modify and use trafo_matrix_fit.C
        PndLmdDimPropMat(double mom_init){
                mom = mom_init;
                if (mom == 1.5){
                        double matrix_init[49] = {
                                        1.012,  1.22,   5.216e-05,      0.2599,         1.096,  -0.1589,        -0.0001459,
                                        -0.5851,        1.042,  0.001092,       -3.991e-06,     0.3103,         1.064,  -3.991e-06,
                                        -0.04447,       -0.04501,       0.001764,       11.24,  -0.04385,       0.007135,       0.0001733,
                                        -1.36e-15,      1.816e-14,      3.031e-15,      1,      1.513e-16,      -2.272e-16,     7.904e-17,
                                        0.006668,       1.356,  -0.0001811,     0.000127,       0.9794,         -0.1215,        0.4001,
                                        -0.6638,        0.04068,        0.001168,       -7.627e-06,     0.2952,         0.9388,         -7.627e-06,
                                        -1.003e-15,     -3.248e-15,     9.487e-17,      5.594e-18,      2.41e-16,       -1.086e-16,     1
                        };
                        Set(matrix_init);
                } else {
                        if (mom == 15){
                                double matrix_init[49] = {
                                                0.9633,         0.2086,         -0.01173,       0.2598,         1.123,  -0.03217,       -0.0002463,
                                                -0.104,         0.97,   -0.006443,      5.119e-06,      0.06427,        1.12,   5.119e-06,
                                                -0.03828,       0.005853,       0.004478,       11.24,  -0.0447,        0.001794,       2.613e-06,
                                                2.376e-06,      -0.0003806,     -2.583e-07,     1,      -1.647e-05,     -1.871e-05,     -3.862e-07,
                                                -0.01744,       0.1777,         -0.01384,       0.0001758,      0.9997,         -0.02581,       0.4002,
                                                -0.09394,       -0.05056,       -0.002047,      3.624e-05,      0.06091,        0.994,  3.624e-05,
                                                2.376e-06,      -0.0003806,     -2.583e-07,     -3.862e-07,     -1.647e-05,     -1.871e-05,     1
                                };
                                Set(matrix_init);
                        } else { // identity matrix
                                cout << " Warning in PndLmdDimPropMat: loading identity matrix only " << endl;
                                double matrix_init[49] = {
                                                1, 0, 0, 0, 0, 0, 0,
                                                0, 1, 0, 0, 0, 0, 0,
                                                0, 0, 1, 0, 0, 0, 0,
                                                0, 0, 0, 1, 0, 0, 0,
                                                0, 0, 0, 0, 1, 0, 0,
                                                0, 0, 0, 0, 0, 1, 0,
                                                0, 0, 0, 0, 0, 0, 1
                                };
                                Set(matrix_init);
                        }
                }
        }
        // transform the position and momentum at the IP
        // to a position and momentum at the first plane of the LMD
        void Propagate(TVector3& pos, TVector3& momdir){
                //TVector3 dir = mom.Unit();
                double mommag = momdir.Mag();
                // position transformation was evaluated for meter and angles * 10.
                double xip[7] = {pos.X()/100., pos.Y()/100., pos.Z()/100., 1, momdir.X()/momdir.Z()*10., momdir.Y()/momdir.Z()*10., 1};
                double xlmd[7] = {0,0,0,0,0,0,0};
                for (int i = 0; i < 7; i++){
                        for (int j = 0; j < 7; j++){
                                xlmd[i] += M[i][j] * xip[j];
                        }
                }
                pos.SetXYZ(xlmd[0]*100., xlmd[1]*100., xlmd[2]*100.);
                momdir.SetXYZ(xlmd[4]/10., xlmd[5]/10., 1.); // not an exact calculation ;)
                momdir = momdir.Unit()*mommag;
        }
};

class Tkey {
public:
        signed char half;
        signed char plane;
        signed char module;
        signed char side;
        signed char die;
        signed char sensor;
        bool operator < (const Tkey & comp) const{
                if (half < comp.half) return true;
                if (half > comp.half) return false;
                if (plane < comp.plane) return true;
                if (plane > comp.plane) return false;
                if (module < comp.module) return true;
                if (module > comp.module) return false;
                if (side < comp.side) return true;
                if (side > comp.side) return false;
                if      (die < comp.die) return true;
                if      (die > comp.die) return false;
                if      (sensor < comp.sensor) return true;
                if      (sensor >= comp.sensor) return false;
                return false;
        }

        bool operator == (const Tkey & comp) const{
                return (half == comp.half) &&
                                (plane == comp.plane) &&
                                (module == comp.module) &&
                                (side == comp.side) &&
                                (die == comp.die) &&
                                (sensor == comp.sensor);
        }

        Tkey (const Tkey& copy){
                half = copy.half;
                plane = copy.plane;
                module = copy.module;
                side = copy.side;
                die = copy.die;
                sensor = copy.sensor;
        }

        Tkey (int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
                half = ihalf;
                plane = iplane;
                module = imodule;
                side = iside;
                die = idie;
                sensor = isensor;
        }

        Tkey (string key){
                int sign (1);
                int ikey(0); // 0,1,2,3,4,5 = ihalf, iplane, imodule, iside, idie, isensor
                int number;
                for (unsigned int ichar = 0; ichar < key.size(); ichar++){
                        if (key[ichar] == '-'){
                                sign = -1;
                        }
                        if (isdigit(key[ichar])){
                                number = (key[ichar]-'0')*sign;
                                sign = 1;
                                if (ikey == 0) half = number;
                                if (ikey == 1) plane = number;
                                if (ikey == 2) module = number;
                                if (ikey == 3) side = number;
                                if (ikey == 4) die = number;
                                if (ikey == 5) sensor = number;
                                ikey++;
                        }
                }
                if (ikey != 6) cout << " Error in Generate_Tkey: key string " << key << " is not valid " << endl;
        }
        Tkey (){

        }
};

class PndLmdDim {
private:
        // in case you change anything in the geometry
        // you must increase this version number in the corresponding .cxx file!
        // it may affect the consistency between the geometry and the transformation matrices
        static int geometry_version;

        static PndLmdDim* pinstance;
        TGeoManager* fgGeoMan;
        PndLmdDim();
        PndLmdDim(const PndLmdDim &instance);
        PndLmdDim& operator=(const PndLmdDim& instance){
                this->box_size_x = instance.box_size_x; // to get rid from pedantic warnings
                return *this;
        }
        ~PndLmdDim();
        // the navigation paths for the detector geometry are stored here
        vector<string> nav_paths;

public:
        static PndLmdDim& Get_instance();
        static PndLmdDim* Instance();

        // this offset in the sensor id is introduced by
        // pandaroot when assembling several detector components
        // the lmd alone would have a sensor id range between 0 and n total sensors
        // A detector loaded before the lmd will introduce an offset in
        // that ID which can be determined from a loaded root geometry by calling
        // Set_sensIDoffset(-1);
        // TODO: add everywhere that sensID offset
        unsigned int sensIDoffset;

        // set the sensIDoffset
        // in case offset is < 0 the offset is
        // tried to be determined from a possibly loaded root geometry
        bool Set_sensIDoffset(int offset = -1);

        //      FairRun* ana;
        //      FairRuntimeDb* rtdb;
        // pi
        double pi;
        // number of detector planes
        unsigned int n_planes;
        // number of modules per plane half
        unsigned int nmodules;
        // position of planes where the first plane defines the origin
        double* plane_pos_z;
        // ****************************** plane half array *****************************
        // x in lmd reference frame
        double half_offset_x;
        // y in lmd reference frame
        double half_offset_y;
        // z in lmd reference frame
        double half_offset_z;
        // phi rotation in lmd reference frame
        double half_tilt_phi;
        // theta rotation in lmd reference frame
        double half_tilt_theta;
        // psi rotation in lmd reference frame
        double half_tilt_psi;
        // ****************************** plane half supports ******************************
        // x in the detector half reference frame
        double plane_half_offset_x;
        // y in the detector half reference frame
        double plane_half_offset_y;
        // z in the detector half reference frame
        double plane_half_offset_z;
        // phi rotation in the detector half reference frame
        double plane_half_tilt_phi;
        // theta rotation in the detector half reference frame
        double plane_half_tilt_theta;
        // psi rotation in the detector half reference frame
        double plane_half_tilt_psi;

        // ****************************** cvd cooling support discs ************************
        // cvd_diamond is cut out of 79.5 mm discs of 200 micron thickness
        // inner min. radius due to beam pipe + a safety margin
        double inner_rad;
        // not used yet but should be the outer acceptance
        double outer_rad;
        // number of CVD diamond discs per plane
        int n_cvd_discs;
        // radius of a CVD diamond disc
        double cvd_disc_rad;
        // the half of the diamond thickness
        double cvd_disc_thick_half;
        // even and odd discs in a plane will be shifted in z in order to prevent
        // mechanical damage during assembly
        double cvd_disc_even_odd_offset;
        // angle from the division of a circle into n_cvd_discs
        double delta_phi;
        // a polygon of n_cvd_discs sides fitting a radius of inner_rad
        // has a side length pol_side_lg of
        double pol_side_lg_half;
        // the minimum distance to the center of the polygone is given by
        double pol_side_dist_min;
        // the cvd disc has to be placed such that the disc crosses
        // the inner ring at an angle of 0 and delta_phi
        // this defines the distance to the center according to pythagoras
        double cvd_disc_dist;

        // the mechanical alignment precision is defined as an offset of and tilt around
        // the middle of the cvd diamond.
        // Values are standard deviation.
        // first comes translation than rotation

        // x is radial to the beam pipe when support is aligned
        double cvd_offset_x;
        // y is tangent to the beam pipe when support is aligned
        double cvd_offset_y;
        // z is along the beam pipe when support is aligned
        double cvd_offset_z;
        // phi rotation in the reference frame of the plane half support
        double cvd_tilt_phi;
        // theta rotation in the reference frame of the plane half support
        double cvd_tilt_theta;
        // psi rotation in the reference frame of the plane half support
        double cvd_tilt_psi;
        // thickness of the kapton flexible circuits to the sensors
        // (chosen to be thicker to simulate also the influence of the
        //  printed circuits themselves)
        double kapton_disc_thick_half;
        // *********************************** one side on a CVD disc *************************************
        // x is radial to the beam pipe when cvd is aligned
        double side_offset_x;
        // y is tangent to the beam pipe when cvd is aligned
        double side_offset_y;
        // z is along the beam pipe when cvd is aligned
        double side_offset_z; // should not be used due to clashing volumes with the support!
        // phi rotation in the reference frame of the cvd support
        double side_tilt_phi; // should not be used due to clashing volumes with the support!
        // theta rotation in the reference frame of the cvd support
        double side_tilt_theta;
        // psi rotation in the reference frame of the cvd support
        double side_tilt_psi; // should not be used due to clashing volumes with the support!
        // *********************************** HV-MAPS *************************************
        //
        //            left   right
        //
        //  |---------|----------|----------| top
        //  ||------|-||-------|-||-------|-|
        //  ||      | ||       | ||       | |
        //  ||  1   | ||   2   | ||   3   | |     row 1
        //  ||      | || active| ||       | |
        //  ||------|-||-------|-||-------|-|
        //  |         | passive  |          |
        //  |---------|----------|----------| bottom
        //  gap
        //            |----------|----------| top
        //            | passive  |          |
        //            ||-------|-||-------|-|
        //            ||       | ||       | |
        //            ||   2   | ||   3   | |     row 2
        //            || active| ||       | |
        //            ||-------|-||-------|-|
        //            |----------|----------| bottom
        //
        //            left   right
        //
        //                A   y; maps_n_row; height
        //                |
        //                |
        //                --> x; maps_n_col; width
        //
        // the current design foresees a rotation of the first parameters

        // even when several maps are placed on one die
        // those will be placed in the simulation next to
        // each other as separate detectors
        // assumption: sensor 1 and 2 are on one die, the rest is separated
        int maps_n_col;
        int maps_n_row;
        // enabled [row][col]
        bool** enabled;
        // number of sensors per side
        int n_sensors;
        // NOTE: MOST of the following VARIABLES are HALF of it
        // due to geometry construction in GEANT
        double maps_thickness;
        double maps_passive_top;
        double maps_passive_bottom;
        double maps_passive_left;
        double maps_passive_right;
        double maps_active_width;
        double maps_active_height;

        double maps_active_pixel_size; // pixel size NOT half of it

        double maps_width;
        double maps_height;

        double maps_active_offset_x;
        double maps_active_offset_y;


        double die_gap; // (cm)

        double maps_die_width;
        double maps_die_height;

        // the mechanical alignment precision is defined as an offset of and tilt around
        // the middle of the cvd diamond.
        // Values are standard deviation.
        // first comes translation than rotation
        // translations along z as well as rotations around x and y are
        // negligible for dies glued on a cvd diamond
        // rotation around z is not working yet

        // x is along the edge of the cvd disc
        double die_offset_x;
        // y is orthogonal to the edge of the cvd disc
        double die_offset_y;
        // z is along the beam pipe
        double die_offset_z; //should not be used -> crashing volumes;
        // x is a rotation around the edge of the cvd disc
        double die_tilt_phi; // please do not use yet
        // y is a rotation around the orthogonal component of the edge of the cvd disc
        double die_tilt_theta; // please do not use yet
        // z is a rotation around an axis parallel to the along the beam pipe
        double die_tilt_psi;// please do not use yet;
        // *********************************** one sensor *************************************
        // x is mostly radial to the beam pipe
        double sensor_offset_x;
        // y is mostly tangent to the beam pipe
        double sensor_offset_y;
        // z is along the beam pipe
        double sensor_offset_z; // should not be used due to clashing volumes with cvd and flex prints
        // phi rotation in the reference frame of the sensor
        double sensor_tilt_phi;
        // theta rotation in the reference frame of the sensor
        double sensor_tilt_theta; // should not be used due to clashing volumes with cvd and flex prints
        // psi rotation in the reference frame of the sensor
        double sensor_tilt_psi; // should not be used due to clashing volumes with cvd and flex prints
        //*********************************** lumi box parameters ***********************************
        // see CAD files for details
        // https://edms.cern.ch/nav/P:FAIR-000000719:V0/P:FAIR-000000726:V0/TAB3
        // width
        double box_size_x;
        // height
        double box_size_y;
        // length
        double box_size_z;
        // thickness of the V2A steel plates
        double box_thickness;
        // position of the inner rib
        double pos_rib;
        // beam pipe radius at entrance
        double rad_entrance;
        // beam pipe radius at exit
        double rad_exit;
        // beam pipe separating non interacting paricles
        double rad_pipe;
        // beam pipe thickness;
        double pipe_thickness;
        // cone height of the transition region
        double length_transision;
        // length of the inner pipe
        double length_pipe;
        // position of the first detector plane
        double pos_plane_0;
        //*********************************** global parameters *************************************
        // where bending starts with
        double end_seg_upstream;
        // the bending radius
        double r_bend;
        // and the angle of the circle path
        double phi_bend;
        // the point where both tangents of the straight beam pipe tubes meet is
        double pos_rot_z;
        // z position of the lmd
        double pos_z; //(cm)
        double end_seg_bend;
        // x position of the lmd
        double pos_x;
        // y position of the lmd
        double pos_y;
        double rot_x;
        double rot_y;
        double rot_z;

        // returns false when one of the variables exceeds design values
        bool Is_valid_idcall(int ihalf, int iplane = 0, int imodule = 0, int iside = 0, int idie = 0, int isensor = 0){
                if (ihalf   < 0 || ihalf   >= 2)
                        return false;
                if (iplane  < 0 || (unsigned)iplane  >= n_planes)
                        return false;
                if (imodule < 0 || imodule >= n_cvd_discs)
                        return false;
                if (iside   < 0 || iside   >= 2)
                        return false;
                if (idie    < 0 || idie    >= 2)
                        return false;
                // allow to count the non existing inner sensor at die 2
                if (isensor < 0 || isensor >= n_sensors + 1)
                        return false;
                return true;
        }

        // get the sensor id for a sensor on a given side, module and plane
        int Get_sensor_id(int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
                if (idie == 1) isensor += 2; // the parallel sensor to sensor 0 is not there!
                int result = isensor + (iside + (imodule + (iplane + ihalf * n_planes) * nmodules) * 2 ) * n_sensors;
                return result + sensIDoffset;
        }

        // get the sensor position by it's id in terms of plane module and side
        void Get_sensor_by_id(const int sensor_id, int& ihalf, int& iplane, int& imodule, int& iside, int& idie, int& isensor){
                int _sensor_id = sensor_id - sensIDoffset;
                isensor = _sensor_id % n_sensors;
                idie = 0;
                if (isensor > 2) {
                        idie = 1;
                        isensor -= 2; // add the first but non existing sensor at die 2
                }
                _sensor_id /= n_sensors;
                iside = _sensor_id % 2;
                _sensor_id /= 2;
                imodule = _sensor_id % nmodules;
                _sensor_id /= nmodules;
                iplane = _sensor_id % n_planes;
                _sensor_id /= n_planes;
                ihalf = _sensor_id % 2;
                if (!Is_valid_idcall(ihalf, iplane, imodule, iside, idie, isensor)){
                        ihalf = 0;
                        iplane = 0;
                        imodule = 0;
                        iside = 0;
                        idie = 0;
                        isensor = 0;
                        cout << "Error in PndLmdDim::Get_sensor_by_id: "<< sensor_id <<" is not a valid sensor id!" << endl;
                }
        }

        // when calling for a displaced module, the generated
        // value must be kept same for the cvd diamond as well
        // as the sensors sitting on the diamond it self
        // the storage is realized by a map
        // the vector contains offsets in x, y, z, rotphi, rottheta, rotpsi;
        map<Tkey, vector<double> > offsets;
        map<Tkey, vector<double> >::iterator itoffset;
        /* replaced by an ugly but faster version below
        string Generate_key(int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
                stringstream keystream;
                keystream << ihalf << iplane << imodule << iside << idie << isensor;
                return keystream.str();
        }*/

        char* itoa(int value, char* result, int base) {
                // check that the base if valid
                char* last_char;
                if (base < 2 || base > 36) { *result = '\0'; return result; }
                char* ptr = result, *ptr1 = result, tmp_char;
                int tmp_value;
                do {
                        tmp_value = value;
                        value /= base;
                        *ptr++ = "zyxwvutsrqponmlkjihgfedcba9876543210123456789abcdefghijklmnopqrstuvwxyz" [35 + (tmp_value - value * base)];
                } while ( value );
                // Apply negative sign
                if (tmp_value < 0) *ptr++ = '-';
                last_char = ptr;
                *ptr-- = '\0';
                //cout << last_char << endl;
                while(ptr1 < ptr) {
                        tmp_char = *ptr;
                        *ptr--= *ptr1;
                        *ptr1++ = tmp_char;
                }
                return last_char;
        }

        string Generate_key(int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
                char key[100];
                char* ptr;
                ptr = itoa(ihalf, key, 10);
                ptr = itoa(iplane, ptr, 10);
                ptr = itoa(imodule, ptr, 10);
                ptr = itoa(iside, ptr, 10);
                ptr = itoa(idie, ptr, 10);
                ptr = itoa(isensor, ptr, 10);
                string result(key);
                //stringstream keystream;
                //keystream << ihalf << iplane << imodule << iside << idie << isensor;
                return result;
        }

        // generate a unique integer key not same as the string above since there negative numbers are allowed
        // so in case of adding -1's key is not unique!
        // moreover numbers should be kept 1 digit long
        int Generate_keynumber(unsigned int ihalf = 0, unsigned int iplane = 0,unsigned int imodule = 0,
                        unsigned int iside = 0,unsigned int idie = 0,unsigned int isensor = 0){
                stringstream keystream;
                keystream << ihalf << iplane << imodule << iside << idie << isensor;
                int key;
                key = atoi(keystream.str().c_str());
                return key;
        }

        // same structure as for offsets is used for the transformation matrices
        // stored are matrix operations
        // global -> local lumi
        //   key: ihalf = -1 iplane = -1 imodule = -1 iside = -1 idie = -1 isensor = -1
        // local lumi -> local side on cvd disc
        //   key: ihalf >=0 iplane >= 0 imodule >=0 iside = -1 idie = -1 isensor = -1
        // local side on cvd disc -> local sensor
        //   key: all variable
        //map<string, TGeoMatrix* > transformation_matrices;
        //map<string, TGeoMatrix* > transformation_matrices_aligned; // alternative aligned detector description
        //map<string, TGeoMatrix* >::iterator it_transformation_matrices;
        // to increase the performance by a factor of 5-6 a struct is used for the key now
        map<Tkey, TGeoMatrix* > transformation_matrices;
        map<Tkey, TGeoMatrix* > transformation_matrices_aligned; // alternative aligned detector description
        map<Tkey, TGeoMatrix* >::iterator it_transformation_matrices;

        // cleanup some maps containing only references
        void Cleanup();

        // read transformation matrices from a given file
        // aligned and not aligned are two separate maps
        // containing the description of the detector positions
        // if not filename is specified matrices are searched in
        // VMCWORKDIR/input/matrices.txt
        // you may overwrite the version number if necessary
        // VMCWORKDIR/geometry folder is used if no filename is specified
        // ATTENTION! aligned matrices means perfect geometry, set to false if using misaligned geometry
        void Read_transformation_matrices(string filename = "", bool aligned = true, int version_number = geometry_version);

        // write transformation matrices from a given file
        // aligned and not aligned are two separate maps
        // containing the description of the detector positions
        // you may overwrite the version number if necessary
        // version == 0 is reserved for backward compatibility!
        // VMCWORKDIR/geometry folder is used if no filename (e.g. "") is specified
        // warning overwrites existing trafo_matrices_lmd.dat!
        void Write_transformation_matrices(string filename, bool aligned = true, int version_number = geometry_version);

        // read transformation matrices from a loaded geometry
        // aligned and not aligned are two separate maps
        // containing the description of the detector positions
        // the geometry must be loaded otherwise matrices cannot be read
        // version number will be set according to the geometry version number
        // ToDo: multiply also matrices on the way to the key matrices
        //        in case those are not unity matrices
        bool Read_transformation_matrices_from_geometry(bool aligned = true);

        // apply transformation matrices to a loaded geometry
        // aligned and not aligned are two separate maps
        // containing the description of the detector positions
        // the geometry must be loaded otherwise matrices cannot be read
        // version number will be set according to the geometry version number
        // IMPORTANT: you may choose which PndLmdDim matrices you want to use
        // but a ROOT Geometry can always be only aligned. The original
        // matrix stays untouched!
        // TODO: multiply also matrices on the way to the key matrices
        //        in case those are not unity matrices
        // TODO: Find out how to store the aligned geometry as a default
        //        one to pandaroot parameter files
        bool Write_transformation_matrices_to_geometry(bool aligned = true);

        // get a list of sensor paths in the geometry navigation model
        // returns the path to the lmd top volume
        // It is a recursive search, call it once with the default found_lmd variable
        // in case first_call then gGeoManager->CdTop(); is executed to get to the top node
        // of a geometry
        // The geometry must be loaded
        string Get_List_of_Sensors(vector <string>& list_of_sensors, bool found_lmd = false, bool first_call = true);

        // check a list of sensor paths for validity
        // result is true if tests were sucessful
        // offset is the offset in the sensor number which may be
        // not 0 if the geometry was not created in the first place
        bool Test_List_of_Sensors(vector <string> list_of_sensors, int& offset);

        // Get an offset for a volume, if not existent and random
        // a random offset is generated and stored
        // Is used during generation of geometries when calling
        void Get_offset(int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
                        double& x, double& y, double& z,
                        double& rotphi, double& rottheta, double& rotpsi, bool random = false);

        // set an offset for example for the case of existent alignment
        // values, to generate with Get_offset a geometry that matches
        // those offsets
        void Set_offset(int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
                        double x, double y, double z,
                        double rotphi, double rottheta, double rotpsi);

        //read alignment constants from DB 
        //Feb 2013: currently DB is ASCII file
        void Read_DB_offsets(PndLmdAlignPar *lmdalignpar);

        //correct transformation matrices by align
        void Correct_transformation_matrices();
        void reCreate_transformation_matrices();
        // several functions returning the position and orientation of
        // the luminosity detector

        // the global system is defined by PANDA
        // returns the position and rotation of the lmd detector
        // the origin is PANDA interaction point
        // first comes translation than rotation
        void Get_pos_lmd_global(double& x, double& y, double& z,
                        double& rotx, double& roty, double& rotz,
                        bool misaligned = false){
                rotx = rot_x;
                roty = rot_y;// phi_bend;
                rotz = rot_z;
                x = pos_x;
                y = pos_y;
                z = pos_z;
                if (misaligned) rotx += 0; // pedantic compiler fix
        }

        // decoding a digital hit by the position of the sensor given by the sensor ID
        // and the row and column entry
        // parameters of MC geometry input are used to determine the position of
        // the active area within the sensor volume
        // The hit point is returned in the panda global reference frame
        // needs transformation matrices to determine the position
        // column and row can be also the mean from a cluster and therefore
        // not an integer
        TVector3 Decode_hit(const int sensorID, const double column, const double row, const bool aligned = true, bool newVersion=false);

        // matrices for a fast prediction of a particle track to the lmd
        // depending on the momentum setting of the hesr which is the key
        map<double, PndLmdDimPropMat> propagation_matrices;

        // propagate a particle at the IP to the first plane of the lmd
        // the propagation is based on a transformation matrix which was
        // fit to GEANT4 propagated data from January 2015
        // only 1.5 GeV/c and 15 GeV/c are implemented properly
        // TODO: interpolation between matrices
        // back propagation not implemented yet since
        // Mathematica inverted matrices did not work
        void Propagate_fast_ip_to_lmd(TVector3& pos, TVector3& mom, double pbeam);

        // Find the corresponding sensor on the opposite side of a module
        // point is in the panda frame
        // point is transformed into the local lmd frame and
        // the projection is tested 1. to lie on the given sensor
        // 2. on all other sensors on the opposite side for intersection
        // returns true if the corresponding sensor was found
        // in that case input parameters are set to the corresponding sensor
        // otherwise those are set to -1!
        // do not use in time critical cases!
        bool Get_overlapping_sensor(const TVector3& point,int& ihalf, int& iplane, int& imodule, int& iside, int& idie, int& isensor, bool aligned = true);

        // Test if a point lies on a specific sensor
        // The point in the panda frame is transformed into the sensor frame
        // the projection is tested to be within the boundaries of the active sensor area
        bool Is_on_Sensor(const TVector3& point,int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // get the overlapping sensors by geometric constraints
        // hard coded overlapping areas are given sorted by the overlapping area
        // returned is the number of overlapping sensors
        // jdie and jsensor are vectors with the length of the returned number
        int Get_overlapping_sensor(int idie, int isensor, vector<int> &jdie, vector<int> &jsensor);

        // get the transformation matrix from the PANDA global reference frame to the
        // Luminosity reference frame
        TGeoHMatrix Get_transformation_global_to_lmd_local(bool aligned = true);

        // get the transformation matrix from lmd local reference frame to the
        // cvd disc surface reference frame
        TGeoHMatrix Get_transformation_lmd_local_to_module_side(
                        int ihalf, int iplane, int imodule, int iside, bool aligned = true);

        // get the transformation matrix from lmd cvd disc surface frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_module_side_to_sensor(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the transformation matrix from PANDA global reference frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_global_to_sensor(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the transformation matrix from lumi reference frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_lmd_local_to_sensor(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the inverse transformation matrix from the PANDA global reference frame to the
        // Luminosity reference frame
        TGeoHMatrix Get_transformation_lmd_local_to_global(bool aligned = true);

        // get the inverse transformation matrix from lmd local reference frame to the
        // cvd disc surface reference frame
        TGeoHMatrix Get_transformation_module_side_to_lmd_local(
                        int ihalf, int iplane, int imodule, int iside, bool aligned = true);

        // get the inverse transformation matrix from lmd cvd disc surface frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_sensor_to_module_side(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the inverse transformation matrix from PANDA global reference frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_sensor_to_global(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the inverse transformation matrix from lmd local reference frame to the
        // sensor reference frame
        TGeoHMatrix Get_transformation_sensor_to_lmd_local(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);


        // get the transformation matrix from a ideal sensor to the aligned one
        TGeoHMatrix Get_transformation_sensor_to_sensor_aligned(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor);

        // get the transformation matrix from a aligned sensor to the ideal one
        TGeoHMatrix Get_transformation_sensor_aligned_to_sensor(
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor);

        //************************** 3d vector transformations ******************************

        // Transform from the PANDA global reference frame to the
        // Luminosity reference frame
        TVector3 Transform_global_to_lmd_local(const TVector3& point, bool isvector = false, bool aligned = true);

        // Transform from lmd local reference frame to the
        // cvd disc surface reference frame
        TVector3 Transform_lmd_local_to_module_side(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, bool isvector = false, bool aligned = true);

        // Transform from lmd cvd disc surface frame to the
        // sensor reference frame
        TVector3 Transform_module_side_to_sensor(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from PANDA global reference frame to the
        // sensor reference frame
        TVector3 Transform_global_to_sensor(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from lmd local reference frame to the
        // sensor reference frame
        TVector3 Transform_lmd_local_to_sensor(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from the PANDA global reference frame to the
        // Luminosity reference frame
        TVector3 Transform_lmd_local_to_global(const TVector3& point, bool isvector = false, bool aligned = true);

        // Transform from lmd local reference frame to the
        // cvd disc surface reference frame
        TVector3 Transform_module_side_to_lmd_local(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, bool isvector = false, bool aligned = true);

        // Transform from lmd cvd disc surface frame to the
        // sensor reference frame
        TVector3 Transform_sensor_to_module_side(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from PANDA global reference frame to the
        // sensor reference frame
        TVector3 Transform_sensor_to_global(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from lmd local reference frame to the
        // sensor reference frame
        TVector3 Transform_sensor_to_lmd_local(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false, bool aligned = true);

        // Transform from a ideal sensor to the aligned one
        TVector3 Transform_sensor_to_sensor_aligned(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false);

        // Transform from a aligned sensor to the ideal one
        TVector3 Transform_sensor_aligned_to_sensor(const TVector3& point,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool isvector = false);

        // Transform from one sensor to an other
        TVector3 Transform_sensor_to_sensor(const TVector3& point,
                        int ihalf_from, int iplane_from, int imodule_from, int iside_from, int idie_from, int isensor_from,
                        int ihalf_to, int iplane_to, int imodule_to, int iside_to, int idie_to, int isensor_to,
                        bool isvector = false, bool aligned = true);

        // ************************* 3d matrix rotations ***************************************

        // Transform from the PANDA global reference frame to the
        // Luminosity reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_global_to_lmd_local(const TMatrixD& matrix, bool aligned = true);

        // Transform from lmd local reference frame to the
        // cvd disc surface reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_lmd_local_to_module_side(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, bool aligned = true);

        // Transform from lmd local reference frame to the
        // sensor reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_lmd_local_to_sensor(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);


        // Transform from lmd cvd disc surface frame to the
        // sensor reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_module_side_to_sensor(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // Transform from PANDA global reference frame to the
        // sensor reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_global_to_sensor(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // Transform from the PANDA global reference frame to the
        // Luminosity reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_lmd_local_to_global(const TMatrixD& matrix, bool aligned = true);

        // Transform from lmd local reference frame to the
        // cvd disc surface reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_module_side_to_lmd_local(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, bool aligned = true);

        // Transform from lmd cvd disc surface frame to the
        // sensor reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_to_module_side(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // Transform from sensor reference frame to the
        // lmd local reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_to_lmd_local(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // Transform from PANDA global reference frame to the
        // sensor reference frame
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_to_global(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        // Transform from a ideal sensor to the aligned one
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_to_sensor_aligned(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor);

        // Transform from a aligned sensor to the ideal one
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_aligned_to_sensor(const TMatrixD& matrix,
                        int ihalf, int iplane, int imodule, int iside, int idie, int isensor);

        // Transform from one sensor to an other
        // treats only 3 x 3 matrices representing the space
        TMatrixD Transform_sensor_to_sensor(const TMatrixD& matrix,
                        int ihalf_from, int iplane_from, int imodule_from, int iside_from, int idie_from, int isensor_from,
                        int ihalf_to, int iplane_to, int imodule_to, int iside_to, int idie_to, int isensor_to, bool aligned = true);


        // *************************

        // get a pointer to the requested matrices with checks
        // returns NULL if no matrices available
        // do not delete!
        map<Tkey, TGeoMatrix* >* Get_matrices(bool aligned = true);

        // get a pointer to the requested matrix with checks
        // returns NULL if no matrix available
        // do not delete or modify unless you know why you want to
        // load the two files for matrices first
        TGeoMatrix* Get_matrix(int ihalf, int iplane, int imodule, int iside, int idie, int isensor,  bool aligned = true);

        // get the transformation matrix for an path in an existing root geometry
        // no checks are performed in advance
        // result may be 0!
        // if (aligned) the matrix after a possible alignment is returned
        // in that case details to the matrix must be provided in form
        // of ihalf ... isensor
        // TODO: get rid of path and do it only on the basis of ihalf ... isensor
        // if (!aligned) the original matrix is returned
        TGeoHMatrix* Get_matrix(string path, bool aligned = true,
                        int ihalf = -1, int iplane= -1, int imodule = -1, int iside = -1, int idie = -1, int isensor = -1);

        // set the transformation matrix for an path in an existing root geometry
        // A matrix can be only aligned,
        // therefore by default original matrices are not touched!
        // since a key must be created for a node
        // details to it must be provided in form of
        // ihalf ... isensor
        // TODO: get rid of path and do it only on the basis of ihalf ... isensor
        bool Set_matrix(string path, TGeoHMatrix* matrix,
                        int ihalf = -1, int iplane= -1, int imodule = -1, int iside = -1, int idie = -1, int isensor = -1);//, bool aligned = true);

        // get the difference between two matrices of the aligned and misaligned branch
        // in terms of displacement and euler angles
        // all 0 in case of troubles and result is false
        // load the two files for matrices first
        // NOTE: output values are formatted to µm and µrad!!!
        bool Get_matrix_difference(int ihalf, int iplane, int imodule, int iside, int idie, int isensor,
                        double& dx, double& dy, double& dz, double& dphi, double& dtheta, double& dpsi);

        // calculates the differences of matrices between loaded
        // aligned and misaligned ones
        // As an output a table and some histograms are generated
        // load the two files for matrices first
        void Calc_matrix_offsets();

        // see get matrix
        TGeoMatrix* Get_matrix_global_to_lmd_local(bool aligned = true);

        // see get matrix
        TGeoMatrix* Get_matrix_lmd_local_to_module_side(int ihalf, int iplane, int imodule, int iside, bool aligned = true);

        // see get matrix
        TGeoMatrix* Get_matrix_module_side_to_sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true);

        //      TGeoMatrix* Get_matrix_global_to_sector_local(int ihalf, int iplane, int imodule, bool aligned=true);

        // x, y, z coordinate transformation from the PANDA global reference frame to the
        // local reference frame of the luminosity monitor
        void Transform_global_to_lmd_local(double& x, double& y, double& z, bool aligned = true);

        // x, y, z vector transformation from the PANDA global reference frame to the
        // local reference frame of the luminosity monitor
        void Transform_global_to_lmd_local_vect(double& x, double& y, double& z, bool aligned = true);

        // x, y, z coordinate transformation from the local sensor reference frame to the
        // local reference frame of the luminosity monitor
        //void transform_sensor_local_to_lmd_local(const int sensor_id, double& x, double& y, double& z, bool misaligned = false);

        // x, y, z coordinate transformation from the local luminosity reference frame to the
        // global reference frame of panda
        //void transform_local_lmd_to_global(const int sensor_id, double& x, double& y, double& z, bool misaligned = false);

        // x, y, z coordinate transformation between the local sensor reference frames
        // of misaligned and aligned sensors
        //void transform_local_sensor();


        // Generates the luminosity monitor geometry into the mother volume
        // Please make sure that mother volume is large enough or that it
        // is an assembly volume
        // in addition rotation and translation matrices are calculated
        // to store those into a file please use
        // Write_transformation_matrices(filename, false);
        // in case of misaligned, offsets are retrieved via get_offset
        // and included into the geometry
        void Generate_rootgeom(TGeoVolume& mothervol, bool misaligned = false);

        // returns false
        // if version number of the geometry could not be retrieved from a loaded geometry
        // geometry must be available via the root geometry manager
        bool Retrieve_version_number();

        // small function to test some transformation matrices and methods
        void Test_matrices();

        // Draw the Sensors as overlays to an active root pad
        // Projection in XY is used
        // lmd_frame == true : The Lumi reference frame is used
        void Draw_Sensors(int iplane, bool aligned = true, bool lmd_frame = true, int glside = 2); //by default draw both sides
        // Get one Sensor as a polyline to be drawn to an active root pad
        // Projection in XY is used
        // Important -> Do not delete the PolyLine until you are sure that you are
        // finished with displaying it!
        // lmd_frame == true : The Lumi reference frame is used
        TPolyLine* Get_Sensor_Shape(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true, bool lmd_frame = true);

        // Get one Sensor as a vector of graphs
        // Projection in XY is used
        // lmd_frame == true : The Lumi reference frame is used
        // else the sensor reference frame
        // if pixel_subdivision: several graphs are returned
        // representing the passive area and the
        // active area subdivided into several pixels
        // please do not call several times with same
        // parameters if previous result and it's contents
        // are not deleted
        vector<TGraph*> Get_Sensor_Graph(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true, bool lmd_frame = true, bool pixel_subdivision = true);


        // Get a TH2Poly histogram where bins matching the x_y projected shapes
        // of the sensors in one plane
        // The active area can be subdivided further
        // the active pixels: NOT RECOMMENDED due to heavy ram usage
        // Please rename it if you intend to call the function
        // several times with same parameters
        TH2Poly* Get_histogram_Plane(int iplane, int iside, bool aligned =true, bool lmd_frame = true, bool pixel_subdivision = false);

        // Get a TH2Poly histogram for one module side
        // in the reference frame of the lumi
        // on request in the panda reference frame
        // The active area can be subdivided further
        // into the active pixels: WARNING heavy ram usage
        // Please rename it if you intend to call the function
        // several times with same parameters
        TH2Poly* Get_histogram_Moduleside(int ihalf, int iplane, int imodule, int iside, bool aligned = true, bool lmd_frame = true, bool pixel_subdivision = true);

        // Get a TH2Poly histogram for one sensor
        // in the reference frame of the lumi
        // on request in the panda reference frame
        // The active area can be subdivided into
        // and the active pixels
        // Please rename it if you intend to call the function
        // several times with same parameters
        TH2Poly* Get_histogram_Sensor(int ihalf, int iplane, int imodule, int iside, int idie, int isensor, bool aligned = true, bool lmd_frame = true);

        //gets all ids that correspond to an overlapping area
        std::vector<int> getAvailableOverlapIDs();
        int makeOverlapID(int firstSensorId, int secondSensorId);
        int getID1fromOverlapID(int overlapID);
        int getID2fromOverlapID(int overlapID);

        //this is a legacy function and should not be called anymore
        int makeModuleID(int overlapID);

};

#endif /* PNDLMDDIM_H_ */