FairModule.cxx

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairModule.cxx

// -------------------------------------------------------------------------
// -----                     FairModule source file                    -----
// -----            Created 06/01/04  by M. Al-Turany                  -----
// -------------------------------------------------------------------------
/* Generated by Together */

#include "FairModule.h"


#include "FairVolume.h"
#include "FairVolumeList.h"
#include "FairBaseParSet.h"
#include "FairRun.h"

#include "FairGeoNode.h"
#include "FairRuntimeDb.h"

#include "FairGeoInterface.h"
#include "FairGeoLoader.h"
#include "FairGeoNode.h"
#include "FairGeoRootBuilder.h"
#include "FairGeoMedia.h"


#include "TString.h"
#include "TObjArray.h"
#include "TGeoVolume.h"
#include "TFile.h"
#include "TList.h"
#include "TKey.h"
#include "TGeoManager.h"
#include "TGeoVoxelFinder.h"
#include "TGeoMatrix.h"


#include "TSystem.h"
#include <fstream>
#include <iostream>


TArrayI* FairModule::volNumber=0;
Int_t FairModule::fNbOfVolumes=0;
FairVolumeList*  FairModule::vList=0;
TRefArray*    FairModule::svList=0;



//__________________________________________________________________________
void FairModule::ConstructGeometry()
{
  fLogger->Warning(MESSAGE_ORIGIN,"FairModule::ConstructGeometry() : this method has to be implimented in detector class");

}
//__________________________________________________________________________
void FairModule::ConstructOpGeometry()
{
  fLogger->Debug2(MESSAGE_ORIGIN,"FairModule::ConstructOpGeometry : this method has to be implimented in detector class");

}
//__________________________________________________________________________
FairModule::~FairModule()
{

}
//__________________________________________________________________________
FairModule::FairModule(const char* Name, const char* title ,Bool_t Active)
  :TNamed(Name, title),
   fMotherVolumeName(""),
   fgeoVer("Not defined"),
   fgeoName("Not defined"),
   fModId(-1),
   fActive(Active),
   fNbOfSensitiveVol(0),
   fVerboseLevel(0),
   flGeoPar(0),
   kGeoSaved(kFALSE),
   fRootMaterialImportFlag(kFALSE),
   fLogger(FairLogger::GetLogger())
{
  if(!svList) { svList=new TRefArray(); }
  if(!vList) { vList=new FairVolumeList(); }

}

//__________________________________________________________________________

FairModule::FairModule()
  : TNamed(),
    fMotherVolumeName(""),
    fgeoVer("Not defined"),
    fgeoName("Not defined"),
    fModId(-1),
    fActive(kFALSE),
    fNbOfSensitiveVol(0),
    fVerboseLevel(0),
    flGeoPar(0),
    kGeoSaved(kFALSE),
    fRootMaterialImportFlag(kFALSE),
    fLogger(FairLogger::GetLogger())
{

}

