old-macros/create_Lmd_HV_MAPS.C

Go to the documentation of this file.

// compile by adding the lmd subdirectory to the CMakesLists of the macro folder

#include<TROOT.h>
#include<FairGeoLoader.h>
#include<FairGeoInterface.h>
#include<FairGeoBuilder.h>
#include<FairGeoPcon.h>
#include<FairGeoMedia.h>
#include<TGeoCompositeShape.h>
#include<TGeoMatrix.h>
#include<TGeoVolume.h>
#include<TGeoTube.h>
#include<TGeoArb8.h>
#include<TGeoPara.h>
#include<TSystem.h>
#include<TFile.h>
#include<TGeoPcon.h>
#include<TGeoManager.h>
#include<TGeoTorus.h>
#include<TGeoSphere.h>
#include<TMath.h>
#include<TRandom.h>
#include<PndLmdDim.h>
#include<TGeoPhysicalNode.h>
#include<TGeoCone.h>
#include<TApplication.h>

#include<iostream>
#include<cmath>
#include<sstream>
#include <vector>

//#include<TH2.h>
//#include<TPad.h>

//using namespace LmdDim;
// to check the position and dimensions
// one may enable a beam pipe dummy to be displayed
const bool show_beam_pipe_dummy = true;

const bool include_box = false;

void Navigate_into(int depth = 0){
        //for (int idepth = 0; idepth < depth; idepth++){
        //      cout << "        ";
        //}
        cout << gGeoManager->GetPath() << " ";
        cout << "current node is " << gGeoManager->GetCurrentNavigator()->GetCurrentNode()->GetName() << endl;
        int nnodes = gGeoManager->GetCurrentNode()->GetNdaughters();
        for (int i = 0; i < nnodes; i++){
                cout << " navigating into node " << i << endl;
                gGeoManager->CdDown(i);
                Navigate_into(depth+1);
                gGeoManager->CdUp();
        }
}


