AgataGeometryTransformer.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christian Finck        *
 *   stezow(AT)ipnl.in2p3.fr, cfinck(AT)ires.in2p3.fr                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <Riostream.h>
#include <map>

#ifndef Gw_AgataGeometryTransformer
#include "AgataGeometryTransformer.h"
#endif

#ifndef ROOT_TGeoBBox
#include "TGeoBBox.h"
#endif

#ifndef ROOT_TColor
#include "TColor.h"
#endif

#ifndef ROOT_TGeoCone
#include "TGeoCone.h"
#endif 

#ifndef ROOT_TGeoManager
#include "TGeoManager.h"
#endif 

#ifndef ROOT_TGeoPcon
#include "TGeoPcon.h"
#endif

#ifndef ROOT_TGeoShape
#include "TGeoShape.h"
#endif

#ifndef ROOT_TGeoTube
#include "TGeoTube.h"
#endif

#ifndef ROOT_TGeoXtru
#include "TGeoXtru.h"
#endif

#ifndef ROOT_TGeoMatrix
#include "TGeoMatrix.h"
#endif

#ifndef ROOT_TList
#include "TList.h"
#endif

#ifndef ROOT_TMath
#include "TMath.h"
#endif

#ifndef ROOT_TObjArray
#include "TObjArray.h"
#endif

#ifndef ROOT_TVector3
#include "TVector3.h"
#endif

#ifndef ROOT_TGeoCompositeShape
#include "TGeoCompositeShape.h"
#endif


using namespace Gw;
using namespace std;

const Int_t   AgataGeometryTransformer::fgkMaxDetectorID = 180;
const Int_t   AgataGeometryTransformer::fgkMaxDepthID    = 6;
const Float_t AgataGeometryTransformer::fgkDepth[6]      = {0, 8, 21, 36, 54, 72}; 
const Char_t* AgataGeometryTransformer::fgkGeomName      = "Agata";
const Char_t* AgataGeometryTransformer::fgkGeomTitle     = "Geometry for Agata";
      Char_t  AgataGeometryTransformer::fgDefaultTransp  = 50;

//______________________________________________________________________________
AgataGeometryTransformer::AgataGeometryTransformer()
  : TObject(),
    fMatrixList(new TObjArray(fgkMaxDetectorID)),
    fLog("AgataGeometryTransformer")
{
// Standard constructor
  fMatrixList->SetOwner(true);
}

//______________________________________________________________________________
AgataGeometryTransformer::~AgataGeometryTransformer()
{
// Destructor
  delete fMatrixList;
}

//______________________________________________________________________________
Bool_t AgataGeometryTransformer::ReadTransformations(const TString fileName)
{
  fLog.SetProcessMethod("ReadTransformations(const TString )");

  ifstream in(fileName, ios::in);
  if (!in) {
    fLog << error << Form("file %s not found\n", fileName.Data()) << dolog;
    return false;
  }
  
  Int_t nd, k;
  Double_t x,y,z;
  
  Double_t trans[3];
  Double_t rot[9];
  Char_t line[255];
  
  while ( !in.eof() ) {
    
    in.getline(line,255);
    
    if(sscanf(line, "%d %d %lf %lf %lf", &nd, &k, &x, &y, &z) != 5) {
      fLog << error << "Problem reading translation factors" << dolog;
      return false;
    }
    trans[0] = x; trans[1] = y, trans[2] = z;

    if (nd > fgkMaxDetectorID) {
     fLog << error <<   Form("Id %d too large, max: %d\n", nd, fgkMaxDetectorID) << dolog;
         return false;
    }
    
    in.getline(line,255);
    if(sscanf(line, "%d  %lf %lf %lf", &k, &x, &y, &z) != 4)
    {
      fLog << error << Form("Problem reading 1st rotation factors for id %d\n", nd) << dolog;
      return false;
    }
    rot[0] = x; rot[1] = y; rot[2] = z;
    
    in.getline(line,255);
    if(sscanf(line, "%d  %lf %lf %lf", &k, &x, &y, &z) != 4)
    {
      fLog << error << Form("Problem reading 2nd rotation factors for id %d\n", nd) << dolog;
      return false;
    }
    rot[3] = x; rot[4] = y; rot[5] = z;
    
    in.getline(line,255);
    if(sscanf(line, "%d  %lf %lf %lf", &k, &x, &y, &z) != 4)
    {
      fLog << error << Form("Problem reading 3rd rotation factors for id %d\n", nd) << dolog;
      return false;
    }
    rot[6] = x; rot[7] = y; rot[8] = z;
    
    // add matrix for the given id
    if (fMatrixList->At(nd)) {
      fLog << error << Form("Matrix transformation already exist for id %d\n", nd) << dolog;
      return false;
    } else
      fMatrixList->AddAt( FillTransMatrix(trans, rot), nd );
  }
  
  in.close( );
        
  fLog << dolog;
  
  return true;    
}  