//__________________________________________________________________________
void FairModule::Streamer(TBuffer& b)
{
  TNamed::Streamer(b);


  if (b.IsReading()) {
    fgeoVer.Streamer(b);
    fgeoName.Streamer(b);
    b >> fActive;
    b >> fModId;
  } else {
    fgeoVer.Streamer(b);
    fgeoName.Streamer(b);
    b << fActive;
    b << fModId;
  }

}
//__________________________________________________________________________
void FairModule::SetGeometryFileName(TString fname, TString geoVer)
{
  fgeoVer=geoVer;
  TString FileName = fname;
  TString work = getenv("VMCWORKDIR");
  TString userwork = getenv("GEOMPATH");
  if(userwork != "") {
    fgeoName=userwork;
    if (!fgeoName.EndsWith("/")) { fgeoName+="/"; }
    //fgeoName+=fname;
    if (TString(gSystem->FindFile(fgeoName.Data(),fname)) != TString("")) {
      fgeoName=fname;
      fLogger->Info(MESSAGE_ORIGIN, "User path for detector geometry : %s ", fgeoName.Data());
    } else {
      fLogger->Warning(MESSAGE_ORIGIN, "Detector geometry was not found in user path : %s ", FileName.Data());
      fgeoName=work+"/geometry/";
      fLogger->Info(MESSAGE_ORIGIN, "Try the standard path : %s ",fgeoName.Data());
      if (TString(gSystem->FindFile(fgeoName.Data(),FileName)) != TString("")) {
        fgeoName=FileName;
        fLogger->Info(MESSAGE_ORIGIN, "Reading detector geometry from : %s ", FileName.Data());
      } else {
        fLogger->Fatal(MESSAGE_ORIGIN,"Detector geometry not found.");
      }
    }
  } else {
    fgeoName=work+"/geometry/";
    fgeoName+=fname;
  }
}
//__________________________________________________________________________
void FairModule::ProcessNodes(TList* aList)
{

  TListIter iter(aList);
  FairGeoNode* node   = NULL;
  FairGeoNode* MotherNode =NULL;
  FairVolume*  volume = NULL;
  FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb();
  FairBaseParSet* par=(FairBaseParSet*)(rtdb->getContainer("FairBaseParSet"));
  TObjArray* fNodes = par->GetGeoNodes();
  while( (node = (FairGeoNode*)iter.Next()) ) {

    node->calcLabTransform();
    MotherNode=node->getMotherNode();
    volume = new FairVolume( node->getTruncName(), fNbOfVolumes++);
    volume->setRealName(node->GetName());
    vList->addVolume(volume);
    volume->setGeoNode(node);
    volume->setCopyNo(  node->getCopyNo());

    if(MotherNode!=0) {
      volume->setMotherId(node->getMCid());
      volume->setMotherCopyNo(MotherNode->getCopyNo());
    }
    FairGeoVolume* aVol=NULL;

    if ( node->isSensitive() && fActive ) {
      volume->setModId(fModId);
      volume->SetModule(this);
      svList->Add(volume);
      aVol = dynamic_cast<FairGeoVolume*> ( node );
      fNodes->AddLast( aVol );
      fNbOfSensitiveVol++;
    }

  }
}
//__________________________________________________________________________
void  FairModule::AddSensitiveVolume(TGeoVolume* v)
{

  fLogger->Debug2(MESSAGE_ORIGIN, "FairModule::AddSensitiveVolume", v->GetName());
  // Only register volumes which are not already registered
  // Otherwise the stepping will be slowed down
  if( ! vList->findObject(v->GetName() ) ) {
    FairVolume*  volume = NULL;
    volume = new FairVolume(v->GetName(), fNbOfVolumes++);
    vList->addVolume(volume);
    volume->setModId(fModId);
    volume->SetModule(this);
    svList->Add(volume);
    fNbOfSensitiveVol++;
  }
}


