KalFitReadGdml.cxx

Go to the documentation of this file.

#include "G4Geo/MdcG4Geo.h"
#include "G4Geo/BesG4Geo.h"
#include "G4Geo/CgemG4Geo.h"
#include "KalFitAlg/KalFitAlg.h"
#include "KalFitAlg/KalFitTrack.h"
#include "G4Material.hh"
#include "G4Tubs.hh"
#include "GDMLProcessor.hh"
#include "G4NistManager.hh"
#include "math.h"
 
void KalFitAlg::setBesFromGdml(void){
    // --- gas Volumes
    int nGas=3;
    TString gasVolumes[3]={"logicalMdc", "logicalWorld", "Cgem_logic"};
    if(useNCGem_>0) nGas=2;
 
    // --- inner DC
    int nIDC=3;
    TString IDCVolumes[3]={"LogicalMdcInnerFilm1", "logicalMdcSegment2", "LogicalMdcInnerFilm0"};
    if(useNCGem_>0) nIDC=0;
    
    // --- CGEM-IT
    int nCgemVolumes=114;
    TString CgemVolumes[114]={
        "Separator_logic_Al2", "Separator_logic_CarFib", "Separator_logic_Al1",  // 3
 
        "Anode_Cu4_logic2", "Anode_Kapton6_logic2", "Anode_Kapton5_logic2", "Anode_Cu3_logic2", "Anode_Epoxy5_logic2", "Anode_Carbonf2_logic2", "Anode_Epoxy4_logic2", "Anode_Honeycomb_logic2", 
        "Anode_Epoxy3_logic2", "Anode_Kapton3_logic2", "Anode_Epoxy2_logic2", "Anode_Kapton2_logic2", "Anode_Epoxy1_logic2", "Anode_Cu2_logic2", "Anode_Kapton1_logic2", "Anode_Cu1_logic2", //16
        "Gap_I_logic2",
        "GemFoil_Cu2_logic2foil2", "GemFoil_Kapton_logic2foil2", "GemFoil_Cu1_logic2foil2", 
        "Gap_T2_logic2",
        "GemFoil_Cu2_logic2foil1", "GemFoil_Kapton_logic2foil1", "GemFoil_Cu1_logic2foil1",
        "Gap_T1_logic2",
        "GemFoil_Cu2_logic2foil0", "GemFoil_Kapton_logic2foil0", "GemFoil_Cu1_logic2foil0",
        "Gap_D_logic2",//13
        "Cathode_Cu2_logic2", "Cathode_Kapton3_logic2", "Cathode_Epoxy3_logic2", "Cathode_Kapton2_logic2", "Cathode_Epoxy2_logic2", "Cathode_Honeycomb_logic2", "Cathode_Epoxy1_logic2", "Cathode_Carbonf_logic2", "Cathode_Epoxy0_logic2",//9
 
        "Anode_Kapton5_logic1", "Anode_Cu3_logic1", "Anode_Epoxy5_logic1", "Anode_Rohacell2_logic1", "Anode_Epoxy4_logic1", "Anode_Kapton4_logic1", "Anode_Epoxy3_logic1", 
        "Anode_Rohacell1_logic1", "Anode_Epoxy2_logic1", "Anode_Kapton2_logic1", "Anode_Epoxy1_logic1", "Anode_Cu2_logic1", "Anode_Kapton1_logic1", "Anode_Cu1_logic1",//14
        "Gap_I_logic1",
        "GemFoil_Cu2_logic1foil2", "GemFoil_Kapton_logic1foil2", "GemFoil_Cu1_logic1foil2", 
        "Gap_T2_logic1",
        "GemFoil_Cu2_logic1foil1", "GemFoil_Kapton_logic1foil1", "GemFoil_Cu1_logic1foil1",
        "Gap_T1_logic1",
        "GemFoil_Cu2_logic1foil0", "GemFoil_Kapton_logic1foil0", "GemFoil_Cu1_logic1foil0",
        "Gap_D_logic1",//13
        "Cathode_Cu2_logic1", "Cathode_Kapton3_logic1", "Cathode_Epoxy4_logic1", "Cathode_Rohacell2_logic1", "Cathode_Epoxy3_logic1", "Cathode_Kapton2_logic1", "Cathode_Epoxy2_logic1", "Cathode_Rohacell1_logic1", "Cathode_Epoxy1_logic1", "Cathode_Kapton1_logic1",//10
 
        "Anode_Kapton5_logic0", "Anode_Cu3_logic0", "Anode_Epoxy5_logic0", "Anode_Carbonf2_logic0", "Anode_Epoxy4_logic0", "Anode_Honeycomb_logic0", "Anode_Epoxy3_logic0", "Anode_Carbonf1_logic0", 
        "Anode_Epoxy2_logic0", "Anode_Kapton2_logic0", "Anode_Epoxy1_logic0", "Anode_Cu2_logic0", "Anode_Kapton1_logic0", "Anode_Cu1_logic0",//14
        "Gap_I_logic0",
        "GemFoil_Cu2_logic0foil2", "GemFoil_Kapton_logic0foil2", "GemFoil_Cu1_logic0foil2", 
        "Gap_T2_logic0",
        "GemFoil_Cu2_logic0foil1", "GemFoil_Kapton_logic0foil1", "GemFoil_Cu1_logic0foil1",
        "Gap_T1_logic0",
        "GemFoil_Cu2_logic0foil0", "GemFoil_Kapton_logic0foil0", "GemFoil_Cu1_logic0foil0",
        "Gap_D_logic0",//13
        "Cathode_Cu2_logic0", "Cathode_Kapton3_logic0", "Cathode_Epoxy3_logic0", "Cathode_Kapton2_logic0", "Cathode_Epoxy2_logic0", "Cathode_Honeycomb_logic0", "Cathode_Epoxy1_logic0", "Cathode_Kapton1_logic0", "Cathode_Cu1_logic0"//9
    };
    if(useNCGem_==0) nCgemVolumes=0;
 
    // --- beam pipe
    int nPipeV=4;
    TString beamPipeVolumes[4]={"logicalouterBe", "logicaloilLayer", "logicalinnerBe", "logicalgoldLayer"};
 
    // --- volumes to construct
    TString volumes[124];
    int iV=0; int nV=nGas+nIDC+nCgemVolumes+nPipeV;
    cout<<"KalFitAlg needs to prepare "<<nV<<" volumes"<<endl;
    for(int i=0; i<nGas; i++, iV++) volumes[iV]=gasVolumes[i];
    for(int i=0; i<nIDC; i++, iV++) volumes[iV]=IDCVolumes[i];
    for(int i=0; i<nCgemVolumes; i++, iV++) volumes[iV]=CgemVolumes[i];
    for(int i=0; i<nPipeV; i++, iV++) volumes[iV]=beamPipeVolumes[i];
 
    // --- G4 from gdml
    MdcG4Geo* aMdcG4Geo = new MdcG4Geo();
    if(useNCGem_>0) CgemG4Geo* aCgemG4Geo = new CgemG4Geo();
    BesG4Geo* aBesG4Geo = new BesG4Geo();
 
    // --- fill materials
    double Z(0.),A(0.),Ionization(0.),Density(0.),Radlen(0.);
    for(iV=0; iV<nV; iV++)// fill materials
    {
        G4LogicalVolume* logicVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume(volumes[iV]));
        G4Material* material = logicVolume->GetMaterial();
        Z = 0.;
        A = 0.;
        for(int i=0; i<material->GetElementVector()->size(); i++){    
            Z += (material->GetElement(i)->GetZ())*
                (material->GetFractionVector()[i]); 
            A += (material->GetElement(i)->GetA())* 
                (material->GetFractionVector()[i]);
        }
        Ionization = material->GetIonisation()->GetMeanExcitationEnergy();
        Density = material->GetDensity()/(g/cm3);
        Radlen = material->GetRadlen();
        if(iV==0) KalFitTrack::mdcGasRadlen_ = Radlen/10.;
        KalFitMaterial kalFitMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
        _BesKalmanFitMaterials.push_back(kalFitMaterial);
        cout<<volumes[iV]<<": Z: "<<Z<<", A: "<<(A/(g/mole))<<", Ionization: "<<(Ionization/eV)<<", Density: "<<Density<<", Radlen: "<<Radlen<<", as material "<<_BesKalmanFitMaterials.size()<<endl;
    }
 
    // --- fill cylinders
    double lastRmin=-100;
    double radius, Rmax, thick, length , z0;
    for(iV=nGas; iV<nV; iV++)// fill cylinders
    {
        G4LogicalVolume* logicVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume(volumes[iV]));
 
        G4Tubs* tub = dynamic_cast<G4Tubs*>(logicVolume->GetSolid());
        radius = tub->GetInnerRadius()/(cm);
        Rmax   = tub->GetOuterRadius()/(cm);
        thick  = tub->GetOuterRadius()/(cm) - tub->GetInnerRadius()/(cm);
        length = 2.0*tub->GetZHalfLength()/(cm);
        z0     = 0.0;
        if(iV>nGas&&lastRmin!=Rmax) {
            cout<<"  lastRmin!=Rmax from "<<volumes[iV-1]<<" to "<<volumes[iV]<<std::endl;
            if(lastRmin-Rmax>0.0001) // 1 micon
            {
                KalFitCylinder aCylinder(&_BesKalmanFitMaterials[1], Rmax, lastRmin-Rmax, length , z0);
                _BesKalmanFitWalls.push_back(aCylinder);
                cout<<"  add a cylinder with "<<lastRmin-Rmax<<" cm air"<<endl;
            }
        }
        KalFitCylinder aCylinder(&_BesKalmanFitMaterials[iV], radius, thick, length , z0);
        _BesKalmanFitWalls.push_back(aCylinder);
        std::cout<<volumes[iV]<<": "<<"radius outer to inner: "<<radius+thick<<" "<<radius<<", length: "<<length<<std::endl;
        cout<<" wall "<<_BesKalmanFitWalls.size()<<"  Z: "<<_BesKalmanFitMaterials[iV].Z()<<endl;
        lastRmin=radius;
    }
 
}//end of setBesFromGdml