void rad_length() {
    
    gROOT->LoadMacro("$VMCWORKDIR/example/gconfig/basiclibs.C");
    basiclibs();
    
    gSystem->Load("libFairTools");
    gSystem->Load("libFairDB");
    gSystem->Load("libGeoBase");
    gSystem->Load("libParBase");
    gSystem->Load("libBase");
    gSystem->Load("libMCStack");
    gSystem->Load("libGen");
    gSystem->Load("libfield");
    
    gSystem->Load("libEICPASSIVE");
    gSystem->Load("libFairSimpletracker");
    gSystem->Load("libFairDircDet");
    gSystem->Load("libFairFsceDet");
    gSystem->Load("libFairEmlfDet");
    gSystem->Load("libFairFgtDet");
    gSystem->Load("libFairHcalDet");
    gSystem->Load("libFairHtckDet");
    gSystem->Load("libFairSitDet");
    gSystem->Load("libFairTpcDet");
    gSystem->Load("libFairTrsDet");
    
    TFile* f = new TFile("$VMCWORKDIR/example/eicroot/simpletracker/macros/data/eic_rad_length.mc.root");
    f->Get("FairBaseParSet");
    
    Int_t nbound = 0;
    Float_t length = 0;
    Float_t safe   = 0;
    Float_t rad    = 0;
    
    vector<string> subs;
    //    subs.push_back("/cave_1/pipe");
    //subs.push_back("/CAVE/SIT");
    //subs.push_back("/CAVE/TPC");
    //subs.push_back("/CAVE/TRS");
    //    subs.push_back("/cave_1/dirc");
    
    vector<string> subn;
    //    subn.push_back("PIPE");
    //subn.push_back("SIT");
    //subn.push_back("TPC");
    //subn.push_back("TRS");
    //    subn.push_back("DIRC");
    
    const int NSUB = subs.size();
    
    vector<void*> h;
    for (int i = 0; i < NSUB; i++) {
        string name = "h[";
        name.append(toString(i));
        name.append("]");
        cout << "name: " << name << endl;
        TH1F* tmp = new TH1F(name.c_str(), subs[i].c_str(), 180, 0, 180);
        tmp->SetLineWidth(1.0);
        tmp->SetLineColor(1);
        tmp->SetMarkerStyle(21);
        tmp->SetMarkerColor(i+2);
        tmp->SetMarkerSize(1.5);
        tmp->SetFillColor(i+2);
        h.push_back(tmp);
    }
    
    int ns = NSUB;
    int i = 1;
    Double_t theta = 0;
    Double_t eta = 0;
    
    for (Int_t i = 0; i <= 180; i++) {
        std::cout << "theta: " << i << "\n";
        for (int j = 0; j < NSUB; j++) {
            ((TH1F*)h[j])->SetBinContent(i, 0);
        }
        //eta = -3.0 + i * 0.025;
        //theta = atan( exp( (-eta) ) ) * 2;
        theta = Float_t(i)*TMath::Pi()/180.0;
        FindRad(subs, NSUB, h, i, 0.0, 0.0, 0.1, theta, 0.0, nbound, length, safe, rad, kFALSE); // xy scan
    }
    
    THStack *hs = new THStack("hs","EIC Detector Geometry: Radiation Length Scan");
    for (int i = 0; i < NSUB; i++) {
        hs->Add( (TH1F*)h[i] );
    }
    
    TCanvas* c1 = new TCanvas("c1","c1", 800, 600);
    c1->SetFrameBorderSize(0);
    hs->Draw();
    
    gStyle->SetLabelSize(0.025, "Y"); 
    
    hs->GetYaxis()->SetTitle("radiation length, X / X0");
    hs->GetXaxis()->SetTitle("Theta");
    
    TLegend *legend=new TLegend(0.7,0.7,0.88,0.88);
    legend->SetTextFont(72);
    legend->SetTextSize(0.04);
    for (int i = 0; i < NSUB; i++) {
        legend->AddEntry( (TH1F*)h[ NSUB - i - 1 ], subn[NSUB - i - 1].c_str(),"p");
    }
    legend->Draw();
    
    c1->SaveAs("VMCWORKDIR/example/eicroot/images/radlen.gif");
    
}