//__________________________________________________________________________
FairVolume* FairModule::getFairVolume(FairGeoNode* fN)
{
  FairVolume* fv;
  FairVolume*  fvol=0;
  for(Int_t i=0; i<vList->getEntries(); i++) {
    fv=vList->At(i);
    if((fv->getGeoNode())==fN) {
      fvol=fv;
      return fvol;
      break;
    }
  }
  return fvol;
}
//__________________________________________________________________________
void FairModule::ConstructRootGeometry()
{
  TGeoManager* OldGeo=gGeoManager;
  TGeoManager* NewGeo=0;
  TGeoVolume* volume=0;;
  TFile* f=new TFile(GetGeometryFileName().Data());
  TList* l= f->GetListOfKeys();
  TKey* key;
  TIter next( l);
  TGeoNode* n=0;
  TGeoVolume* v1=0;
  while ((key = (TKey*)next())) {
    //std::cout << key->GetClassName() << std::endl;
    if (strcmp(key->GetClassName(),"TGeoManager") != 0) { continue; }
    gGeoManager=0;
    NewGeo = (TGeoManager*)key->ReadObj();
    break;
  }
  if (NewGeo!=0) {
    NewGeo->cd();
    volume=(TGeoVolume*)NewGeo->GetNode(0)->GetDaughter(0)->GetVolume();
    v1=volume->MakeCopyVolume(volume->GetShape());
    // n=NewGeo->GetTopNode();
    n=v1->GetNode(0);
    //  NewGeo=0;
    // delete NewGeo; //move it to the end of the method

  } else {
    //key=(TKey*) l->At(0);  //Get the first key in the list
    //volume=dynamic_cast<TGeoVolume*> (key->ReadObj());
    //if(volume!=0) { n=volume->GetNode(0); }
    //if(n!=0) { v1=n->GetVolume(); std::cout << n->GetName() << " " << GetName() << std::endl; }

    //
    //   - AYK, 2013/09/29 -> consider to perform a bit cleaner search (loop 
    //   through all keys and request name match); this way may use MC files
    //   to import detector geometry;
    //
    TKey* qkey;
    TIter qnext( l);
    
    while ((qkey = (TKey*)qnext())) {
      volume=dynamic_cast<TGeoVolume*> (qkey->ReadObj());
      if (!volume) continue;

      n = volume->GetNode(0);
      if (!n) continue;

      // Well, may want to do a better check later (remove "_N" suffix and request
      // exact match); later;
      if (TString(n->GetName()).BeginsWith(GetName()))
      {
  v1 = n->GetVolume();

  std::cout << n->GetName() << " " << GetName() << std::endl; 
  break;
      } //if
    }
  }

  if(v1==0) {
    fLogger->Fatal(MESSAGE_ORIGIN, "\033[5m\033[31mFairModule::ConstructRootGeometry(): could not find any geometry in File!!  \033[0m", GetGeometryFileName().Data());
  }
  gGeoManager=OldGeo;
  gGeoManager->cd();
  // If AddToVolume is empty add the volume to the top volume Cave
  // If it is defined check i´f the volume exists and if it exists add the volume from the root file
  // to the already existing volume
  TGeoVolume* Cave=NULL;
  if ( 0 == fMotherVolumeName.Length() ) {
    Cave= gGeoManager->GetTopVolume(); //printf("%s\n", Cave->GetName()); exit(0);
  } else {
    Cave = gGeoManager->GetVolume(fMotherVolumeName); //exit(0);
  }
  if(Cave!=NULL) {
    gGeoManager->AddVolume(v1); 
    TGeoVoxelFinder* voxels = v1->GetVoxels();
    if (voxels) { voxels->SetNeedRebuild(); }

    // else { fLogger->Fatal(MESSAGE_ORIGIN, "\033[5m\033[31mFairModule::ConstructRootGeometry(): could not find voxels  \033[0m"); }

    TGeoMatrix* M = n->GetMatrix();
    SetDefaultMatrixName(M);

    gGeoManager->GetListOfMatrices()->Remove(M);
    TGeoHMatrix* global = gGeoManager->GetHMatrix();
    gGeoManager->GetListOfMatrices()->Remove(global); //Remove the Identity matrix
    Cave->AddNode(v1,0, M);
    AssignMediumAtImport(v1);
    ExpandNode(n);
    if(NewGeo!=0) { delete NewGeo; }
    delete f;
  } else {
    fLogger->Fatal(MESSAGE_ORIGIN,"\033[5m\033[31mFairModule::ConstructRootGeometry(): could not find the given mother volume \033[0m   %s \033[5m\033[31m where the geomanger should be added. \033[0m", fMotherVolumeName.Data());
  }
}
//__________________________________________________________________________
void FairModule::ConstructASCIIGeometry()
{
  fLogger->Warning(MESSAGE_ORIGIN,"FairModule::ConstructASCIIGeometry() : this method has to be implimented in detector class");

}
//__________________________________________________________________________
Bool_t FairModule::CheckIfSensitive(std::string name)
{
  fLogger->Warning(MESSAGE_ORIGIN,"FairModule::CheckIfSensitive(std::string name): this method has to be implimented in detector class");
  return kFALSE;
}
//__________________________________________________________________________