//______________________________________________________________________________
TGeoCombiTrans* AgataGeometryTransformer::FillTransMatrix(Double_t* trans, Double_t* rot)
{
  
  TGeoRotation rotation;  rotation.SetMatrix(rot);
  
  TGeoTranslation translation;  
  translation.SetTranslation(trans[0], trans[1], trans[2]);
  
  return new TGeoCombiTrans(translation, rotation);
}

                 
//_____________________________________________________________________________
void AgataGeometryTransformer::Global2Local(Int_t detID,
                                            Double_t xg, Double_t yg, Double_t zg, 
                                            Double_t& xl, Double_t& yl, Double_t& zl) const
{  
  if (detID < 0 || detID > fgkMaxDetectorID) {
          fLog.SetProcessMethod("Global2Local()");
          fLog << error << "Wrong detector id number" << dolog; 
          return ;
  }
  
  TGeoCombiTrans* mat = static_cast<TGeoCombiTrans*> ( fMatrixList->At(detID) );
  Double_t local[3]  = {0., 0., 0.};
  Double_t global[3] = {xg, yg, zg};
  
  mat->MasterToLocal(global, local);
  xl = local[0];
  yl = local[1];
  zl = local[2];
  
}   

//_____________________________________________________________________________
void AgataGeometryTransformer::Global2Local(Int_t detID,
                                            Float_t xg, Float_t yg, Float_t zg, 
                                            Float_t& xl, Float_t& yl, Float_t& zl) const
{
  if (detID < 0 || detID > fgkMaxDetectorID) {
        fLog.SetProcessMethod("Global2Local()");
    fLog << error << "Wrong detector id number" << dolog; 
        return;
  }
  
  TGeoCombiTrans* mat = static_cast<TGeoCombiTrans*> ( fMatrixList->At(detID) );
  Double_t local[3]  = {0., 0., 0.};
  Double_t global[3] = {xg, yg, zg};
  
  mat->MasterToLocal(global, local);
  xl = local[0];
  yl = local[1];
  zl = local[2];
  
}   


//_____________________________________________________________________________
void AgataGeometryTransformer::Local2Global(Int_t detID,
                                            Double_t xl, Double_t yl, Double_t zl, 
                                            Double_t& xg, Double_t& yg, Double_t& zg) const
{
  if (detID < 0 || detID > fgkMaxDetectorID) {
        fLog.SetProcessMethod("Local2Global()");
        fLog << error << "Wrong detector id number" << dolog; 
        return;
  }
  
  TGeoCombiTrans* mat = static_cast<TGeoCombiTrans*> ( fMatrixList->At(detID) );
  Double_t local[3]  = {xl, yl, zl};
  Double_t global[3] = {0., 0., 0.};
  
  mat->LocalToMaster(local, global);
  xg = global[0];
  yg = global[1];
  zg = global[2];

}   

//_____________________________________________________________________________
void AgataGeometryTransformer::Local2Global(Int_t detID,
                                            Float_t xl, Float_t yl, Float_t zl, 
                                            Float_t& xg, Float_t& yg, Float_t& zg) const
{
  if (detID < 0 || detID > fgkMaxDetectorID) {
        fLog.SetProcessMethod("Local2Global()");
    fLog << error << "Wrong detector id number" << dolog; 
        return ;
  }
  
  TGeoCombiTrans* mat = static_cast<TGeoCombiTrans*> ( fMatrixList->At(detID) );
  Double_t local[3]  = {xl, yl, zl};
  Double_t global[3] = {0., 0., 0.};
  
  mat->LocalToMaster(local, global);
  xg = global[0];
  yg = global[1];
  zg = global[2];
  
}   

//_____________________________________________________________________________
Int_t AgataGeometryTransformer::GetDepthId(Double_t z)
{
  Int_t depthId = 0;
  for (Int_t i = fgkMaxDepthID-1; i >= 0; --i) {
    if (z > fgkDepth[i]) {
      depthId = i;
      break;
    }
  }
  return depthId+1;
}

//_____________________________________________________________________________
Int_t AgataGeometryTransformer::GetSectorId(Double_t x, Double_t y)
{
  // tmp solution
  Double_t rho = TMath::Sqrt(x*x + y*y);
  Double_t phi = TMath::ACos(x/rho);
  
  Double_t angle = TMath::Pi()/3;
  
  Int_t id = Int_t(phi/angle);

  if (y < 0) {
    switch (id) {
      case 0:
        return id+5;
        break;
      case 1:
        return id+3;
        break;
      case 2:
        return id+1;
        break;
    }
  }
  return id;
}

//_____________________________________________________________________________
namespace {
  const Double_t mm = 0.1;
  
  struct ASolid
  {
        Int_t Which;            // Crystal fType defined by an integer
        Int_t Sides;            // Number of faces for this crystal
        
        Double_t HoleR;         // Coaxial hole with radius HoleR   
        Double_t HoleL;         // Hole starts at the depth HoleL
        
        Double_t CylR;          // Radius of the cylinder that is the base of the shaped crystal     
        Double_t CylL;          // Length of the cylinder that is the base of the shaped crystal     
        Double_t CylX;          // Additionnal offset 
        Double_t CylY;          // Additionnal offset 
        Double_t CylZ;          // Additionnal offset
        
        Double_t ThickB;        // Thickness of the passive area (back)   
        Double_t ThickC;        // Thickness of the passive area (coaxial)
        Double_t CapS;          // The crystal-encapsulation spacing is capS
        Double_t CapT;          // The capsule is capT thick
        Double_t Tolerance;     // CapS+CapT+tol
        
        Float_t ColX;           // RGB color components    
        Float_t ColY;           // RGB color components   
        Float_t ColZ;           // RGB color components
        
        TObjArray *Vertex;      // Array of 3D points (TVector3)
        