void FindRad(vector<string>& subs, int NSUB, vector<void*>& h, int bin, Double_t x, Double_t y, Double_t z,Double_t theta, Double_t phi, 
             Int_t &nbound, Float_t &length, Float_t &safe, Float_t &rad, Bool_t verbose) {
    
    Double_t xp  = TMath::Sin(theta)*TMath::Cos(phi);
    Double_t yp  = TMath::Sin(theta)*TMath::Sin(phi);
    Double_t zp  = TMath::Cos(theta);
    Double_t snext;
    TString  path;
    Double_t pt[3];
    Double_t loc[3];
    Double_t epsil = 1.E-2;
    Double_t lastrad = 0.;
    Int_t    ismall = 0;
    nbound = 0;
    length = 0.;
    safe   = 0.;
    rad    = 0.;
    TGeoMedium *med;
    TGeoShape  *shape;
    TGeoNode   *lastnode;
    TGeoManager *gGeoManager;
    gGeoManager->InitTrack(x,y,z,xp,yp,zp);
    if ( verbose ) {
        std::cout << "Track:  ["<< x << ", " << y << ", " << z << "], [" << xp << ", " << yp << ", " << zp << "]" << std::endl;
    }
    path = gGeoManager->GetPath();
    
    TGeoNode *nextnode = gGeoManager->GetCurrentNode();
    safe = gGeoManager->Safety();
    while (nextnode) {
        med = 0;
        if (nextnode) {
            med = nextnode->GetVolume()->GetMedium();
        } else {
            return;
        }
        shape = nextnode->GetVolume()->GetShape();
        lastnode = nextnode;
        nextnode = gGeoManager->FindNextBoundaryAndStep();
        snext  = gGeoManager->GetStep();
        if ( snext < 1.e-8 ) {
            ismall++;
            if ( ( ismall < 3 ) && ( lastnode != nextnode ) ) {
                // First try to cross a very thin layer
                length += snext;
                nextnode = gGeoManager->FindNextBoundaryAndStep();
                snext = gGeoManager->GetStep();
                if ( snext < 1.E-8 ) { 
                    continue;
                }
                // We managed to cross the layer
                ismall = 0;
            } else {
                // Relocate point
                if ( ismall > 3 ) {
                    fprintf(stderr,"ERROR: Small steps in: %s shape=%s\n",gGeoManager->GetPath(), shape->ClassName());
                    return;
                }
                memcpy(pt,gGeoManager->GetCurrentPoint(),3*sizeof(Double_t));
                const Double_t *dir = gGeoManager->GetCurrentDirection();
                for( Int_t i = 0; i < 3; i++ ) { 
                    pt[i] += epsil*dir[i]; 
                }
                snext = epsil;
                length += snext;
                rad += lastrad*snext;
                gGeoManager->CdTop();
                nextnode = gGeoManager->FindNode( pt[0], pt[1], pt[2] );
                if ( gGeoManager->IsOutside() ) { 
                    return;
                }
                TGeoMatrix *mat = gGeoManager->GetCurrentMatrix();
                mat->MasterToLocal(pt,loc);
                if ( !gGeoManager->GetCurrentVolume()->Contains(loc) ) {
                    gGeoManager->CdUp();
                    nextnode = gGeoManager->GetCurrentNode();
                }
                continue;
            }
        } else {
            ismall = 0;
        }
        nbound++;
        length += snext;
        if ( med ) {
            Double_t radlen = med->GetMaterial()->GetRadLen();
            if ( radlen > 1.e-5 && radlen < 1.e10 ) {
                lastrad = 1.0 / radlen;
                rad += (lastrad*snext);
            } else {
                lastrad = 0.;
            }
            
            for (int iii = 0; iii < NSUB; iii++) {
                if ( path.BeginsWith( subs[iii].c_str() ) ) {
                    Double_t bcon = ((TH1F*)h[iii])->GetBinContent(bin);
                    ((TH1F*)h[iii])->SetBinContent(bin, bcon + (snext/med->GetMaterial()->GetRadLen()));
                    break;
                }
            }
            
            if ( verbose ) {
                std::cout << " STEP #" << nbound << ": " << path.Data() << " || " << "    step=" << snext 
                <<"  length=" << length << "  radl_x0s = " << (snext/med->GetMaterial()->GetRadLen()) << " " << med->GetName() << std::endl;
            }
            path =  gGeoManager->GetPath();
        }
    }
}

template <class T>
std::string toString(T& arg)
{
    std::ostringstream os;
    os << arg;
    return os.str();
}