void FairModule::ExpandNode(TGeoNode* fN)
{
  TGeoMatrix* Matrix =fN->GetMatrix();
  if(gGeoManager->GetListOfMatrices()->FindObject(Matrix)) { gGeoManager->GetListOfMatrices()->Remove(Matrix); }
  TGeoVolume* v1=fN->GetVolume();
  TObjArray* NodeList=v1->GetNodes();
  for (Int_t Nod=0; Nod<NodeList->GetEntriesFast(); Nod++) {
    TGeoNode* fNode =(TGeoNode*)NodeList->At(Nod);
    TGeoMatrix* M =fNode->GetMatrix();
    SetDefaultMatrixName(M);

    TGeoVolume* v= fNode->GetVolume();

    // Skip th erest of the loop (including ExpandNode() recursion) once a volume was served once; 
    // NB: this implementation is still not clean since v->RegisterYourself() has a recursion inside
    // (so daughters will be registered there more than one time, depending on level difference
    // in the hierarchy); ideally one would call n->GetVolume()->RegisterYourself() in 
    // ConstructRootGeometry() outside of the ExpandNode() call; this way it works, but 
    // simulation itself (for FEMC) becomes few times slower (why?); the implemented trivial 
    // solution speeds up FEMC init phase from ~2 minutes to <1 second and does not affect 
    // simulation time; leave it like this and be happy; eventually may want to arrange
    // a similar lookup table in AssignMediumAtImport() - see eg eic.patch file in EicRoot 
    // r219 where this stuff is commented out but not removed (media lookup is not really
    // needed, volume lookup suffices there);
    if (fGeoVolumeLut.find(v) != fGeoVolumeLut.end()) continue;
    fGeoVolumeLut.insert(v);

    if(fNode->GetNdaughters()>0) { ExpandNode(fNode); }
    AssignMediumAtImport(v);
    if (!gGeoManager->FindVolumeFast(v->GetName())) {
      fLogger->Debug2(MESSAGE_ORIGIN,"Register Volume : %s ", v->GetName());
      v->RegisterYourself();
    }
    if (CheckIfSensitive(v->GetName())) {
      fLogger->Debug2(MESSAGE_ORIGIN,"Sensitive Volume : %s ", v->GetName());
      AddSensitiveVolume(v);
    }
  }
}
//__________________________________________________________________________
void FairModule::SetDefaultMatrixName(TGeoMatrix* matrix)
{
  // Copied from root TGeoMatrix::SetDefaultName() and modified (memory leak)
  // If no name was supplied in the ctor, the type of transformation is checked.
  // A letter will be prepended to the name :
  //   t - translation
  //   r - rotation
  //   s - scale
  //   c - combi (translation + rotation)
  //   g - general (tr+rot+scale)
  // The index of the transformation in gGeoManager list of transformations will
  // be appended.
  if (!gGeoManager) { return; }
  if (strlen(matrix->GetName())) { return; }
  char type = 'n';
  if (matrix->IsTranslation()) { type = 't'; }
  if (matrix->IsRotation()) { type = 'r'; }
  if (matrix->IsScale()) { type = 's'; }
  if (matrix->IsCombi()) { type = 'c'; }
  if (matrix->IsGeneral()) { type = 'g'; }
  TObjArray* matrices = gGeoManager->GetListOfMatrices();
  Int_t index = 0;
  if (matrices) { index =matrices->GetEntriesFast() - 1; }
  matrix->SetName(Form("%c%i", type, index));
}

//__________________________________________________________________________