        ASolid() : Which(-1),Sides(0),
        HoleR(0.0),HoleL(0.0),CylR(0.0),CylL(0.0),CylX(0.0),CylY(0.0),CylZ(0.0),
        ThickB(0.0),ThickC(0.0),CapS(0.0),CapT(0.0),Tolerance(0.0),
        ColX(0.0),ColY(0.0),ColZ(0.0) { Vertex = new TObjArray(); Vertex->SetOwner() ; }
        ~ASolid() { Vertex->Delete(); delete Vertex; }
        
        void Print() const;
  };
  
  void ASolid::Print() const
  {
        cout << " Crystal Type " << Which << endl;
        cout << " Sides " << Sides << endl;
        cout << " Dimensions crystal " << CylL << " " << CylR << endl;  
        cout << " Dimensions hole " << HoleL << " " << HoleR << endl;           
  }
  
} // anonymous namespace for local declarations

//_____________________________________________________________________________
Int_t AgataGeometryTransformer::BuildAgataCrystals(const char *filename, const char *basecrysname)
{
  if ( gGeoManager == 0x0 ) 
  { cout << " BuildAgataCrystals needs an existing gGeoManager  ... " << endl; return 0 ; }
  
  // try to open the file
  ifstream infil; infil.open(filename); if ( !infil.is_open() ) return 0 ;
  
  cout << " ---> Reading description of crystals from file " << filename << " ..." << endl;
  
  // read the file (see AgataDetectorArray::ReadSolidFile())
  int i1,i2,i3,lline; double x,y,z,X,Y,Z; char line[255];
  
  std::vector<ASolid *> pgons; ASolid *pPg = 0x0; int nPgons = 0; // number of crystals defined in the file
  
  int opgon = -1;
  while( infil.good() ) {
        infil.getline(line,255); lline = strlen(line);
        
        if ( lline < 2 ) continue; if ( line[0] == '#' ) continue;
        
        // decode the line
        if(sscanf(line,"%d %d %d %lf %lf %lf %lf %lf %lf", &i1, &i2, &i3, &x, &y, &z, &X, &Y, &Z) != 9) { 
          nPgons++;
          break;
        } 
        if(opgon != i1) { // a new crystal is being defined, so init a new ASolid structure
          nPgons++; 
          opgon = i1;
          pPg = new ASolid(); pPg->Which = i1; pPg->Sides  = i2; pgons.push_back( pPg );        
        }
        if(i2==0 && i3==0) { // basic shape for the crystal
          pPg->HoleR = x * mm;
          pPg->CylR  = y * mm;
          pPg->CylL  = z * mm;
          pPg->CylX  = X * mm;
          pPg->CylY  = Y * mm;
          pPg->CylZ  = Z * mm;
        }
        else if(i2==0 && i3==1) { // passive and capsule
          pPg->HoleL  = x * mm;
          pPg->ThickB = y * mm;
          pPg->ThickC = z * mm;
          pPg->CapS   = X * mm;
          pPg->CapT   = Y * mm;
          pPg->Tolerance   = Z * mm;    
        }
        else if(i2==0 && i3==2) { // colors 
          pPg->ColX     = x;
          pPg->ColY     = y;
          pPg->ColZ     = z;
        }
        else { // a new point to define the ploygon
          pPg->Vertex->AddAt(new TVector3( x * mm, y * mm,  z * mm ),   i3);
          pPg->Vertex->AddAt(new TVector3( X * mm, Y * mm,  Z * mm ),i3+i2);
        }
  } // while good
  infil.close();
  cout << " ... " << nPgons << " different crystals have been built  " << endl; // end read the file. 
  
  // now it built the different crystals and add them to the TGeoManager
  TGeoMaterial *mat;
  TGeoMedium *med;
  
  if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Ge")) == 0x0 )
        mat = new TGeoMaterial("Ge", 72.59,32,5.32); 
  
  if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Ge")) == 0x0 )
        med = new TGeoMedium("Ge",1,mat);
  
  vector<ASolid*>::iterator itthesolids;
  for (itthesolids = pgons.begin(); itthesolids!=pgons.end(); itthesolids++) { 
        
        pPg = *itthesolids; TString cname, pname;  // to get the current crystal definition and set names
        
        // Base shape for this crystal
        // Creation of a Pcon with a name based on the crystal number written in the asolid file
        cname = "BaseShape"; 
        cname += pPg->Which; 
        TGeoPcon *baseshape = new TGeoPcon(0.0,360.0,4); 
        
        baseshape->SetName(cname.Data());
        baseshape->DefineSection(0,0.0,0.0,pPg->CylR);
        baseshape->DefineSection(1,0.9990*pPg->HoleL,0.0,pPg->CylR);
        baseshape->DefineSection(2,pPg->HoleL,pPg->HoleR,pPg->CylR);
        baseshape->DefineSection(3,pPg->CylL,pPg->HoleR,pPg->CylR);
                
        // Method 3 to build the crystal shape  
        Double_t scale = 0.0, xv[12], yv[12];
        TVector3 *p1, *p2;                
        for (int i = 0; i < pPg->Sides; i++ ) {
          
          p1 = (TVector3 *)pPg->Vertex->At(i); p2 = (TVector3 *)pPg->Vertex->At(i+pPg->Sides);  
          xv[i] = p1->x(); yv[i] = p1->y();
          
          Double_t s = TMath::Sqrt(p2->x()*p2->x()+p2->y()*p2->y())
          / TMath::Sqrt(p1->x()*p1->x()+p1->y()*p1->y()) ;
          
          if ( s > scale ) scale = s;
        }
        
        // now build the union as a composition shape that is the shaper of the raw crystal
        pname = "Shaper"; pname += pPg->Which; 
        TGeoXtru *shaper = new TGeoXtru(2); shaper->SetName(pname.Data());
        
        shaper->DefinePolygon(pPg->Sides,xv,yv);
        shaper->DefineSection(0,-0.001,0.0,0.0,1.0);
        shaper->DefineSection(1,1.001*pPg->CylL,0.0,0.0,scale); 
        
        // Intersection between the raw crystal and the shaper
        cname += " * "; cname += "Shaper"; cname += pPg->Which; 
        
        // methode 2
        // cname += ":"; cname += "TrShaper"; cname += pPg->Which;
        
        pname = basecrysname; pname += pPg->Which;      
        TGeoCompositeShape *crysshape = new TGeoCompositeShape(pname.Data(),cname.Data());
        
        //              TGeoVolume *thecrys = new TGeoVolume(pname.Data(),arb8,med);
        //              TGeoVolume *thecrys = new TGeoVolume(pname.Data(),shaper,med);
        TGeoVolume *thecrys = new TGeoVolume(pname.Data(),crysshape,med);               
        //              TGeoVolume *thecrys = new TGeoVolume(pname.Data(),baseshape,med);               
        thecrys->SetLineColor(TColor::GetColor(pPg->ColX,pPg->ColY,pPg->ColZ)); 
    thecrys->SetTransparency(fgDefaultTransp);
  } // itthesolids
  
  return  nPgons;       
}