int create_HV_MAPS(bool misalign = true) {
        /*
        // ****************************** Parameters of Detector ************************
        //--------------------------------------------------------------------
        gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");

        // Load this libraries
        gSystem->Load("libGeoBase");
        gSystem->Load("libParBase");
        gSystem->Load("libBase");
        gSystem->Load("libPndData");
        gSystem->Load("libPassive");


        TString outfile = "geo/HV_MAPS-Design.root";
        TFile* fi = new TFile(outfile, "RECREATE");

        FairGeoLoader* geoLoad = NULL;// FairGeoLoader::Instance();//
        if (!geoLoad){
                cout << " creating FairGeoLoader instance " << endl;
                geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
        }

        FairGeoInterface *geoFace = geoLoad->getGeoInterface();
        geoFace->setMediaFile("../../geometry/media_pnd.geo");
        geoFace->readMedia();
        geoFace->print();

        FairGeoMedia *Media = geoFace->getMedia();
        FairGeoBuilder *geobuild = geoLoad->getGeoBuilder();

        // retrieve available media
        std::string str_ti = "Aluminum"; // "titanium" is in Panda software not yet implemented
        FairGeoMedium* FairMediumAir = Media->getMedium("air");
        FairGeoMedium* FairMediumSteel = Media->getMedium("steel");
        FairGeoMedium* FairMediumTi = Media->getMedium(str_ti.c_str()); //titanium"); not found in media_pnd.geo !
        FairGeoMedium *FairMediumSilicon = Media->getMedium("silicon");
        FairGeoMedium *FairMediumDiamond = Media->getMedium("HYPdiamond");
        FairGeoMedium *FairMediumVacuum = Media->getMedium("vacuum");
        FairGeoMedium *FairMediumKapton = Media->getMedium("kapton");

        if (!FairMediumAir || !FairMediumSteel || !FairMediumTi || !FairMediumKapton || !FairMediumSilicon || !FairMediumVacuum) {
                std::cout << " warning: not all media found " << std::endl;
        }

        int nmed;
        nmed = geobuild->createMedium(FairMediumAir);
        nmed = geobuild->createMedium(FairMediumSteel);
        nmed = geobuild->createMedium(FairMediumTi);
        nmed = geobuild->createMedium(FairMediumKapton);
        nmed = geobuild->createMedium(FairMediumSilicon);
        nmed = geobuild->createMedium(FairMediumDiamond);
        nmed = geobuild->createMedium(FairMediumVacuum);

        TGeoManager* gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
        if (!gGeoMan) {
                cout << " no FAIR Geometry manager found! " << endl;
                return;
        }

        TGeoVolume *top = new TGeoVolumeAssembly("lmd_HV_MAPS");

        gGeoMan->SetTopVolume(top);
        PndLmdDim& lmddim = PndLmdDim::Get_instance();
         */
        //--------------------------------------------------------------------
        //gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");

        // Load this libraries
        //gSystem->Load("libGeoBase");
        //gSystem->Load("libParBase");
        //gSystem->Load("libBase");
        //gSystem->Load("libPndData");
        //gSystem->Load("libPassive");
        char* dirc = getenv("VMCWORKDIR");
        string dir;
        if(dirc == NULL){
                //env must be set, will crash otherwise
                cout << "warning! environment variable VMCWORKDIR not found! set this variable and try again!\n";
                exit(1);
        }
        else{
                dir = dirc;
        }

        //string moveto = dir + "/input/" + fileName + ".dat";
        //system(("mv "+fileName+".dat"+" "+moveto).c_str());

        TString outfile;
        if (misalign)
                outfile = dir+"/geometry/Luminosity-Detector_misaligned-0u.root";
        else
                outfile = dir+"/geometry/Luminosity-Detector.root";
        TFile* fi = new TFile(outfile,"RECREATE");

        FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
        FairGeoInterface *geoFace = geoLoad->getGeoInterface();
        geoFace->setMediaFile((dir+"/geometry/media_pnd.geo").c_str());
        geoFace->readMedia();
        geoFace->print();

        FairGeoMedia *Media =  geoFace->getMedia();
        FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();

        FairGeoMedium *FairMediumSilicon = Media->getMedium("silicon");

        Int_t nmed=geobuild->createMedium(FairMediumSilicon);

        TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");

        TGeoVolume *top = new TGeoVolumeAssembly("lum");

        gGeoMan->SetTopVolume(top);

        PndLmdDim& lmddim = PndLmdDim::Get_instance();

        // for testing purposes ************************
        // only misaligned sensors, everything else should be perfectly aligned
        lmddim.side_offset_x = 0;
        lmddim.side_offset_y = 0;
        lmddim.side_offset_z = 0;
        lmddim.side_tilt_phi = 0;
        lmddim.side_tilt_theta = 0;
        lmddim.side_tilt_psi = 0;
        lmddim.cvd_offset_x = 0;
        lmddim.cvd_offset_y = 0;
        lmddim.cvd_offset_z = 0;
        lmddim.cvd_tilt_phi = 0;
        lmddim.cvd_tilt_theta = 0;
        lmddim.cvd_tilt_psi = 0;
        lmddim.plane_half_offset_x = 0;
        lmddim.plane_half_offset_y = 0;
        lmddim.plane_half_offset_z = 0;
        lmddim.plane_half_tilt_phi = 0;
        lmddim.plane_half_tilt_theta = 0;
        lmddim.plane_half_tilt_psi = 0;
        lmddim.half_offset_x = 0;
        lmddim.half_offset_y = 0;
        lmddim.half_offset_z = 0;
        lmddim.half_tilt_phi = 0;
        lmddim.half_tilt_theta = 0;
        lmddim.half_tilt_psi = 0;
        // **********************************************
        lmddim.Generate_rootgeom(*top, misalign);
        //lmddim.Write_transformation_matrices("matrices.txt", false);

        if (!misalign)
                lmddim.Write_transformation_matrices("",false); // create default matrices
        else
                //lmddim.Write_transformation_matrices(dir+"/geometry/trafo_matrices_lmd_misaligned.dat", false); // create a misaligned version
                lmddim.Write_transformation_matrices(dir+"/geometry/trafo_matrices_lmd_misaligned-0u.dat", false); // create a misaligned version
        /*
          TH2* hist = lmddim.Get_histogram_Plane(0,0,false,true,false);
          hist->Fill(-5,-5,177);
          hist->Fill(5,5,2);
         // hist->SetMinimum(5);
          hist->Draw("text");
          gPad->Print("hist.pdf");
          gPad->Print("hist.C");*/

        //lmddim.Write_transformation_matrices("", false); // generate standard matrices in VMCWORKDIR/geometry

        //lmddim.Read_transformation_matrices(dir+"geometry/matrices.txt", true);
        //lmddim.Write_transformation_matrices(dir+"geometry/matrices_aligned.txt", true);

        gGeoMan->CloseGeometry();

        gGeoMan->CheckOverlaps(0.001); // [cm]
        //gGeoManager->CheckGeometryFull();
        gGeoMan->PrintOverlaps();
        top->Write();
        fi->Close();
        //gGeoManager->Export(outfile);
        gGeoMan->SetVisLevel(20);
        //lmddim.Retrieve_version_number();
        top->Draw("ogl");

        cout << " found " << gGeoManager->GetListOfPhysicalNodes()->GetEntries() << " entries " << endl;

        if (!misalign){
                cout << " trying to navigate through the geometry " << endl;
                vector<string> list_of_sensors;
                cout << lmddim.Get_List_of_Sensors(list_of_sensors) << endl;

                int offset;

                if (lmddim.Test_List_of_Sensors(list_of_sensors, offset)) cout << " testing of sensors was fine and gave an offset of " << offset << endl;

                //cout << " found " << list_of_sensors.size() << " sensors " << endl;

                //for (int isensor = 0; isensor < list_of_sensors.size() ; isensor++){
                //        cout << list_of_sensors[isensor] << endl;
                //}

                if (1){ // consistency checks
                        cout << " reading matrices from file into the aligned map " << endl;
                        lmddim.Read_transformation_matrices(dir+"/geometry/trafo_matrices_lmd_misaligned.dat", true);
                        cout << " loading transformation matrices from geometry into not aligned map " << endl;
                        lmddim.Read_transformation_matrices_from_geometry(false);
                        cout << " writing matrices from aligned map to geometry" << endl;
                        lmddim.Write_transformation_matrices_to_geometry(true);
                        cout << " reading transformation matrices from geometry once again " << endl;
                        lmddim.Read_transformation_matrices_from_geometry(true);
                        //cout << " testing matrices " << endl;
                        //lmddim.Calc_matrix_offsets();
                }

                if (0){ // test setting individual matrices
                        TGeoHMatrix* matrix = lmddim.Get_matrix("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/lmd_vol_half_1", false, 1,-1,-1,-1,-1,-1);
                        matrix->Print();
                        matrix->RotateX(90.);
                        matrix->Print();
                        lmddim.Set_matrix("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/lmd_vol_half_1", matrix, 1,-1,-1,-1,-1,-1);

                        cout << " the original matrix is " << endl;
                        matrix = lmddim.Get_matrix("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/lmd_vol_half_1", false, 1,-1,-1,-1,-1,-1);
                        matrix->Print();

                        cout << " the aligned matrix is " << endl;
                        matrix = lmddim.Get_matrix("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/lmd_vol_half_1", true, 1,-1,-1,-1,-1,-1);
                        matrix->Print();
                }
                cout << " testing matrices " << endl;
                lmddim.Calc_matrix_offsets();

                gGeoMan->RefreshPhysicalNodes(); // should be called but is not a must
                //gGeoMan->CloseGeometry();
                top->Draw("ogl"); // an already drawn geometry will be not updated according to changes in the matrices

        }
        //gGeoManager->CdTop();
        //Navigate_into();

        /*
          gGeoManager->cd("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/"
                "lmd_vol_half_1/lmd_vol_plane_3/lmd_vol_module_4/"
                "lmd_vol_side_1/lmd_vol_die_1/LumActivePixelRect_398");
          gGeoManager->GetCurrentNode()->GetMatrix();

//gGeoManager->cd("/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0");
          //cout << gGeoManager->GetPath() << " " << gGeoManager->GetCurrentNode()->GetName() << endl;
          //gGeoManager->GetCurrentNode()->GetMatrix()->Print();
          //gGeoManager->GetCurrentMatrix()->Print();


          gGeoManager->cd("/lum_1/lmd_vol_vac_3");
          TGeoVolume* topvol = gGeoManager->GetCurrentNode()->GetVolume();
          if (topvol){
          gGeoManager->SetTopVolume(topvol);
          gGeoManager->cd("/lmd_vol_ref_sys_0");
          cout << gGeoManager->GetPath() << " " << gGeoManager->GetCurrentNode()->GetName() << endl;
          gGeoManager->GetCurrentNode()->GetMatrix()->Print();
          gGeoManager->GetCurrentMatrix()->Print();
         */

        // according to http://personalpages.to.infn.it/~puccio/htmldoc/src/AliTRDalignment.cxx.html
        /*
          TGeoPNEntry *pne;
          TGeoMatrix* ideSm[400];
          TGeoMatrix* misSm[400];
          for (int isensor = 0; isensor < 399; isensor++){
                  stringstream path;
                  path << "/lum_1/lmd_vol_vac_3/lmd_vol_ref_sys_0/lmd_vol_half_1/lmd_vol_plane_3/lmd_vol_module_4/lmd_vol_side_1/lmd_vol_die_1/LumActivePixelRect_" << isensor;
                  gGeoManager->cd(path.str().c_str());
                  cout << gGeoManager->GetCurrentNode()->GetName() << endl;
                  pne = gGeoManager->GetAlignableEntry(path.str().c_str());
                  if (pne){
                          TGeoPhysicalNode* node = pne->GetPhysicalNode();
                          if (!node) {
                                  cout << "physical node entry has no physical node " << path.str() << endl;
                                  continue;
                          }
                          misSm[isensor] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
                          ideSm[isensor] = new TGeoHMatrix(*node->GetOriginalMatrix());
                  } else {
                          cout << " sorry could not get alignable entry from sensor " << isensor << endl;
                  }
          }*/

        //   gGeoManager->Export(outfile);
        //gGeoMan->SetVisLevel(20);
        //top->Draw("ogl");
        //top->RandomPoints(100000);
  return 0;
}

int main() {
        TApplication myapp("myapp", 0, 0);
        cout << "before\n";
        create_HV_MAPS();
        cout << "after\n";
        myapp.Run();
        return 0;
}