CgemGeomSvc.cxx

Go to the documentation of this file.

#include "GaudiKernel/Algorithm.h"
#include "CgemGeomSvc/CgemGeomSvc.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentListener.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "CalibData/CalibModel.h"
#include "GaudiKernel/MsgStream.h"
#include "EventModel/EventHeader.h"
#include "EventModel/Event.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include <math.h>
#include <fstream>
#include <string>
#include <float.h>
 
 
using namespace std;
 
CgemGeomSvc::CgemGeomSvc(const std::string& name, ISvcLocator* pSvcLocator)
     :Service(name,pSvcLocator)
{
     declareProperty("UseIniGeomFile", m_iFile = 6); // CHECK it must be read from the file!
     declareProperty("PrintDebugInfo", m_debugCout = false); 
     declareProperty("TestReadoutPlane", m_testReadoutPlane = false); 
     declareProperty("AlignFile", m_alignFile);
     declareProperty("UseEffectiveDensityHoles", m_effDensityHoles = false);
     declareProperty("UseEffectiveDensityStrips", m_effDensityStrips = false);
}
 
CgemGeomSvc::~CgemGeomSvc(){
    for(int i=0; i<m_N_CgemLayer; i++){
        for(int j=0; j < m_N_GemFoil; j++) delete m_foil[i][j];
        int nSheet=m_CgemLayer[i]->getNumberOfSheet();
        for(int j=0; j<nSheet; j++) delete m_ReadoutPlane[i][j];
        delete m_CgemLayer[i];
    }
    delete m_CgemSeparator;
    delete m_align;
    delete m_midDrift;
 
} 
 
StatusCode CgemGeomSvc::queryInterface (const InterfaceID& riid, void** ppvInterface ){
 
    if ( IID_ICgemGeomSvc.versionMatch(riid) ) { 
        *ppvInterface = static_cast<ICgemGeomSvc*> (this); 
    } else {
        return Service::queryInterface(riid, ppvInterface) ; 
    }   
    return StatusCode::SUCCESS;
}   
 
 
StatusCode CgemGeomSvc::initialize(){
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "initializing......" << endreq;
 
    StatusCode sc = Service::initialize();
    if ( sc.isFailure() ) return sc;
    
    initGeom();
    m_align = new CgemGeoAlign();
    m_align->initAlignPar(m_alignFile);
    for(int layer=0; layer<m_N_CgemLayer; layer++){
        double r = m_CgemLayer[layer]->getMiddleROfGapD();
        m_align->setLayerR(layer, r);
    }
 
    cout << "CgemGeomSvc::initialize(): alignfile: " << m_alignFile << "  --------------------------------------" << endl;
 
    m_midDrift = new CgemMidDriftPlane();
    m_midDrift->setAlignment(m_align);
    for(int layer=0; layer<m_N_CgemLayer; layer++){
        //double rInner = m_CgemLayer[layer]->getInnerROfGapD();
        //double rOuter = m_CgemLayer[layer]->getOuterROfGapD();
        //double r = (rInner + rOuter)*0.5;
        double  r   =   m_CgemLayer[layer]->getMiddleROfGapD();
        m_midDrift->setR(layer, r);
 
        int nSheet=m_CgemLayer[layer]->getNumberOfSheet();
        for(int i=0; i<nSheet; i++){
            m_midDrift->setReadoutPlane(layer, i, m_ReadoutPlane[layer][i]);
        }
 
        // double vAngle = m_CgemLayer[layer]->getAngleOfStereo();
        // m_midDrift->setVangle(layer, vAngle);
    }
 
    if(m_testReadoutPlane) TestReadoutPlane();
 
    return StatusCode::SUCCESS;
}
 
 
StatusCode CgemGeomSvc::finalize(){
 
    MsgStream log(msgSvc(),name());
 
    log<<MSG::INFO << "finalizing....   "<<endreq;
 
    log<<MSG::INFO << "End of Run   "<<endreq;
    return StatusCode::SUCCESS;
 
}
 
void CgemGeomSvc::handle(const Incident& inc){
    MsgStream log( messageService(), name() );
    log << MSG::DEBUG << "handle: " << inc.type() << endreq;
 
    IDataProviderSvc* m_eventSvc;
    Gaudi::svcLocator()->service("EventDataSvc", m_eventSvc, true);
    SmartDataPtr<Event::EventHeader> eventHeader(m_eventSvc,"/Event/EventHeader");
    if (!eventHeader) { 
        log << MSG::FATAL << "Could not find Event Header" << endreq;
    }
//      if (m_updatealign) return;
//      if (inc.type() == "NewRun" ){
//    log << MSG::DEBUG << "Begin Event" << endreq;
//    clean();
//    m_updatealign = true;
//    if(m_nomcalignment&&m_mindex==0) {  
//         int   RunNo=eventHeader->runNumber();
//         if(RunNo<0)  m_readAlignParDataBase=false ;
//         else   m_readAlignParDataBase=true;
//         m_mindex+=1;
//         cout<<"m__RunNo="<<RunNo<<"m_mindex="<<m_mindex<<endl;
//    }
//    //std::cout<<"############"<<m_readAlignParDataBase<<std::endl;
//    ReadAlignPar();          
//      }
}
 
void CgemGeomSvc::readAlignPar(){
}
 