void FairModule::AssignMediumAtImport(TGeoVolume* v)
{
  FairGeoMedia* Media       = FairGeoLoader::Instance()->getGeoInterface()->getMedia();
  FairGeoBuilder* geobuild  = FairGeoLoader::Instance()->getGeoBuilder();

  TGeoMedium* med1=v->GetMedium();
  if(med1 && strcmp(med1->GetName(), "dummy")) {
    TGeoMaterial* mat1=v->GetMaterial();
    // FairGeoBuilder does not distinguish material and media names (I guess); STAR
    // geometry database does (say media CAVE_AIR is done from material AIR); so 
    // switch to using media names here and below;
#if _OLD_
    TGeoMaterial* newMat = gGeoManager->GetMaterial(mat1->GetName());
#else
    TGeoMaterial* newMat = gGeoManager->GetMaterial(med1->GetName());
#endif
    if( newMat==0) {
#if _OLD_
      FairGeoMedium* FairMedium=Media->getMedium(mat1->GetName());
#else
      FairGeoMedium* FairMedium=Media->getMedium(med1->GetName());
#endif
      if (!FairMedium) {
  if (fRootMaterialImportFlag) {
    // Create medium a la FairRoot by hand; sort of mimic FairGeoMedium::read() call 
    // sequence; this may cause conflicts if more than one ROOT file have material 
    // with the same name; FIXME: if this ever becomes a problem; for now I just need 
    // to import STAR geometry database;
    FairMedium = new FairGeoMedium(med1->GetName());
    
    FairMedium->setNComponents(mat1->GetNelements());
    //printf("%3d\n", mat1->GetNelements());
    for(unsigned iq=0; iq<mat1->GetNelements(); iq++) {
      double a, z, w;
      
      mat1->GetElementProp(a, z, w, iq);
      FairMedium->setComponent(iq, a, z, w);
    } //for iq
    
    FairMedium->setDensity(mat1->GetDensity());
    
    // Well, FairGeoMedium::read() also calculates rad.length for pure elements only;
    if (mat1->GetNelements() == 1) FairMedium->calcRadiationLength();
    
    // Parameter order is *ucked up compared to TGeo and parameters renamed -> 
    // carelessly use FairGeoMedium::setMediumPar() sequence; 
    double params[8];
    for(unsigned iq=0; iq<8; iq++)
      params[iq] = med1->GetParam(iq);
    FairMedium->sensFlag = params[0];//=sensFlag;
    FairMedium->fldFlag  = params[1];//=fldFlag;
    FairMedium->fld      = params[2];//=fld;
    FairMedium->madfld   = params[3];//=madfld;
    FairMedium->maxstep  = params[4];//=maxstep;
    FairMedium->maxde    = params[5];//=maxde;
    FairMedium->epsil    = params[6];//=epsil;
    FairMedium->minstep  = params[7];//=minstep;
    
    // And do not care about Cerenkov, sorry;
  }
  else
    fLogger->Fatal(MESSAGE_ORIGIN,"Material %s is not defined in ASCII file nor in Root file we Stop creating geometry", 
       mat1->GetName());
      } //if

      Int_t nmed=geobuild->createMedium(FairMedium);
      //FairMedium->print();
      v->SetMedium(gGeoManager->GetMedium(nmed));
      gGeoManager->SetAllIndex();
    } else {
#if _OLD_
      TGeoMedium* med2= gGeoManager->GetMedium(mat1->GetName());
#else
      TGeoMedium* med2= gGeoManager->GetMedium(med1->GetName());
#endif
      v->SetMedium(med2);
    }
  } else {
    if (strcmp(v->ClassName(),"TGeoVolumeAssembly") != 0) {
      //[R.K.-3.3.08]  // When there is NO material defined, set it to avoid conflicts in Geant
      fLogger->Fatal(MESSAGE_ORIGIN,"The volume  %s  Has no medium information and not an Assembly so we have to quit", v->GetName());
    }
  }
}

//__________________________________________________________________________
ClassImp(FairModule)