//_____________________________________________________________________________
Int_t AgataGeometryTransformer::BuildAgataClusters(const char *filename, const char *baseclusname, const char *basecrysname)
{
  // crystals are not yet defined so the cluster cannot be built
  if ( gGeoManager == 0x0 ) { 
        cout << " BuildAgataClusters needs an existing gGeoManager  ... " << endl; 
        return 0;
  }
  
  // try to open the file
  ifstream infil; infil.open(filename); if ( !infil.is_open() ) return 0 ;
  
  cout << "\nReading description of clusters from file " << filename << " ..." << endl;
  
  Char_t    line[256];
  Int_t     lline, i1, i2, i3, oclust, nclust; Double_t  ps, th, ph, x, y, z;
  
  TString crysname, clusname; TGeoVolume *clu = 0x0, *cris = 0x0;
  
  nclust = 0; oclust  = -1;     
  while(infil.good()) {
        infil.getline(line,255); lline = strlen(line);
        if(lline < 2) continue;
        if(line[0] == '#') continue;
        if(sscanf(line,"%d %d %d %lf %lf %lf %lf %lf %lf", &i1, &i2, &i3, &ps, &th, &ph, &x, &y, &z) != 9) {
          break;
        }
        if(oclust != i1) {
          oclust = i1; 
          nclust++;
          
          // a new assembly is created to define the cluster
          clusname = baseclusname; clusname += oclust; clu = new TGeoVolumeAssembly(clusname.Data());
        }
        // look for the corresponding crystal in the list of existing volumes
        crysname = basecrysname; crysname += i2; cris = gGeoManager->GetVolume(crysname.Data());
        if ( cris == 0x0 ) {
          cout << "Cannot build correctly the cluster " << clusname << ", " << crysname << " does not exist " << endl;
        }
        else { // add this crystal after transformation to the assembly
          TGeoHMatrix H; TGeoRotation R; TGeoTranslation T(x*mm,y*mm,z*mm);
          R.RotateZ(ps);
          R.RotateY(th);
          R.RotateZ(ph);        
          
          H = T * R; 
          TGeoHMatrix *h = new TGeoHMatrix(H); clu->AddNode(cris,i3,h);
        }
  }
  infil.close();
  
  cout << nclust << " different cluster(s) have been build " << endl; return nclust;    
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::BuildAgataDetector(const char *filename, const char *baseclusname, const char *agataName)
{
  // crystals are not yet defined so the cluster cannot be built
  if ( gGeoManager == 0x0 ) {  
        cout << " BuildAgata needs an existing gGeoManager  ... " << endl; 
        return 0x0; 
  }
  
  // try to open the file
  ifstream infil; infil.open(filename); if ( !infil.is_open() ) return 0x0 ;
  
  cout << "\n Reading description of AGATA from file " << filename << " ..." << endl;
  
  TString clusname;
  Char_t    line[256];
  Int_t     lline, i1, i2, nclust = 0; Double_t ps, th, ph, x, y, z;
  
  TGeoVolume *agata = 0x0;
  while(infil.good())  {
        infil.getline(line,255); lline = strlen(line);
        if(lline < 2) continue;
        if(line[0] == '#') continue;
        if(sscanf(line,"%d %d %lf %lf %lf %lf %lf %lf", 
                          &i1, &i2, &ps, &th, &ph, &x, &y, &z) != 8) break;
        
        clusname = baseclusname; clusname += i2;
        TGeoVolume *clus = gGeoManager->GetVolume(clusname.Data());     
        if ( clus ) { // this cluster is known
          TGeoMedium   *med;
          TGeoMaterial *mat;
          if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
                mat = new TGeoMaterial("Vacuum",0,0,0);                         
          if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
                med = new TGeoMedium("Vacuum",1,mat);
          
          if ( agata == 0x0 )
                agata = gGeoManager->MakeBox(agataName,med,60,60,60); // volume corresponding to agata
          
          TGeoHMatrix H; TGeoRotation R; TGeoTranslation T(x*mm,y*mm,z*mm);
          R.RotateZ(ps);
          R.RotateY(th);
          R.RotateZ(ph);        
          
          H = T * R; 
          TGeoHMatrix *h = new TGeoHMatrix(H); agata->AddNode(clus,i1,h); nclust++;
        }
  }
  infil.close();
  
  cout << nclust << " cluster(s) composed this geometry " << endl; return agata;        
}       

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::ImportAgata(const char *asolid, const char *aclust, const char *aeuler, 
                                                                                                  const char *baseclusname, const char *agataName)
{  
  if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
        new TGeoManager(fgkGeomName, fgkGeomTitle);
  }
  
  // try to build the different part of the detector
  TGeoVolume *agata = 0x0;
  if ( BuildAgataCrystals(asolid) > 0 ) { // at least one crystal has been build
        if ( BuildAgataClusters(aclust) > 0 ) { // at least one cluster has been build
          agata = BuildAgataDetector(aeuler, baseclusname, agataName);
        }
  }

  return agata;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::BuildAgata(const char* agataPath, const char *asolid, const char *aclust, const char *aeuler, 
                                                                                                  const char *baseclusname, const char *agataName)
{     
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   TString path(agataPath);
   if (!path.EndsWith("/"))
          path += "/";
   TString geoSolid(asolid);
   TString geoCluster(aclust);
   TString geoEuler(aeuler);
   
   geoSolid.Prepend(path);
   geoCluster.Prepend(path);
   geoEuler.Prepend(path);
   
   // try to build the different part of the detector
   TGeoVolume *agata = 0x0;
   if ( BuildAgataCrystals(geoSolid.Data()) > 0 ) { // at least one crystal has been build
          if ( BuildAgataClusters(geoCluster.Data()) > 0 ) { // at least one cluster has been build
                 agata = BuildAgataDetector(geoEuler.Data(), baseclusname, agataName);
          }
   }
      
   return agata;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::BuildDante(const char* danteConf, const char* basemoduleName, const char *danteName)
{
   std::ifstream in(danteConf); std::string tmp;
   
   if ( in.is_open() == false ) {
          cout << " ---> Problem opening file " << danteConf << " ..." << endl;
          return 0x0;
   }
   
   cout << " ---> Reading description of Dante from file " << danteConf << " ..." << endl;

   getline(in,tmp);
   tmp += " ";
   TGeoVolume* dante = 0x0;
   while ( in.good() && !in.eof() ) { // read input stream line by line
          
          if ( tmp[0] == '#' ) {
                 getline(in,tmp); 
                 continue; 
          } // this line is a comment 
          
          Float_t dz, rz1, ry, rz2, rx, tx1, tx2,ty1, ty2, ax, bx, cx, ay, by, cy ;
          Int_t which; 
          std::istringstream decode(tmp);
          
          decode >> which;
          
          if (which != -1) {
                 
                 decode >> dz >> rz1 >> ry >> rz2 >> rx >> tx1 >> tx2 >> ty1 >> ty2 >> ax >> bx >> cx >> ay >> by >> cy ;
                 
                 dante = AgataGeometryTransformer::AddDanteModule(dz/10., rz1, ry, rz2, rx, basemoduleName, danteName);
                 
                 // add 0x0 character at the end since getline
                 // does not store \n and it may give troubles when decoding
          }      
          getline(in,tmp);
          tmp += " ";
          
   }    
   in.close();
      

   return dante;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddDanteModule(const Float_t dZ, const Float_t rotZ1, const Float_t rotY, const Float_t rotZ2, 
                                                                                                         const Float_t rotX, const char* basemoduleName, const char *danteName)
{
  if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
        new TGeoManager(fgkGeomName, fgkGeomTitle);
  }
  
  TGeoVolume* dante = gGeoManager->FindVolumeFast(danteName);
  if ( dante == 0x0 ) {
        TGeoMedium   *med;
        TGeoMaterial *mat;
        if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
          mat = new TGeoMaterial("Vacuum",0,0,0);                       
        if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
          med = new TGeoMedium("Vacuum",1,mat);
        
        dante = gGeoManager->MakeBox(danteName,med,60,60,60); // volume corresponding to Dante
  } 

  // create module
  TGeoMaterial* matMod;
  TGeoMedium*   medMod;
  
  if ( (matMod = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Si")) == 0x0 )
        matMod = new TGeoMaterial("Si", 28.09, 14, 2.3);
  if ( (medMod = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Si")) == 0x0 )
        medMod = new TGeoMedium("Si",2,matMod);

  TGeoBBox *box = new TGeoBBox(Form("%s_Box",basemoduleName), 2., 3., 0.2);
  
  // Apply transfo
  TGeoRotation rot;
  TGeoHMatrix  transfo;
  TGeoHMatrix* hm = 0x0;
  TGeoTranslation trans(0., 0., dZ);
  
  rot.RotateZ(rotZ1);
  rot.RotateY(rotY);
  rot.RotateZ(rotZ2);
  rot.RotateX(rotX);
  
  transfo = rot * trans;
  hm = new TGeoHMatrix(transfo);

  Int_t nbModule = 0;

  TGeoVolume *danteMod = new TGeoVolume(Form("%s_Dante",basemoduleName),box, medMod);
  danteMod->SetLineColor(kAzure-5);
  danteMod->SetTransparency(fgDefaultTransp);

   TObjArray* list = dante->GetNodes();
  if (list) {
        for (Int_t i = 0; i < list->GetEntries(); ++i) {
          TGeoVolume* vol = (TGeoVolume*)list->At(i);
          if (vol) {
                TString name(vol->GetName());
                if ( name.Contains(Form("%s_Dante",basemoduleName)) )
                        nbModule++;
          }
        }
  }
  
  dante->AddNode(danteMod, nbModule, hm);
  return dante;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddBeamPipe(const Float_t dZ, const Float_t minR, const Float_t maxR, const char *pipeName )
{
  if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
        new TGeoManager(fgkGeomName, fgkGeomTitle);
  }
  
  // create module
  TGeoMaterial* matPipe;
  TGeoMedium*   medPipe;
  
  if ( (matPipe = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Al")) == 0x0 )
        matPipe = new TGeoMaterial("Al", 26.98,13,2.7);
  if ( (medPipe = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Al")) == 0x0 )
        medPipe = new TGeoMedium("Al", 2, matPipe);
  
   
  TGeoVolume* tube = gGeoManager->MakeTube(pipeName, medPipe, minR, maxR, dZ);
  tube->SetVisibility(true);
  tube->SetLineColor(17);
  tube->SetTransparency(fgDefaultTransp);
  
  return tube;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddTarget(const Float_t dx, const Float_t dy, const Float_t dz, 
                                                                                                const char *targetName, TGeoMaterial* mat)
{
  if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
        new TGeoManager(fgkGeomName, fgkGeomTitle);
  }
  
  // create module
  TGeoMaterial* matTarget;
  TGeoMedium*   medTarget;
  
  if (mat != 0) {
        matTarget = mat;
        medTarget = new TGeoMedium("Medium", 1, matTarget);
  } else {
        if ( (matTarget = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
          matTarget = new TGeoMaterial("Vacuum", 0., 0., 0.);
        if ( (medTarget = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
          medTarget = new TGeoMedium("Vacuum", 1, matTarget);
  }
  
  TGeoVolume* target = gGeoManager->MakeBox(targetName, medTarget, dx, dy, dz);
  target->SetVisibility(true);
  
  return target;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddFloor(Float_t dx, const Float_t dy, const Float_t dz, const char *floorName )
{
  if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
        new TGeoManager(fgkGeomName, fgkGeomTitle);
  }
  
  // create module
  TGeoMaterial* matFloor;
  TGeoMedium*   medFloor;
  
  if ( (matFloor = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Al")) == 0x0 )
        matFloor = new TGeoMaterial("Al", 26.98,13,2.7);
  if ( (medFloor = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Al")) == 0x0 )
        medFloor = new TGeoMedium("Al", 2, matFloor);
  
  
  TGeoVolume* floor = gGeoManager->MakeBox(floorName, medFloor, dx, dy, dz);
  floor->SetVisibility(true);
  floor->SetLineColor(kWhite);
  floor->SetTransparency(fgDefaultTransp);
  return floor;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddPrismaMCP(const char *prismaName, Float_t dx, Float_t dy, Float_t dz, Float_t phi, TGeoMaterial* mat)
{
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   // cmedMWPCeate module
   TGeoMaterial* matMCP;
   TGeoMedium*   medMCP;
   
   if (mat != 0) {
          matMCP = mat;
          medMCP = new TGeoMedium("Medium", 1, matMCP);
   } else {
          if ( (matMCP = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
                 matMCP = new TGeoMaterial("Vacuum", 0., 0., 0.);
          if ( (medMCP = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
                 medMCP = new TGeoMedium("Vacuum", 1, matMCP);
   }
   
   
   TGeoVolume* mcp = gGeoManager->MakeBox(prismaName,medMCP,dx,dy,dz); 

   TGeoVolume* start = gGeoManager->MakeBox("start", medMCP, dx, dy, dz);
   start->SetVisibility(true);
   start->SetLineColor(kBlue);
   start->SetTransparency(fgDefaultTransp);
   
   TGeoRotation rot;
   rot.RotateY(phi);
   TGeoHMatrix *hm = new TGeoHMatrix(rot);
   mcp->AddNode(start, 0, hm);
   
   
   return mcp;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddPrismaDipoleOld(const char *prismaName, Float_t dz, Float_t rmin, Float_t rmax, Float_t phi1, Float_t phi2)
{
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   TGeoMaterial* matPrismaDipole;
   TGeoMedium*   medPrismaDipole;
   
   if ( (matPrismaDipole = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Fe")) == 0x0 )
          matPrismaDipole = new TGeoMaterial("Fe", 55.9, 26, 7.87);
   if ( (medPrismaDipole = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Fe")) == 0x0 )
          medPrismaDipole = new TGeoMedium("Fe", 3, matPrismaDipole);
   
   TGeoMedium   *med;
   TGeoMaterial *mat;
   if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
          mat = new TGeoMaterial("Vacuum",0,0,0);                       
   if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
          med = new TGeoMedium("Vacuum",1,mat);
   
   TGeoRotation rot;
   rot.RotateX(90);  
   rot.RotateY(65);  
   
   TGeoHMatrix  transfo;
   TGeoTranslation trans(-(rmax+rmin)/2*(1+(1-TMath::Cos(20.*TMath::DegToRad()))), 0., (rmax+rmin)/2);
   transfo = trans*rot;
   TGeoHMatrix *hm = new TGeoHMatrix(transfo);
   
   TGeoVolume* dipole = gGeoManager->MakeCons(Form("%s_%s", prismaName,"Dipole"), medPrismaDipole, dz, rmin, rmax, rmin, rmax, phi1, phi2);
   dipole->SetVisibility(true);
   dipole->SetLineColor(55);
   dipole->SetTransparency(fgDefaultTransp);
   
   TGeoConeSeg* shape = (TGeoConeSeg*)dipole->GetShape();
   shape->ComputeBBox();
   Float_t boxSizeX = shape->GetDX();
   Float_t boxSizeY = shape->GetDY();
   Float_t boxSizeZ = shape->GetDZ();
   
   TGeoVolume* prisma = gGeoManager->MakeBox(prismaName, med, boxSizeX, boxSizeY, boxSizeZ); // volume corresponding to prisma dipole
   prisma->AddNode(dipole, 1, hm);
   
   return prisma;
}


//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddPrismaDipole(const char *prismaName, Float_t dx, Float_t dy, Float_t dz)
{
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }

   TGeoMaterial* matPrismaDipole;
   TGeoMedium*   medPrismaDipole;
   
   if ( (matPrismaDipole = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Fe")) == 0x0 )
          matPrismaDipole = new TGeoMaterial("Fe", 55.9, 26, 7.87);
   if ( (medPrismaDipole = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Fe")) == 0x0 )
          medPrismaDipole = new TGeoMedium("Fe", 3, matPrismaDipole);
   
   TGeoMedium   *med;
   TGeoMaterial *mat;
   if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
          mat = new TGeoMaterial("Vacuum",0,0,0);                       
   if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
          med = new TGeoMedium("Vacuum",1,mat);

   // Dipole entrance
   TGeoRotation rot1;  
   rot1.RotateY(20);  
   
   TGeoHMatrix  transfo1;
   transfo1 = rot1;
   TGeoHMatrix *hm1 = new TGeoHMatrix(transfo1);
   
   TGeoVolume* dipoleIn = gGeoManager->MakeBox(Form("%s_%s", prismaName,"Dipole"), medPrismaDipole, dx, dy, dz);
   dipoleIn->SetVisibility(true);
   dipoleIn->SetLineColor(kViolet);
   dipoleIn->SetTransparency(fgDefaultTransp);
   
   // Dipole exit
   TGeoRotation rot2;  
   rot2.RotateY(-60);  
   Float_t posX = 60.;
   Float_t posZ = 60.*TMath::Tan(60.*TMath::DegToRad());
   TGeoTranslation trans2(-posX, 0, posZ);
   
   TGeoHMatrix  transfo2;
   transfo2 = trans2*rot2;
   TGeoHMatrix *hm2 = new TGeoHMatrix(transfo2);
   
   TGeoVolume* dipoleOut = gGeoManager->MakeBox(Form("%s_%s", prismaName,"Dipole"), medPrismaDipole, dx, dy, dz);
   dipoleOut->SetVisibility(true);
   dipoleOut->SetLineColor(kViolet);
   dipoleOut->SetTransparency(fgDefaultTransp);
   
   TGeoVolume* prisma = gGeoManager->MakeBox(prismaName, med, 200, 200, 200); // volume corresponding to prisma dipole
   prisma->AddNode(dipoleIn,  1, hm1);
   prisma->AddNode(dipoleOut, 2, hm2);

   return prisma;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddPrismaQpole(const char *prismaName, Float_t dz, Float_t rmin, Float_t rmax)
{
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   TGeoMaterial* matPrismaQpole;
   TGeoMedium*   medPrismaQpole;
   
   if ( (matPrismaQpole = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Fe")) == 0x0 )
          matPrismaQpole = new TGeoMaterial("Fe", 55.9, 26, 7.87);
   if ( (medPrismaQpole = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Fe")) == 0x0 )
          medPrismaQpole = new TGeoMedium("Fe", 3, matPrismaQpole);
   
   TGeoVolume* prisma = gGeoManager->MakeTube(prismaName, medPrismaQpole, rmin, rmax, dz);
   prisma->SetVisibility(true);
   prisma->SetLineColor(kViolet);
   prisma->SetTransparency(fgDefaultTransp);
   
   return prisma;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::AddPrismaMWPC(const char *prismaName, Float_t dzFcIc, Float_t dx, Float_t dy, Float_t dz, TGeoMaterial* mat)
{
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   // cmedMWPCeate module
   TGeoMaterial* matMWPC;
   TGeoMedium*   medMWPC;
   
   if (mat != 0) {
          matMWPC = mat;
          medMWPC = new TGeoMedium("Medium", 1, matMWPC);
   } else {
          if ( (matMWPC = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
                 matMWPC = new TGeoMaterial("Vacuum", 0., 0., 0.);
          if ( (medMWPC = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
                 medMWPC = new TGeoMedium("Vacuum", 1, matMWPC);
   }
   
 
   TGeoVolume* mwpc = gGeoManager->MakeBox(prismaName,medMWPC,dx,dy,100); 

   TGeoVolume* fc = gGeoManager->MakeBox("FC", medMWPC, dx, 10, 1);
   fc->SetVisibility(true);
   fc->SetLineColor(kBlue);
   fc->SetTransparency(fgDefaultTransp);
   mwpc->AddNode(fc,0);
   
   TGeoVolume* ic = gGeoManager->MakeBox("IC", medMWPC, dx, 10, dz);
   ic->SetVisibility(true);
   ic->SetLineColor(kBlue);
   ic->SetTransparency(fgDefaultTransp);
   TGeoTranslation trans1(0., 0., dzFcIc+dz);
   TGeoHMatrix *hm1 = new TGeoHMatrix(trans1);
   mwpc->AddNode(ic, 1, hm1);
   
 
   return mwpc;
}

//_____________________________________________________________________________
TGeoVolume* AgataGeometryTransformer::BuildPrisma(const char* prismaConf, Float_t dzMCP, Float_t dzQpole, Float_t dzDipole, Float_t dzDipMWPC, 
                                                                                                  const char* nameMCP, const char* nameQPole, const char* nameDipole, const char* nameMWPC, 
                                                                                                  const char* prismaName)
{     
   if ( gGeoManager == 0x0 ) { // a new Geo Manager is created if needed
          new TGeoManager(fgkGeomName, fgkGeomTitle);
   }
   
   TGeoMedium   *med;
   TGeoMaterial *mat;
   if ( (mat = (TGeoMaterial *)gGeoManager->GetListOfMaterials()->FindObject("Vacuum")) == 0x0 )
          mat = new TGeoMaterial("Vacuum",0,0,0);                       
   if ( (med = (TGeoMedium *)gGeoManager->GetListOfMedia()->FindObject("Vacuum")) == 0x0 )
          med = new TGeoMedium("Vacuum",1,mat);
   

   Float_t dzMCPLoc; Float_t dzQpoleLoc; Float_t dzDipoleLoc; Float_t dzDipMWPCLoc;
   if ( ReadPrismaConf(dzMCPLoc, dzQpoleLoc, dzDipoleLoc, dzDipMWPCLoc, prismaConf) ) {
          dzMCP     = dzMCPLoc; 
          dzQpole   = dzQpoleLoc; 
          dzDipole  = dzDipoleLoc; 
          dzDipMWPC = dzDipMWPCLoc;
          cout << " ---> Reading description of Prisma from file " << prismaConf << " ..." << endl;
   } else {
          cout << " ---> Using default values for Prisma geometry" << " ..." << endl;  
   }
   
   // global volume
   Float_t dyRot = 30.;
   TGeoVolume* prisma = gGeoManager->MakeBox(prismaName,med,450,450,450); 

   // MCP
   TGeoVolume* mcp = AgataGeometryTransformer::AddPrismaMCP(nameMCP);
   TGeoTranslation trans0(0., 0., dzMCP);
   TGeoHMatrix *hm0 = new TGeoHMatrix(trans0);
   
   prisma->AddNode(mcp, 0, hm0);
   
   // Qpole   
   TGeoVolume* Qpole = AgataGeometryTransformer::AddPrismaQpole(nameQPole);
   TGeoTube* shape = (TGeoTube*)Qpole->GetShape(); 
   shape->ComputeBBox();
   Float_t dzQp   = shape->GetDZ();
   Float_t posZQp = dzQpole + dzQp;

   TGeoTranslation trans1(0., 0., posZQp);
   TGeoHMatrix *hm1 = new TGeoHMatrix(trans1);

   prisma->AddNode(Qpole, 1, hm1);
   
   // Dipole
   Float_t posZDip = dzDipole;
   
   TGeoVolume* dipole = AgataGeometryTransformer::AddPrismaDipole(nameDipole);
   TGeoHMatrix  transfo2;
   TGeoTranslation trans2(0., 0.,posZDip);
   transfo2 = trans2;
   TGeoHMatrix *hm2 = new TGeoHMatrix(transfo2);
   prisma->AddNode(dipole, 2, hm2);
   
   //MWPC     
   Float_t posXmwpc = 60 + dzDipMWPC * TMath::Cos(dyRot*TMath::DegToRad());
   Float_t posZmwpc = posZDip + 60.*TMath::Sqrt(3.) + dzDipMWPC * TMath::Sin(dyRot*TMath::DegToRad());
   
   TGeoRotation rot3;
   TGeoHMatrix  transfo3;
   Float_t dphi = dyRot-90;
   TGeoTranslation trans3(-posXmwpc, 0, posZmwpc);
   rot3.RotateY(dphi);
   transfo3 = trans3*rot3;
   TGeoHMatrix *hm3 = new TGeoHMatrix(transfo3);
   TGeoVolume* mwpc = AgataGeometryTransformer::AddPrismaMWPC(nameMWPC);        
   prisma->AddNode(mwpc, 3, hm3);
   
   return prisma;
}

//_____________________________________________________________________________
Bool_t AgataGeometryTransformer::ReadPrismaConf(Float_t& dzMCP, Float_t& dzQpole, Float_t& dzDipole, Float_t& dzDipMWPC, const char* prismaConf)
{
   std::map<TString, float> map;
   ifstream in(prismaConf);
   if ( !in ) {
          return false;
   }
   
   Char_t line[255];
   while ( in.getline(line,255) ) {
          
          TString tmp(line);

          if (tmp[0] == '#') continue;
          if (tmp[0] == ' ') continue;
          if (tmp.Length() < 5)  continue;

          Int_t pos = tmp.First('=');
          
          TString sValue = tmp(pos+1, tmp.Length()-pos);
          
          TString sToken = tmp(0, pos-1);
          
          Float_t value = sValue.Atof();
          map[sToken] = value;
   }

   std::map<TString, float>::iterator iter;
   
   iter = map.find(TString("target_mcp_distance"));  
   dzMCP = iter->second/10.;
   
   iter = map.find(TString("target_quad_distance"));  
   dzQpole = iter->second/10.;
   
  iter = map.find(TString("target_dipole_distance"));  
   dzDipole = iter->second/10.;
   
   iter = map.find(TString("dip_fp_dist"));  
   dzDipMWPC= iter->second/10.;
   
   return true;
}