bool CgemGeomSvc::initGeom(){
    cout.precision(10);
    // cout<<"precision: "<<cout.precision()<<endl;
    /* Get the geometry parameter file from asigned path and value */
    string geoFilePath = getenv("CGEMGEOMSVCROOT");
    // int iFile=iPar
    cout << "IFILE " << m_iFile << endl;
    //      if(m_iFile==0) m_iFile=ReadBoostRoot::GetCgem();
    switch (m_iFile){
        case 1  :
            geoFilePath += "/dat/cgem_model_1_default.txt";
            break;
        case 2 :
            geoFilePath += "/dat/cgem_model_2_default.txt";
            break;
        case 3 :
            geoFilePath += "/dat/cgem_model_3_default.txt";
            break;
        case 4 :
            geoFilePath += "/dat/cgem_model_4_default.txt";
            break;
        case 5 :
            geoFilePath += "/dat/cgem_model_5_default.txt";
            break;
        case 6 :
            geoFilePath += "/dat/cgem_model_6_default.txt";
            break;
 
        case 7 :
            geoFilePath += "/dat/cgem_model_oldgeom.txt";
            break;
 
        default :
            cout << "Warning in CgemGeomSvc::initGeom(): The geometry will be constructed through default model!" << endl;
            geoFilePath += "/dat/cgem_model_6_default.txt";
            break;
    }
    cout << "INFO in CgemGeomSvc::initGeom(): Geometry will be constructed through File: " 
        << geoFilePath << endl;
 
    /* Open and read the geometry parameter file */
    ifstream fin(geoFilePath.c_str(), ios::in);
    if (!fin.good())
    {
        cout << "Error in CgemGeomSvc::initGeom(): Fail to open geomtry file: " << endl;
        cout << geoFilePath << endl;
        return false;
    }
    string strcom; /* One input line of fin file */
 
 
    if(m_debugCout) cout << left << "=======get BesMdcGeoParameter from CgemGeomSvc=======" << endl;
 
    /* Information of CgemLayer  and Readout Strip */
    std::getline(fin, strcom); // Line 2
    std::getline(fin, strcom); // Line 3
    fin  >> m_N_CgemLayer      
        >> m_N_GemFoil;      // Line 3
    if(m_debugCout) cout << "N_CgemLayer : " << left  << setw(8) << m_N_CgemLayer 
        << "N_GemFoil : "   << left  << setw(8) << m_N_GemFoil   << endl;
    std::getline(fin, strcom); // Line 4
    std::getline(fin, strcom); // Line 5
    std::getline(fin, strcom); // Line 6
    double rLayer, lenLayer, numSheet, wSheet, aStereo, wPitchX, wPitchV, wStripX, wStripV, nChannelX, nChannelV;
    double phi_start[3][2], dX_start_strip[3][2], dV_start_strip[3][2];
    double gap_d[3];
    for (int i=0; i < m_N_CgemLayer; i++)
    {
        fin  >> rLayer >> lenLayer >> numSheet >> wSheet >> aStereo >> wPitchX >> wPitchV 
            >> wStripX >> wStripV >> nChannelX >> nChannelV 
            >> phi_start[i][0] >> phi_start[i][1] 
            >> dX_start_strip[i][0] >> dX_start_strip[i][1] 
            >> dV_start_strip[i][0] >> dV_start_strip[i][1] 
            >> gap_d[i] >> m_Orientation;  // Line 6 7 8
        //cout<<"rLayer = "<<rLayer<<endl;
        //nChannelX=nChannelX/numSheet;
        //nChannelV=nChannelV/numSheet;
        m_T_Gap_D[i] = gap_d[i];
        m_CgemLayer[i] = new CgemGeoLayer();
        m_CgemLayer[i]->setInnerROfCgemLayer(rLayer);
        m_CgemLayer[i]->setLengthOfCgemLayer(lenLayer);
        m_CgemLayer[i]->setNumberOfSheet(numSheet);
        m_CgemLayer[i]->setWidthOfSheet(wSheet);
        m_CgemLayer[i]->setAngleOfStereo(aStereo);
        m_CgemLayer[i]->setWidthOfPitchX(wPitchX);
        m_CgemLayer[i]->setWidthOfPitchV(wPitchV);
        m_CgemLayer[i]->setWidthOfStripX(wStripX);
        m_CgemLayer[i]->setWidthOfStripV(wStripV);
        m_CgemLayer[i]->setNumberOfChannelX(nChannelX);
        m_CgemLayer[i]->setNumberOfChannelV(nChannelV);
        m_CgemLayer[i]->setOrientation(m_Orientation);
        std::getline(fin, strcom); // Line 7 8 9
    } 
 
    /* Information of Cathode*/
    double ncMaterials[3];
    double tc_cu1[3],tc_kapton1[3], tc_epoxy0[3],tc_carbonf[3],tc_epoxy1[3],tc_honeycomb[3],tc_rohacell1[3],tc_epoxy2[3],tc_kapton2[3],tc_epoxy3[3],tc_rohacell2[3],tc_epoxy4[3],tc_kapton3[3],tc_cu2[3];
    std::getline(fin, strcom); // Line 10
 
 
    //******** L3 new modification from here ******///
 
    fin >>nmc_Materials[0]>>nmc_Materials[1]>>nmc_Materials[2]>>nmc_Materials[3]>>nmc_Materials[4]>>nmc_Materials[5]
        >>nmc_Materials[6]>>nmc_Materials[7]>>nmc_Materials[8]>>nmc_Materials[9]>>nmc_Materials[10]>>nmc_Materials[11]
        >>nmc_Materials[12]>>nmc_Materials[13]>>nmc_Materials[14];  //Line 10
    for (int i=0; i < m_N_CgemLayer; i++)
    {
        std::getline(fin,strcom);   //Line 11   12  13
        fin >>ncMaterials[i]>>tc_cu1[i]>>tc_kapton1[i]>>tc_epoxy0[i]>>tc_carbonf[i]>>tc_epoxy1[i]>>tc_honeycomb[i]>>tc_rohacell1[i]
            >>tc_epoxy2[i]>>tc_kapton2[i]>>tc_epoxy3[i]>>tc_rohacell2[i]>>tc_epoxy4[i]>>tc_kapton3[i]>>tc_cu2[i];
        if ( ncMaterials[i] > N_MATERIAL_MAX) {
            cout << "Error : CgemGeomSvc::initGeom(), Number of materials in Cathode[Layer "<<i<<"] is wrong!" << endl;
            cout << "It must be less than " << N_MATERIAL_MAX << " and is " << ncMaterials[i] << endl; 
            cout << "please, either increase N_MATERIAL_MAX or decrease m_N_Cathode1_Materials" << endl;
            return false;
        }
    }
    /* Information of GemFoil*/
    double ngMaterials[3];
    double tg_cu1[3],tg_kapton[3],tg_cu2[3],r_i_hole[3],r_o1_hole[3],r_o2_hole[3],l_hole[3];
    double tg_gt1[3],tg_gt2[3],tg_gi[3];
    std::getline(fin,strcom); //Line 14
    std::getline(fin,strcom); //Line 15
    for (int i=0;i<m_N_CgemLayer;i++)
    {
        std::getline(fin,strcom);//Line 16  17  18
        fin>>tg_gt1[i]>>tg_gt2[i]>>tg_gi[i];
    }
    std::getline(fin,strcom);   //line 19
    std::getline(fin,strcom);   //line 20
    fin >>nmg_Materials[0]>>nmg_Materials[1]>>nmg_Materials[2]>>nmg_Materials[3];
    for (int i=0;i<m_N_CgemLayer;i++)
    {
        std::getline(fin,strcom);   //line 21   22  23
        fin>>ngMaterials[i]>>tg_cu1[i]>>tg_kapton[i]>>tg_cu2[i]>>r_i_hole[i]>>r_o1_hole[i]>>r_o2_hole[i]>>l_hole[i];
    }
 
 
    /* Information of Anode*/
    double naMaterials[3];
    double ta_cu1[3],ta_kapton1[3],ta_cu2[3],ta_epoxy1[3],ta_kapton2[3],ta_epoxy2[3],ta_rohacell1[3];
    double ta_carbonf1[3], ta_kapton3[3], ta_epoxy3[3], ta_honeycomb[3],ta_kapton4[3],ta_epoxy4[3],ta_rohacell2[3];
    double ta_carbonf2[3],ta_epoxy5[3],ta_cu3[3],ta_kapton5[3], ta_kapton6[3], ta_cu4[3];
    std::getline(fin,strcom); //line 24
    std::getline(fin,strcom); //line 25
    fin >>nma_Materials[0]>>nma_Materials[1]>>nma_Materials[2]>>nma_Materials[3]>>nma_Materials[4]>>nma_Materials[5]
        >>nma_Materials[6]>>nma_Materials[7]>>nma_Materials[8]>>nma_Materials[9]>>nma_Materials[10]>>nma_Materials[11]
        >>nma_Materials[12]>>nma_Materials[13]>>nma_Materials[14]>>nma_Materials[15]>>nma_Materials[16]>>nma_Materials[17]>>nma_Materials[18]
        >>nma_Materials[19]>>nma_Materials[20];
    for (int i=0;i<m_N_CgemLayer;i++)
    {
        std::getline(fin,strcom); //line 26 27 28
        fin >>naMaterials[i]>>ta_cu1[i]>>ta_kapton1[i]>>ta_cu2[i]>>ta_epoxy1[i]>>ta_kapton2[i]>>ta_epoxy2[i]>>ta_rohacell1[i]
            >>ta_carbonf1[i]>>ta_kapton3[i]>>ta_epoxy3[i]>>ta_honeycomb[i]>>ta_kapton4[i]>>ta_epoxy4[i]>>ta_rohacell2[i]
            >>ta_carbonf2[i]>>ta_epoxy5[i]>>ta_cu3[i]>>ta_kapton5[i]>>ta_kapton6[i]>>ta_cu4[i]
            >>m_B_strip_x_on>>m_B_strip_v_on;
        if (naMaterials[i] > N_MATERIAL_MAX) {
            cout << "Error : CgemGeomSvc::initGeom(), Number of materials in Anode[L"<<i<<"] is wrong!" << endl;
            cout << "It must be less than " << N_MATERIAL_MAX << " and is " << naMaterials[i] << endl; 
            cout << "please, either increase N_MATERIAL_MAX or decrease N_Materials" << endl;
            return false;
        }
    }
 
    /** cable density (g/cm3) **/
    std::getline(fin, strcom); // Line 29
    std::getline(fin, strcom); // Line 30 
    fin >> m_Rho_Cable; // Line 30                                                                                                
    if(m_debugCout) cout << "cable density " << m_Rho_Cable << endl;
 
    /** CGEM/MDC separator **/
    std::getline(fin, strcom); // Line 31
    std::getline(fin, strcom); // Line 32
    std::getline(fin, strcom); // Line 33
    double r_i_separator,r_o_separator,l_separator;
    double t_al1,t_al2,t_carbonf; 
    fin >> r_i_separator >> r_o_separator >> l_separator; // Line 33
 
    std::getline(fin, strcom); // Line 34
    std::getline(fin, strcom); // Line 35
    fin >> t_al1 >>  t_carbonf>>  t_al2; // Line 35
 
    m_CgemSeparator = new CgemGeoSeparator();
    m_CgemSeparator->setInnerR(r_i_separator);
    m_CgemSeparator->setOuterR(r_o_separator);
    m_CgemSeparator->setLength(l_separator);
    m_CgemSeparator->setThickOfInnerAluminum(t_al1); 
    m_CgemSeparator->setThickOfCarbonFiber(t_carbonf);
    m_CgemSeparator->setThickOfOuterAluminum(t_al2);
 
    /* Information of Total volume in simulation */
 
    std::getline(fin, strcom); // Line 36
    std::getline(fin, strcom); // Line 37
    std::getline(fin, strcom); // Line 38
    fin >> m_R_i_Cgem >> m_R_o_Cgem >> m_L_Cgem; // Line 38
 
    std::getline(fin, strcom); // Line 39
    std::getline(fin, strcom); // Line 40
    std::getline(fin, strcom); // Line 41
    fin >> m_passive; // Line 41
    if(m_debugCout) cout << "Cgem passive elements are switched " << m_passive << endl;
 
    double T_Cathode,T_Anode, T_GemFoil;
    for (int i=0; i < m_N_CgemLayer; i++) {
        T_Cathode   =   tc_cu1[i]+tc_kapton1[i]+tc_epoxy0[i]+tc_carbonf[i]+tc_epoxy1[i]+tc_honeycomb[i]+tc_rohacell1[i]+tc_epoxy2[i]+tc_kapton2[i]+tc_epoxy3[i]+tc_rohacell2[i]+tc_epoxy4[i]+tc_kapton3[i]+tc_cu2[i];
        m_T_Cathode[i] = T_Cathode;
 
        T_Anode =   ta_cu1[i]+ta_kapton1[i]+ta_cu2[i]+ta_epoxy1[i]+ta_kapton2[i]+ta_epoxy2[i]+ta_rohacell1[i]+ta_carbonf1[i]+ta_kapton3[i]+ta_epoxy3[i]+ta_honeycomb[i]+ta_kapton4[i]+ta_epoxy4[i]+ta_rohacell2[i]+ta_carbonf2[i]+ta_epoxy5[i]+ta_cu3[i]+ta_kapton5[i]+ta_kapton6[i]+ta_cu4[i];
        m_T_Anode[i] = T_Anode;
 
        T_GemFoil   =   tg_cu1[i]+tg_kapton[i]+tg_cu2[i];
        m_T_CgemLayer[i] = T_Cathode+ gap_d[i] + tg_gt1[i] + tg_gt2[i] + tg_gi[i] + T_GemFoil*3 + T_Anode;
        /* output the thickness information of each material layer */
        if(m_debugCout) 
        {
            cout    << left << "R_CgemLayer " << left << "L_CgemLayer "
                << left << "N_Sheet " << left << "W_Sheet "
                << left << "A_Stereo " << left << "W_Pitch_X "
                << left << "W_Pitch_V "<< left << "W_Strip_X " 
                << left << "W_Strip_V "<< left << "N_Channel_Phi " 
                << left << "N_Channel_V " << left << "Orientation " << endl;
            cout    << left << setw(12) << m_CgemLayer[i]->getInnerROfCgemLayer()
                << left << setw(12) << m_CgemLayer[i]->getLengthOfCgemLayer()
                << left << setw(8)  << m_CgemLayer[i]->getNumberOfSheet()    
                << left << setw(8)  << m_CgemLayer[i]->getWidthOfSheet()       
                << left << setw(8)  << m_CgemLayer[i]->getAngleOfStereo()     
                << left << setw(10)  << m_CgemLayer[i]->getWidthOfPitchX()
                << left << setw(10)  << m_CgemLayer[i]->getWidthOfPitchV()
                << left << setw(10) << m_CgemLayer[i]->getWidthOfStripX()     
                << left << setw(10) << m_CgemLayer[i]->getWidthOfStripV()     
                << left << setw(14) << m_CgemLayer[i]->getNumberOfChannelX()
                << left << setw(12) << m_CgemLayer[i]->getNumberOfChannelV() 
                << left << setw(11) << m_CgemLayer[i]->getOrientation() << endl;
 
            cout    << left << "T_Cathode L"<<i << left << "T_GemFoil L"<<i << left << "T_Anode L"<<i
                << left << "T_CgemLayer "<<i << endl;
            cout    << left << setw(10) << T_Cathode << left << setw(10) << T_GemFoil
                << left << setw(8)  << T_Anode  << left << setw(9) << m_T_CgemLayer[i] << endl;
 
            cout    << "from r= " << m_CgemLayer[i]->getInnerROfCgemLayer() 
                << " to r="    << m_CgemLayer[i]->getInnerROfCgemLayer() + m_T_CgemLayer[i] << endl;
            cout    << "thickness " <<  m_T_CgemLayer[i] << endl;
        }
 
        /* Information of CgemLayer */
        m_CgemLayer[i]->setIDOfCgemLayer(i);
        m_CgemLayer[i]->setOuterROfCgemLayer(m_CgemLayer[i]->getInnerROfCgemLayer() + m_T_CgemLayer[i]);
 
        /* Information of Cathode */
        bool fReversed = m_CgemLayer[i]->getOrientation();
        if(fReversed == false) {
            m_CgemLayer[i]->setInnerROfCathode( m_CgemLayer[i]->getInnerROfCgemLayer() );
            m_CgemLayer[i]->setOuterROfCathode( m_CgemLayer[i]->getInnerROfCgemLayer() + T_Cathode); 
        } else {
            m_CgemLayer[i]->setOuterROfCathode( m_CgemLayer[i]->getOuterROfCgemLayer() );
            m_CgemLayer[i]->setInnerROfCathode( m_CgemLayer[i]->getOuterROfCgemLayer() - T_Cathode);
        }
 
        if(m_debugCout) cout << "CATHODE from " <<  m_CgemLayer[i]->getInnerROfCathode() << " to " <<  m_CgemLayer[i]->getOuterROfCathode() << " thickness " <<  m_CgemLayer[i]->getOuterROfCathode() -  m_CgemLayer[i]->getInnerROfCathode() << endl;
 
        int counter = 0;
        double inner_r_cathode[NMATERIALS_CATHODE];
        inner_r_cathode[0] =   m_CgemLayer[i]->getInnerROfCathode();
        if(fReversed == true) counter = NMATERIALS_CATHODE-1;
 
        for(int imat = 0; imat < NMATERIALS_CATHODE; imat++) {
 
            switch(counter) {
                case 0 :
                    m_CgemLayer[i]->setInnerROfCathodeCu1(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeCu1(inner_r_cathode[imat]+tc_cu1[i]);
                    inner_r_cathode[imat]+=tc_cu1[i];
                    break;
                case 1 :
                    m_CgemLayer[i]->setInnerROfCathodeKapton1(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeKapton1(inner_r_cathode[imat]+tc_kapton1[i]);
                    inner_r_cathode[imat]+=tc_kapton1[i];
                    break;
                case 2 :
                    m_CgemLayer[i]->setInnerROfCathodeEpoxy0(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeEpoxy0(inner_r_cathode[imat]+tc_epoxy0[i]);
                    inner_r_cathode[imat]+=tc_epoxy0[i];
                    break;
                case 3 :
                    m_CgemLayer[i]->setInnerROfCathodeCarbonf(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeCarbonf(inner_r_cathode[imat]+tc_carbonf[i]);
                    inner_r_cathode[imat]+=tc_carbonf[i];
                    break;
                case 4 :
                    m_CgemLayer[i]->setInnerROfCathodeEpoxy1(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeEpoxy1(inner_r_cathode[imat]+tc_epoxy1[i]);
                    inner_r_cathode[imat]+=tc_epoxy1[i];
                    break;
                case 5 :
                    m_CgemLayer[i]->setInnerROfCathodeHoneycomb(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeHoneycomb(inner_r_cathode[imat]+tc_honeycomb[i]);
                    inner_r_cathode[imat]+=tc_honeycomb[i];
                    break;
                case 6 :
                    m_CgemLayer[i]->setInnerROfCathodeRohacell1(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeRohacell1(inner_r_cathode[imat]+tc_rohacell1[i]);
                    inner_r_cathode[imat]+=tc_rohacell1[i];
                    break;
                case 7 :
                    m_CgemLayer[i]->setInnerROfCathodeEpoxy2(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeEpoxy2(inner_r_cathode[imat]+tc_epoxy2[i]);
                    inner_r_cathode[imat]+=tc_epoxy2[i];
                    break;
                case 8 :
                    m_CgemLayer[i]->setInnerROfCathodeKapton2(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeKapton2(inner_r_cathode[imat]+tc_kapton2[i]);
                    inner_r_cathode[imat]+=tc_kapton2[i];
                    break;
                case 9 :
                    m_CgemLayer[i]->setInnerROfCathodeEpoxy3(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeEpoxy3(inner_r_cathode[imat]+tc_epoxy3[i]);
                    inner_r_cathode[imat]+=tc_epoxy3[i];
                    break;
                case 10 :
                    m_CgemLayer[i]->setInnerROfCathodeRohacell2(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeRohacell2(inner_r_cathode[imat]+tc_rohacell2[i]);
                    inner_r_cathode[imat]+=tc_rohacell2[i];
                    break;
                case 11 :
                    m_CgemLayer[i]->setInnerROfCathodeEpoxy4(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeEpoxy4(inner_r_cathode[imat]+tc_epoxy4[i]);
                    inner_r_cathode[imat]+=tc_epoxy4[i];
                    break;
                case 12 :
                    m_CgemLayer[i]->setInnerROfCathodeKapton3(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeKapton3(inner_r_cathode[imat]+tc_kapton3[i]);
                    inner_r_cathode[imat]+=tc_kapton3[i];
                    break;
                case 13 :
                    m_CgemLayer[i]->setInnerROfCathodeCu2(inner_r_cathode[imat]);
                    m_CgemLayer[i]->setOuterROfCathodeCu2(inner_r_cathode[imat]+tc_cu2[i]);
                    inner_r_cathode[imat]+=tc_cu2[i];
 
                    if(m_debugCout) cout << " CATHODE CU INNER " << inner_r_cathode[imat] << " OUTER " << inner_r_cathode[imat]+tc_cu2[i] << endl;
 
                    break;
                default:
                    break;
            }
 
            //  cout << imat << " " << counter << " inner " << inner_r_cathode[imat] << " outer " << outer_r_cathode[imat] << endl;
            if(imat == NMATERIALS_CATHODE - 1) {
                if(m_debugCout) cout << "breaking" << endl;
                break;
            }
            inner_r_cathode[imat + 1] = inner_r_cathode[imat];
            if(fReversed == false) counter++;
            else counter--;
        }
        if(m_debugCout) 
        {
            cout << "T_cathode: " << setw(15) <<T_Cathode; 
            cout << setw(20) << m_CgemLayer[i]->getInnerROfCathode() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeCu1()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeCu1()
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeKapton1()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeKapton1()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeEpoxy0()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeEpoxy0()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeCarbonf()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeCarbonf()
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeEpoxy1()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeEpoxy1()<<endl
                <<  setw(20) << m_CgemLayer[i]->getInnerROfCathodeHoneycomb()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeHoneycomb()
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeRohacell1()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeRohacell1()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeEpoxy2()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeEpoxy2() 
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeKapton2()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeKapton2()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeEpoxy3()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeEpoxy3() 
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeRohacell2()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeRohacell2()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeEpoxy4()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeEpoxy4() 
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeKapton3()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeKapton3()<<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfCathodeCu2()
                << setw(20) << m_CgemLayer[i]->getOuterROfCathodeCu2() << endl
                << setw(20) << m_CgemLayer[i]->getOuterROfCathode() << endl;
        }
 
 
        /* Check over the radius */
        if(fReversed == false) { 
            if (fabs((m_CgemLayer[i]->getOuterROfCathode() / m_CgemLayer[i]->getOuterROfCathodeCu2()) - 1) > FLT_EPSILON){
                cout << fabs(m_CgemLayer[i]->getOuterROfCathode() - m_CgemLayer[i]->getOuterROfCathodeCu2()) << endl;
                cout << setw(20) << m_CgemLayer[i]->getOuterROfCathode()
                    << setw(20) << m_CgemLayer[i]->getOuterROfCathodeCu2()
                    << ",  FLT_EPSILON = " << FLT_EPSILON << endl;
                cout << "Error : CgemGeomSvc::initGeom(), Cathode size is wrong!" << endl;
                return false;
            }
        }
        else {
            if (fabs((m_CgemLayer[i]->getInnerROfCathode() / m_CgemLayer[i]->getInnerROfCathodeCu1()) - 1) > FLT_EPSILON){//?
                cout << fabs(m_CgemLayer[i]->getInnerROfCathode() - m_CgemLayer[i]->getInnerROfCathodeCu1()) << endl;
                cout << setw(20) << m_CgemLayer[i]->getInnerROfCathode()
                    << setw(20) << m_CgemLayer[i]->getInnerROfCathodeCu1()
                    << ",  FLT_EPSILON = " << FLT_EPSILON << endl;
                cout << "Error : CgemGeomSvc::initGeom(), Cathode size is wrong!" << endl;
                return false;
            }
        }
 
        /* Information of Gap and GemFoil */
        if(fReversed == false) {
            m_CgemLayer[i]->setInnerROfGapD( m_CgemLayer[i]->getOuterROfCathode() );
            m_CgemLayer[i]->setOuterROfGapD( m_CgemLayer[i]->getOuterROfCathode() + gap_d[i] );
            m_CgemLayer[i]->setInnerROfGapT1( m_CgemLayer[i]->getOuterROfGapD() + tg_cu1[i]+tg_kapton[i]+tg_cu2[i] );
            m_CgemLayer[i]->setOuterROfGapT1( m_CgemLayer[i]->getInnerROfGapT1() + tg_gt1[i] );
            m_CgemLayer[i]->setInnerROfGapT2( m_CgemLayer[i]->getOuterROfGapT1() + tg_cu1[i]+tg_kapton[i]+tg_cu2[i]);
            m_CgemLayer[i]->setOuterROfGapT2( m_CgemLayer[i]->getInnerROfGapT2() + tg_gt2[i] );
            m_CgemLayer[i]->setInnerROfGapI( m_CgemLayer[i]->getOuterROfGapT2() + tg_cu1[i]+tg_kapton[i]+tg_cu2[i] );
            m_CgemLayer[i]->setOuterROfGapI( m_CgemLayer[i]->getInnerROfGapI() + tg_gi[i]);
        }
        else {//?
            m_CgemLayer[i]->setOuterROfGapD( m_CgemLayer[i]->getInnerROfCathode() );
            m_CgemLayer[i]->setInnerROfGapD( m_CgemLayer[i]->getInnerROfCathode() - gap_d[i]);
            m_CgemLayer[i]->setOuterROfGapT1( m_CgemLayer[i]->getInnerROfGapD() - tg_cu1[i]-tg_kapton[i]-tg_cu2[i]);
            m_CgemLayer[i]->setInnerROfGapT1( m_CgemLayer[i]->getOuterROfGapT1() - tg_gt1[i]);
            m_CgemLayer[i]->setOuterROfGapT2( m_CgemLayer[i]->getInnerROfGapT1() - tg_cu1[i]-tg_kapton[i]-tg_cu2[i] );
            m_CgemLayer[i]->setInnerROfGapT2( m_CgemLayer[i]->getOuterROfGapT2() - tg_gt2[i]);
            m_CgemLayer[i]->setOuterROfGapI( m_CgemLayer[i]->getInnerROfGapT2() - tg_cu1[i]-tg_kapton[i]-tg_cu2[i]);
            m_CgemLayer[i]->setInnerROfGapI( m_CgemLayer[i]->getOuterROfGapI() - tg_gi[i]);
        }
        
        if(m_debugCout) 
        {
            cout << "Geom foil&gap " << setw(15) << gap_d[i] << setw(15)
                << tg_cu1[i]+tg_kapton[i]+tg_cu2[i] << setw(15) << tg_gt1[i]
                << setw(15) << tg_gt2[i] << setw(15) << tg_gi[i] << endl;
 
            cout << "gaps: -----------------------------" << endl;
            cout << setw(20) << m_CgemLayer[i]->getInnerROfGapD()
                << setw(20) << m_CgemLayer[i]->getOuterROfGapD()
                << setw(20) << m_CgemLayer[i]->getInnerROfGapT1()
                << setw(20) << m_CgemLayer[i]->getOuterROfGapT1() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfGapT2()
                << setw(20) << m_CgemLayer[i]->getOuterROfGapT2()
                << setw(20) << m_CgemLayer[i]->getInnerROfGapI()
                << setw(20) << m_CgemLayer[i]->getOuterROfGapI() << endl;
        }
 
        /* The number of GemFoil can only be 3 in the code! */
        if (m_N_GemFoil != 3){
            cout << "Error : CgemGeomSvc::initGeom(), GemFoil number in code is 3! If m_N_GemFoil != 3,you should modify the code,especially the gap information!" << endl;
            return false;
        }
 
        double R_i_GemFoil[3] = {0};
        if(fReversed == false) {
            R_i_GemFoil[0] = m_CgemLayer[i]->getOuterROfGapD();
            R_i_GemFoil[1] = m_CgemLayer[i]->getOuterROfGapT1();
            R_i_GemFoil[2] = m_CgemLayer[i]->getOuterROfGapT2();
        }
        else {
            R_i_GemFoil[0] = m_CgemLayer[i]->getOuterROfGapT1();
            R_i_GemFoil[1] = m_CgemLayer[i]->getOuterROfGapT2();
            R_i_GemFoil[2] = m_CgemLayer[i]->getOuterROfGapI();           
        }
 
        for(int j=0; j < m_N_GemFoil; j++){
            m_foil[i][j] = new CgemGeoFoil();
 
            m_foil[i][j]->setIDOfCgemFoil(j);
 
 
            double inner_r_foil[N_MATERIAL_MAX];
            double outer_r_foil[N_MATERIAL_MAX];
            double m_t_foil_list[N_MATERIAL_MAX] = {tg_cu1[i], tg_kapton[i], tg_cu2[i]};
 
            inner_r_foil[0] = R_i_GemFoil[j];
 
            if(fReversed == true) counter = ngMaterials[i] - 1;
            else counter = 0;
 
            m_foil[i][j]->setInnerROfCgemFoilHole(r_i_hole[j]);
            m_foil[i][j]->setOuterR1OfCgemFoilHole(r_o1_hole[j]);
            m_foil[i][j]->setOuterR2OfCgemFoilHole(r_o2_hole[j]);
            m_foil[i][j]->setLengthOfCgemFoilHole(l_hole[j]);
 
            for(int imat = 0; imat < ngMaterials[i]; imat++) {
                outer_r_foil[imat] = inner_r_foil[imat] + m_t_foil_list[counter];
 
                switch(counter) {
                    case 0 :
                        m_foil[i][j]->setInnerROfCgemFoilCu1(inner_r_foil[imat]);
                        m_foil[i][j]->setOuterROfCgemFoilCu1(outer_r_foil[imat]);
                        break;
                    case 1 :
                        m_foil[i][j]->setInnerROfCgemFoilKapton(inner_r_foil[imat]);
                        m_foil[i][j]->setOuterROfCgemFoilKapton(outer_r_foil[imat]);
                        break;
                    case 2 :
                        m_foil[i][j]->setInnerROfCgemFoilCu2(inner_r_foil[imat]);
                        m_foil[i][j]->setOuterROfCgemFoilCu2(outer_r_foil[imat]);
                        break;
                    default:
                        break;
                }
 
                // cout << imat << " " << counter << " inner " << inner_r_foil[imat] << " outer " << outer_r_foil[imat] << endl;
                if(imat == ngMaterials[i] - 1) {
                    if(m_debugCout) cout << "breaking" << endl;
                    break;
                }
                inner_r_foil[imat + 1] = outer_r_foil[imat];
                if(fReversed == false) counter++;
                else counter--;
            }
 
            m_foil[i][j]->setInnerROfCgemFoil( R_i_GemFoil[j] );
            m_foil[i][j]->setOuterROfCgemFoil( R_i_GemFoil[j] + tg_cu1[i]+tg_kapton[i]+tg_cu2[i] );
 
            m_CgemLayer[i]->setCgemFoil(j, m_foil[i][j]);
            
            if(m_debugCout) 
            {
                cout << "Geom foil:  layer " << i << "  " << j << endl;
                cout << setw(15) << tg_cu1[i] << setw(15)
                    << tg_kapton[i] << setw(15) << tg_cu2[i] << endl;
 
                cout << setw(20) << m_foil[i][j]->getIDOfCgemFoil()
                    << setw(20) << m_foil[i][j]->getInnerROfCgemFoil()
                    << setw(20) << m_foil[i][j]->getInnerROfCgemFoilCu1()
                    << setw(20) << m_foil[i][j]->getOuterROfCgemFoilCu1()
                    << setw(20) << m_foil[i][j]->getInnerROfCgemFoilKapton() << endl
                    << setw(20) << m_foil[i][j]->getOuterROfCgemFoilKapton()
                    << setw(20) << m_foil[i][j]->getInnerROfCgemFoilCu2()
                    << setw(20) << m_foil[i][j]->getOuterROfCgemFoilCu2()
                    << setw(20) << m_foil[i][j]->getOuterROfCgemFoil() << endl;
            }
        }
        
        if(m_debugCout) 
        {
            if(fReversed) {
                cout << "GEM" << endl;
                cout << "I " <<  m_CgemLayer[i]->getInnerROfGapI() << " " <<  m_CgemLayer[i]->getOuterROfGapI() << endl;      
                cout << m_foil[i][2]->getInnerROfCgemFoil() << " " <<  m_foil[i][2]->getOuterROfCgemFoil() << endl;
                cout << "T2 " <<  m_CgemLayer[i]->getInnerROfGapT2() << " " <<  m_CgemLayer[i]->getOuterROfGapT2() << endl;      
                cout << m_foil[i][1]->getInnerROfCgemFoil() << " " <<  m_foil[i][1]->getOuterROfCgemFoil() << endl;
                cout << "T1 " <<  m_CgemLayer[i]->getInnerROfGapT1() << " " <<  m_CgemLayer[i]->getOuterROfGapT1() << endl;      
                cout << m_foil[i][0]->getInnerROfCgemFoil() << " " <<  m_foil[i][0]->getOuterROfCgemFoil() << endl;
                cout << "D " <<  m_CgemLayer[i]->getInnerROfGapD() << " " <<  m_CgemLayer[i]->getOuterROfGapD() << endl;      
            }
            else {
                cout << "GEM" << endl;
                cout << "D " <<  m_CgemLayer[i]->getInnerROfGapD() << " " <<  m_CgemLayer[i]->getOuterROfGapD() << endl;      
                cout << m_foil[i][0]->getInnerROfCgemFoil() << " " <<  m_foil[i][0]->getOuterROfCgemFoil() << endl;
                cout << "T1 " <<  m_CgemLayer[i]->getInnerROfGapT1() << " " <<  m_CgemLayer[i]->getOuterROfGapT1() << endl;      
                cout << m_foil[i][1]->getInnerROfCgemFoil() << " " <<  m_foil[i][1]->getOuterROfCgemFoil() << endl;
                cout << "T2 " <<  m_CgemLayer[i]->getInnerROfGapT2() << " " <<  m_CgemLayer[i]->getOuterROfGapT2() << endl;      
                cout << m_foil[i][2]->getInnerROfCgemFoil() << " " <<  m_foil[i][2]->getOuterROfCgemFoil() << endl;
                cout << "I " <<  m_CgemLayer[i]->getInnerROfGapI() << " " <<  m_CgemLayer[i]->getOuterROfGapI() << endl;      
            }
        }
 
        /* Information of Anode */
        if(fReversed == false) {
            m_CgemLayer[i]->setInnerROfAnode( m_CgemLayer[i]->getOuterROfGapI() );
            m_CgemLayer[i]->setOuterROfAnode( m_CgemLayer[i]->getOuterROfGapI() + T_Anode );
        } else {
            m_CgemLayer[i]->setOuterROfAnode( m_CgemLayer[i]->getInnerROfGapI() );
            m_CgemLayer[i]->setInnerROfAnode( m_CgemLayer[i]->getInnerROfGapI() - T_Anode);
        }
 
        if(m_debugCout) cout << "ANODE from " <<  m_CgemLayer[i]->getInnerROfAnode() << " to " <<  m_CgemLayer[i]->getOuterROfAnode() << " thickness " <<  m_CgemLayer[i]->getOuterROfAnode() -  m_CgemLayer[i]->getInnerROfAnode() << endl;
 
 
 
        double inner_r_anode[N_MATERIAL_MAX];
        double outer_r_anode[N_MATERIAL_MAX];
        double m_t_anode_list[N_MATERIAL_MAX] = {
            ta_cu1[i],ta_kapton1[i],ta_cu2[i],ta_epoxy1[i],ta_kapton2[i],ta_epoxy2[i],ta_rohacell1[i],ta_carbonf1[i],ta_kapton3[i],ta_epoxy3[i],ta_honeycomb[i],ta_kapton4[i],ta_epoxy4[i],ta_rohacell2[i],ta_carbonf2[i],ta_epoxy5[i],ta_cu3[i],ta_kapton5[i],ta_kapton6[i],ta_cu4[i]
        };
        inner_r_anode[0] = m_CgemLayer[i]->getInnerROfAnode();
 
        if(fReversed == true) counter = NMATERIALS_ANODE - 1;
        else counter = 0;
 
        for(int imat = 0; imat < NMATERIALS_ANODE; imat++) {
            outer_r_anode[imat] = inner_r_anode[imat] + m_t_anode_list[counter];
 
            switch(counter) {
                case 0 :
                    m_CgemLayer[i]->setInnerROfAnodeCu1(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCu1(outer_r_anode[imat]);
                    break;
                case 1 :
                    m_CgemLayer[i]->setInnerROfAnodeKapton1(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton1(outer_r_anode[imat]);
                    break;
                case 2:
                    m_CgemLayer[i]->setInnerROfAnodeCu2(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCu2(outer_r_anode[imat]);
                    break;
                case 3:
                    m_CgemLayer[i]->setInnerROfAnodeEpoxy1(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeEpoxy1(outer_r_anode[imat]);
                    break;
                case 4:
                    m_CgemLayer[i]->setInnerROfAnodeKapton2(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton2(outer_r_anode[imat]);
                    break;
                case 5:
                    m_CgemLayer[i]->setInnerROfAnodeEpoxy2(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeEpoxy2(outer_r_anode[imat]);
                    break;
                case 6:
                    m_CgemLayer[i]->setInnerROfAnodeRohacell1(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeRohacell1(outer_r_anode[imat]);
                    break;
                case 7:
                    m_CgemLayer[i]->setInnerROfAnodeCarbonf1(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCarbonf1(outer_r_anode[imat]);
                    break;
                case 8 :
                    m_CgemLayer[i]->setInnerROfAnodeKapton3(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton3(outer_r_anode[imat]);
                    break;
                case 9:
                    m_CgemLayer[i]->setInnerROfAnodeEpoxy3(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeEpoxy3(outer_r_anode[imat]);
                    break;
                case 10:
                    m_CgemLayer[i]->setInnerROfAnodeHoneycomb(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeHoneycomb(outer_r_anode[imat]);
                    break;
                case 11:
                    m_CgemLayer[i]->setInnerROfAnodeKapton4(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton4(outer_r_anode[imat]);
                    break;
                case 12:
                    m_CgemLayer[i]->setInnerROfAnodeEpoxy4(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeEpoxy4(outer_r_anode[imat]);
                    break;
                case 13:
                    m_CgemLayer[i]->setInnerROfAnodeRohacell2(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeRohacell2(outer_r_anode[imat]);
                    break;
                case 14:
                    m_CgemLayer[i]->setInnerROfAnodeCarbonf2(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCarbonf2(outer_r_anode[imat]);
                    break;
                case 15:
                    m_CgemLayer[i]->setInnerROfAnodeEpoxy5(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeEpoxy5(outer_r_anode[imat]);
                    break;
                case 16:
                    m_CgemLayer[i]->setInnerROfAnodeCu3(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCu3(outer_r_anode[imat]);
                    break;
                case 17:
                    m_CgemLayer[i]->setInnerROfAnodeKapton5(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton5(outer_r_anode[imat]);
                    break;
                case 18:
                    m_CgemLayer[i]->setInnerROfAnodeKapton6(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeKapton6(outer_r_anode[imat]);
                    break;
                case 19:
                    m_CgemLayer[i]->setInnerROfAnodeCu4(inner_r_anode[imat]);
                    m_CgemLayer[i]->setOuterROfAnodeCu4(outer_r_anode[imat]);
                    break;
                default:
                    break;
            }
 
            //  cout << imat << " " << counter << " inner " << inner_r_anode[imat] << " outer " << outer_r_anode[imat] << endl;
            if(imat == NMATERIALS_ANODE - 1) {
                if(m_debugCout) cout << "breaking" << endl;
                break;
            }
            inner_r_anode[imat + 1] = outer_r_anode[imat];
            if(fReversed == false) counter++;
            else counter--;
        }
        if(m_debugCout) 
        {
            cout << "ANODE: " << endl;
            cout << setw(20) << m_CgemLayer[i]->getInnerROfAnode() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCu1()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCu1()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton1()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton1() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCu2()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCu2()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeEpoxy1()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeEpoxy1() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton2()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton2()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeEpoxy2()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeEpoxy2() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeRohacell1()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeRohacell1()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCarbonf1()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCarbonf1() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton3()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton3() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeEpoxy3()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeEpoxy3()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeHoneycomb()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeHoneycomb() <<endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton4()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton4()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeEpoxy4()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeEpoxy4() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeRohacell2()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeRohacell2()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCarbonf2()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCarbonf2() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeEpoxy5()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeEpoxy5() 
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCu3()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCu3() << endl
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton5()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton5() << endl 
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeKapton6()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeKapton6()
                << setw(20) << m_CgemLayer[i]->getInnerROfAnodeCu4()
                << setw(20) << m_CgemLayer[i]->getOuterROfAnodeCu4() << endl
                << setw(20) << m_CgemLayer[i]->getOuterROfAnode() << endl;
            cout << "X strip description is " << m_B_strip_x_on << " and V strip description is " << m_B_strip_v_on << endl;
        }
 
        /* Check over the radius! */
        if(fReversed == false) { 
            if(fabs((m_CgemLayer[i]->getOuterROfAnode() / m_CgemLayer[i]->getOuterROfCgemLayer())-1) > FLT_EPSILON){
                cout << fabs(m_CgemLayer[i]->getOuterROfAnode() - m_CgemLayer[i]->getOuterROfCgemLayer()) << endl;
                cout << FLT_EPSILON << endl;
                cout << "Error : CgemLayer size is wrong!" << endl;
                return false;
            }
        }
        else {
            if(fabs((m_CgemLayer[i]->getInnerROfAnode() / m_CgemLayer[i]->getInnerROfCgemLayer())-1) > FLT_EPSILON){
                cout << fabs(m_CgemLayer[i]->getInnerROfAnode() - m_CgemLayer[i]->getInnerROfCgemLayer()) << endl;
                cout << FLT_EPSILON << endl;
                cout << "Error : CgemLayer size is wrong!" << endl;
                return false;
            }
 
        }
 
        m_CgemLayer[i]->print();
 
        int nSheet=m_CgemLayer[i]->getNumberOfSheet();
        double rVstrips=m_CgemLayer[i]->getInnerROfAnode();
        double rXstrips = m_CgemLayer[i]->getInnerROfAnodeCu2();
        //double rMidAnode = 0.5*(rVstrips+m_CgemLayer[i]->getOuterROfAnodeCu2());
        if(fReversed) {
            rVstrips=m_CgemLayer[i]->getOuterROfAnode();
            rXstrips = m_CgemLayer[i]->getOuterROfAnodeCu2();
        }
        double width_sheet=rXstrips*CLHEP::twopi/nSheet;
        if(m_debugCout) cout<<"construct CgemReadoutPlanes: "<<endl;
 
        //cout << "!!!!!!!! CHECK DELETE " << rXstrips << " " << width_sheet << " " << endl;
 
        for(int j=0; j<nSheet; j++){
            if(m_debugCout) cout<<"------ layer "<<i<<", sheet "<<j<<" ------ "<<endl;
            double rStrips = 0.5*(rXstrips+rVstrips);
            dX_start_strip[i][j]=(2.*M_PI*rStrips-m_CgemLayer[i]->getNumberOfChannelX()*m_CgemLayer[i]->getWidthOfPitchX())/(2.*nSheet)*rXstrips/rStrips;// dX calculation FIXME
            m_ReadoutPlane[i][j]=new CgemGeoReadoutPlane(i,j,
                    rXstrips, rVstrips, 
                    phi_start[i][j], dX_start_strip[i][j], 
                    dV_start_strip[i][j],
                    width_sheet,
                    -0.5*m_CgemLayer[i]->getLengthOfCgemLayer(),
                    m_CgemLayer[i]->getLengthOfCgemLayer(),
                    m_CgemLayer[i]->getNumberOfChannelX()/nSheet,
                    m_CgemLayer[i]->getNumberOfChannelV()/nSheet,
                    m_CgemLayer[i]->getWidthOfPitchX(),
                    m_CgemLayer[i]->getWidthOfStripX(),
                    m_CgemLayer[i]->getWidthOfPitchV(),
                    m_CgemLayer[i]->getWidthOfStripV(),
                    m_CgemLayer[i]->getAngleOfStereo(),
                    m_CgemLayer[i]->getMiddleROfGapD(), 
                    m_CgemLayer[i]->getOuterROfGapD() 
                    );
        }// loop sheet
 
    } /* end of "for (int i=0; i < m_N_CgemLayer; i++)" */
 
    if(m_debugCout) cout << "INFO : CgemGeomSvc::initGeom(), Successfully input geometry parameters!" << endl;
 
    /** CHECK: now it is the passive dimensions
      m_R_i_Cgem = m_CgemLayer[0]->getInnerROfCgemLayer();
      m_R_o_Cgem = m_CgemSeparator->getOuterROfSeparator(); // m_CgemLayer[m_N_CgemLayer-1]->getOuterROfCgemLayer(); // CHECK
      m_L_Cgem   = m_CgemSeparator->getLengthOfSeparator(); // m_CgemLayer[m_N_CgemLayer-1]->getLengthOfCgemLayer(); // CHECK
     **/
    //if(m_debugCout) 
    {
        cout << "====================CGEM Information=========================" << endl;
        cout << "R_i_Cgem  " << "R_o_Cgem  " << "L_Cgem" << endl;
        cout << left << setw(10) << m_R_i_Cgem  
            << left << setw(10) << m_R_o_Cgem  
            << left << setw(10) << m_L_Cgem << endl;
        cout << "====================CGEM Information End!====================" << endl;
        cout << " " << endl;
 
 
        cout << "************************ FINAL CHECK OF THE DIMENSIONS" << endl;
        for(int ilayer = 0; ilayer < 3; ilayer++) {
 
 
            cout << "inner global radius of layer" << ilayer << ": " << m_CgemLayer[ilayer]->getInnerROfCgemLayer() 
                << " outer "                                        << m_CgemLayer[ilayer]->getOuterROfCgemLayer() << endl; 
            cout << "inner  radius of cathode: "                     << m_CgemLayer[ilayer]->getInnerROfCathode()
                << " outer "                                        << m_CgemLayer[ilayer]->getOuterROfCathode() << endl;
 
            for(int ifoil = 0; ifoil < 3; ifoil++) {
                cout << "inner  radius of gem" << ifoil << ": "        << m_foil[ilayer][ifoil]->getInnerROfCgemFoil()
                    << " outer "                                      << m_foil[ilayer][ifoil]->getOuterROfCgemFoil() << endl;
            }
 
            cout << "inner  radius of anode: "                     << m_CgemLayer[ilayer]->getInnerROfAnode()
                << " outer "                                      << m_CgemLayer[ilayer]->getOuterROfAnode() << endl;
        }
        cout << endl;
    }
 
    return true;
} /* end of "Execute()" */
 
 
void CgemGeomSvc::TestReadoutPlane() {
    const int nSheet[3]={1,2,2};
    //const double phimin(-3.14159265358979312), phimax(3.14159265358979312);
    //const double phimin(-3.14159265358979312), phimax(-3.);
    //const double phimin(-0.1), phimax(0.1);
    const double phimin(0.), phimax(7.);
    const int nphi=20000;
    const double zmin(-850), zmax(850);// L: 532, 690, 847,  0.5*L: 266, 345, 423.5
    const int nz=100000;
 
    double phi(0.), z(0), v(0.);
    //double phi(0.), z(-423.5), v(0.);
 
    int onPlane[3][2];
    int Xid[3][2];
    int Xid2[3][2];
    int Xid3[3][2];
    int Xid4[3][2];
    int Vid[3][2];
    int Vid2[3][2];
    int Vid3[3][2];
    int Vid4[3][2];
    double X[3][2], dX2[3][2], dX3[3][2], PHI[3][2],
           V[3][2], dV2[3][2], dV3[3][2],
           zCross[3][2];
 
    cout.precision(10);
    cout<<setiosflags(ios::left)
        <<setw(14)<<"phi"      <<setw(2)<<"  " 
        <<setw(14)<<"z"        <<setw(2)<<"  "
        <<setw( 6)<<"layer"    <<setw(2)<<"  "
        <<setw( 6)<<"sheet"    <<setw(2)<<"  "
        <<setw( 6)<<"onIt"     <<setw(2)<<"  "
        <<setw(14)<<"X"        <<setw(2)<<"  "
        <<setw( 6)<<"Xid"      <<setw(2)<<"  "
        <<setw( 6)<<"Xid2"     <<setw(2)<<"  "
        <<setw( 6)<<"Xid3"     <<setw(2)<<"  "
        <<setw(18)<<"dX2"      <<setw(2)<<"  "
        <<setw(18)<<"dX3"      <<setw(2)<<"  "
        <<setw(18)<<"PHI"      <<setw(2)<<"  "
        <<setw(18)<<"V"        <<setw(2)<<"  "
        <<setw( 6)<<"Vid"      <<setw(2)<<"  "
        <<setw( 6)<<"Vid2"     <<setw(2)<<"  "
        <<setw( 6)<<"Vid3"     <<setw(2)<<"  "
        <<setw(18)<<"dV2"      <<setw(2)<<"  "
        <<setw(18)<<"dV3"      <<setw(2)<<"  "
        <<setw(18)<<"zCross"   <<setw(2)<<"  "
        <<setw(18)<<"Vnext"    <<setw(2)<<"  "
        <<setw(18)<<"VStripLen"<<setw(2)<<"  "
        <<endl;
 
    //for(int ii=0; ii<nz; ii++) 
    for(int ii=0; ii<nphi; ii++) 
    {
        phi= phimin+(ii+0.5)*(phimax-phimin)/nphi;
        //z  = zmin+(ii+0.5)*(zmax-zmin)/nz;
        //cout<<setiosflags(ios::left)
        for(int i=0; i<m_N_CgemLayer; i++)
        {
            //cout<<"i = "<<i<<", nsheet = "<<nSheet[i]<<endl;
            for(int j=0; j<nSheet[i]; j++)
            {
                //cout<<"j = "<<j<<endl;
                CgemGeoReadoutPlane* anode = getReadoutPlane(i,j);
                onPlane[i][j] =anode->OnThePlane(phi,z);
                X[i][j]       =anode->getX(phi);
                Xid[i][j]     =anode->getXStripID(phi);
                Xid2[i][j]    =anode->getClosestXStripID(phi,dX2[i][j]);
                dX3[i][j]     =anode->getDist2ClosestXStripCenter(phi,Xid3[i][j]);
                V[i][j]       =anode->getVFromPhiZ(phi,z);
                Vid[i][j]     =anode->getVIDFromV(V[i][j]);
                G4ThreeVector pos(cos(phi),sin(phi),z);
                Vid2[i][j]    =anode->getClosestVStripID(pos, dV2[i][j]);
                dV3[i][j]     =anode->getDist2ClosestVStripCenter(pos, Vid3[i][j]);
                anode->getStripID(pos,Xid4[i][j],Vid4[i][j]);
                PHI[i][j]     =anode->getPhiFromXID(Xid[i][j]);
                zCross[i][j]  =anode->getZFromPhiV(phi,V[i][j]);
                int iSheetNext= (j+1)%nSheet[i];
                //cout<<"isheet, inext = "<<j<<", "<<iSheetNext<<endl;
                CgemGeoReadoutPlane* anodeNext = getReadoutPlane(i,iSheetNext);
                double phimin_next = anodeNext->getPhimin();
                double Vnext       = anode->getVInNextSheetFromV(V[i][j],phimin_next);
                double VStripLen   = anode->getVStripLength(Vid[i][j]);
                // getVIDInNextSheetFromVID
                //cout<<phi<<" "<<z<<" "<<i<<" "<<j<<" "<<int(onPlane[i][j])<<" "<<V[i][j]<<" "<<zCross[i][j]<<endl;
                cout<<setiosflags(ios::left)
                    <<setw(14)<<phi                      <<setw(2)<<"  "
                    <<setw(14)<<z                        <<setw(2)<<"  "
                    <<setw( 6)<<i                        <<setw(2)<<"  "
                    <<setw( 6)<<j                        <<setw(2)<<"  "
                    <<setw( 6)<<int(onPlane[i][j])       <<setw(2)<<"  "
                    <<setw(14)<<X[i][j]                  <<setw(2)<<"  "
                    <<setw( 6)<<Xid[i][j]                <<setw(2)<<"  "
                    <<setw( 6)<<Xid2[i][j]               <<setw(2)<<"  "
                    <<setw( 6)<<Xid3[i][j]               <<setw(2)<<"  "
                    <<setw(18)<<dX2[i][j]                <<setw(2)<<"  "
                    <<setw(18)<<dX3[i][j]                <<setw(2)<<"  "
                    //<<setw(18)<<fabs(dX3[i][j])-0.29
                    <<setw(18)<<PHI[i][j]                <<setw(2)<<"  "
                    <<setw(18)<<V[i][j]                  <<setw(2)<<"  "
                    <<setw( 6)<<Vid[i][j]                <<setw(2)<<"  "
                    <<setw( 6)<<Vid2[i][j]               <<setw(2)<<"  "
                    <<setw( 6)<<Vid3[i][j]               <<setw(2)<<"  "
                    <<setw(18)<<dV2[i][j]                <<setw(2)<<"  "
                    <<setw(18)<<dV3[i][j]                <<setw(2)<<"  "
                    <<setw(18)<<zCross[i][j]             <<setw(2)<<"  "
                    <<setw(18)<<Vnext                    <<setw(2)<<"  "
                    <<setw(18)<<VStripLen                <<setw(2)<<"  "
                    <<endl
                    ;
            }// sheet
        }// layer
        //cout<<endl;
    }// loop phi
}