FairModule.cxx

Go to the documentation of this file.

/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
// -------------------------------------------------------------------------
// -----                     FairModule source file                    -----
// -----            Created 06/01/04  by M. Al-Turany                  -----
// -------------------------------------------------------------------------
/* Generated by Together */
#include "FairModule.h"

#include "FairGeoBuilder.h"     // for FairGeoBuilder
#include "FairGeoInterface.h"   // for FairGeoInterface
#include "FairGeoLoader.h"      // for FairGeoLoader
#include "FairGeoMedia.h"       // for FairGeoMedia
#include "FairGeoNode.h"        // for FairGeoNode
#include "FairGeoParSet.h"      // for FairBaseParSet
#include "FairMCApplication.h"
#include "FairRun.h"          // for FairRun
#include "FairRuntimeDb.h"    // for FairRuntimeDb
#include "FairVolume.h"       // for FairVolume
#include "FairVolumeList.h"   // for FairVolumeList

#include <TBuffer.h>           // for TBuffer, operator<<, etc
#include <TCollection.h>       // for TIter
#include <TFile.h>             // for TFile
#include <TGeoManager.h>       // for TGeoManager, gGeoManager
#include <TGeoMaterial.h>      // for TGeoMaterial
#include <TGeoMatrix.h>        // for TGeoMatrix, TGeoHMatrix
#include <TGeoMedium.h>        // for TGeoMedium
#include <TGeoNode.h>          // for TGeoNode
#include <TGeoVolume.h>        // for TGeoVolume
#include <TGeoVoxelFinder.h>   // for TGeoVoxelFinder
#include <TKey.h>              // for TKey
#include <TList.h>             // for TList, TListIter
#include <TObjArray.h>         // for TObjArray
#include <TObject.h>           // for TObject
#include <TRefArray.h>         // for TRefArray
#include <TSystem.h>           // for TSystem, gSystem
#include <TVirtualMC.h>

#ifdef ROOT_HAS_GDML
#include <TGDMLParse.h>
#endif

#include <map>
#include <stdlib.h>   // for getenv
#include <string.h>   // for strcmp, strlen

class FairGeoMedium;

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

void FairModule::ConstructGeometry()
{
    LOG(warn)
        << "The method ConstructGeometry has to be implemented in the detector class which inherits from FairModule";
}

void FairModule::ConstructOpGeometry()
{
    LOG(debug2)
        << "The method ConstructOpGeometry has to be implemented in the detector class which inherits from FairModule";
}

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(nullptr)
    , fGeoSaved(kFALSE)
    , fMC(nullptr)
{
    if (!svList) {
        svList = new TRefArray();
    }
    if (!vList) {
        vList = new FairVolumeList();
    }
}

FairModule::FairModule(const FairModule& rhs)
    : TNamed(rhs)
    , fMotherVolumeName(rhs.fMotherVolumeName)
    , fgeoVer(rhs.fgeoVer)
    , fgeoName(rhs.fgeoName)
    , fModId(rhs.fModId)
    , fActive(rhs.fActive)
    , fNbOfSensitiveVol(rhs.fNbOfSensitiveVol)
    , fVerboseLevel(rhs.fVerboseLevel)
    , flGeoPar(nullptr)
    , fGeoSaved(rhs.fGeoSaved)
    , fMC(nullptr)
{
    if (!svList) {
        svList = new TRefArray();
        for (Int_t i = 0; i < rhs.svList->GetEntries(); i++) {
            svList->Add(rhs.svList->At(i));
        }
    }

    if (!vList) {
        vList = new FairVolumeList();
        for (Int_t i = 0; i < rhs.vList->getEntries(); i++) {
            vList->addVolume(rhs.vList->At(i));
        }
    }

    // Initialize cached pointer to MC (on worker)
    fMC = TVirtualMC::GetMC();

    // TO DO - add when we know what type is the elements of flGeoPar
    // flGeoPar=new TObjArray();
    // TIter it = rhs.flGeoPar->MakeIterator();
    // Copy parameters
    // TObject* obj;
    // while((obj=it->Next())) {
    //  flGeoPar->Add(...);
    //}
}

FairModule::FairModule()
    : TNamed()
    , fMotherVolumeName("")
    , fgeoVer("Not defined")
    , fgeoName("Not defined")
    , fModId(-1)
    , fActive(kFALSE)
    , fNbOfSensitiveVol(0)
    , fVerboseLevel(0)
    , flGeoPar(nullptr)
    , fGeoSaved(kFALSE)
    , fMC(nullptr)
{}

FairModule& FairModule::operator=(const FairModule& rhs)
{
    // check assignment to self
    if (this == &rhs)
        return *this;

    // base class assignment
    TNamed::operator=(rhs);

    // assignment operator
    fMotherVolumeName = rhs.fMotherVolumeName;
    fgeoVer = rhs.fgeoVer;
    fgeoName = rhs.fgeoName;
    fModId = rhs.fModId;
    fActive = rhs.fActive;
    fNbOfSensitiveVol = rhs.fNbOfSensitiveVol;
    fVerboseLevel = rhs.fVerboseLevel;
    flGeoPar = nullptr;
    fGeoSaved = rhs.fGeoSaved;
    fMC = TVirtualMC::GetMC();

    // TO DO - add when we know what type is the elements of flGeoPar
    // flGeoPar=new TObjArray();
    // TIter it = rhs.flGeoPar->MakeIterator();
    // copy parameters
    // TObject* obj;
    // while((obj=it->Next())) {
    //  flGeoPar->Add(...);
    //}

    return *this;
}

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)
{

    // If absolute path is given as argument, try to find it there.
    // If the file don't exist break. Don't look anywhere else.
    if (fname.BeginsWith("/")) {
        if (gSystem->AccessPathName(fname.Data())) {
            LOG(fatal) << fName << ": geometry file " << fname << " not found in absolut path!";
            fgeoName = "";
        }   // file not found
        fgeoName = fname;
        LOG(debug) << fName << ": using geometry file " << fgeoName;
        return;
    }

    // If GEOMPATH is set, try to find the file there.
    // If there is no such file go on to look in the default location.
    TString userPath = getenv("GEOMPATH");
    userPath.ReplaceAll("//", "/");
    if (!userPath.IsNull()) {
        if (!userPath.EndsWith("/")) {
            userPath += "/";
        }
        fgeoName = userPath + fname;
        if (!gSystem->AccessPathName(fgeoName.Data())) {
            LOG(debug) << fName << ": using geometry file " << fgeoName;
            return;
        }
        LOG(debug) << fName << ": geometry file " << fname << " not found in GEOMPATH " << userPath;
    }

    // Look in the standard path
    fgeoName = getenv("VMCWORKDIR");
    fgeoName += "/geometry/" + fname;
    if (!gSystem->AccessPathName(fgeoName.Data())) {
        LOG(debug) << fName << ": using geometry file " << fgeoName;
        return;
    }

    // File not found
    LOG(fatal) << fName << ": geometry file " << fname << " not found in standard path ";
    fgeoName = "";
    return;
}

void FairModule::ProcessNodes(TList* aList)
{
    if (FairMCApplicationState::kConstructGeometry != FairMCApplication::Instance()->GetState()) {
        LOG(fatal) << "Detected call to FairModule::ProcessNodes() \
      while not in FairMCApplication::ConstructGeometry()\n\
      Call templated function FairModule::ConstructASCIIGeometry()\
      from ConstructGeometry() of your detector class. Aborting...";
    }

    if (!svList) {
        svList = new TRefArray();
    }
    if (!vList) {
        vList = new FairVolumeList();
    }

    TListIter iter(aList);
    FairGeoNode* node = nullptr;
    FairGeoNode* MotherNode = nullptr;
    FairVolume* volume = nullptr;
    FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
    FairGeoParSet* par = static_cast<FairGeoParSet*>(rtdb->getContainer("FairGeoParSet"));
    TObjArray* fNodes = par->GetGeoNodes();
    while ((node = static_cast<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 = nullptr;

        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)
{

    LOG(debug2) << "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 = nullptr;
        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(TGeoMatrix* shiftM)
{
    TGeoManager* OldGeo = gGeoManager;
    TGeoManager* NewGeo = 0;
    TGeoVolume* volume = 0;
    ;
    TFile* f = TFile::Open(GetGeometryFileName().Data());
    TList* l = f->GetListOfKeys();
    TKey* key;
    TIter next(l);
    TGeoNode* n = 0;
    TGeoVolume* v1 = 0;
    while ((key = static_cast<TKey*>(next()))) {
        if (strcmp(key->GetClassName(), "TGeoManager") != 0) {
            continue;
        }
        gGeoManager = 0;
        NewGeo = static_cast<TGeoManager*>(key->ReadObj());
        break;
    }
    if (NewGeo != 0) {
        NewGeo->cd();
        volume = static_cast<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 = static_cast<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();
        }
    }

    if (v1 == 0) {
        LOG(fatal) << "Could not find any geometry in file " << GetGeometryFileName().Data();
    } else {
        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 = nullptr;
        if (0 == fMotherVolumeName.Length()) {
            Cave = gGeoManager->GetTopVolume();
        } else {
            Cave = gGeoManager->GetVolume(fMotherVolumeName);
        }
        if (Cave != nullptr) {
            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();

            // very nasty!
            TGeoHMatrix* M2 = new TGeoHMatrix(*M);

            // For debugging
            // M2->Dump();
            // const Double_t* transl1 = M2->GetTranslation();
            // std::cout << transl1[0] << "\t" << transl1[1] << "\t" << transl1[2] << std::endl;
            // if (shiftM) {
            //  shiftM->Dump();
            //  const Double_t* transl2 = shiftM->GetTranslation();
            //  std::cout << transl2[0] << "\t" << transl2[1] << "\t" << transl2[2] << std::endl;
            //}

            if (shiftM) {
                M2->Multiply(shiftM);   // HACK!
            }

            SetDefaultMatrixName(M2);

            // TODO
            // I don't really understand this juggling with the matices,
            // so, please, take care of how it works with my changes.
            // Egor.

            gGeoManager->GetListOfMatrices()->Remove(M2);
            TGeoHMatrix* global = gGeoManager->GetHMatrix();
            gGeoManager->GetListOfMatrices()->Remove(global);   // Remove the Identity matrix
            Cave->AddNode(v1, 0, M2);
            AssignMediumAtImport(v1);
            ExpandNode(n);
            delete NewGeo;
            delete f;
        } else {
            LOG(fatal) << "Could not find the given mother volume " << fMotherVolumeName.Data()
                       << " where the geomanger should be added.";
        }
    }
}

#ifdef ROOT_HAS_GDML

void FairModule::ConstructGDMLGeometry(__attribute__((unused)) TGeoMatrix* posrot)
{
    // Parse the GDML file
    TFile* old = gFile;
    TGDMLParse parser;
    TGeoVolume* gdmlTop;
    gdmlTop = parser.GDMLReadFile(GetGeometryFileName());

    // Change ID of media. TGDMLParse starts allways from 0. Need to shift.
    ReAssignMediaId();

    // Add volume to the cave and go through it recursively
    gGeoManager->GetTopVolume()->AddNode(gdmlTop, 1, posrot);
    ExpandNodeForGDML(gGeoManager->GetTopVolume()->GetNode(gGeoManager->GetTopVolume()->GetNdaughters() - 1));
    gFile = old;
}

void FairModule::ExpandNodeForGDML(TGeoNode* curNode)
{
    // Get pointer to volume and assign medium
    TGeoVolume* curVol = curNode->GetVolume();
    AssignMediumAtImport(curVol);

    // Check if the volume is sensitive
    if ((this->InheritsFrom("FairDetector")) && IsSensitive(curVol->GetName())) {
        LOG(debug2) << "Sensitive Volume " << curVol->GetName();
        AddSensitiveVolume(curVol);
    }

    if (curVol->GetNdaughters() != 0) {
        TObjArray* NodeChildList = curVol->GetNodes();
        TGeoNode* curNodeChild;
        for (Int_t j = 0; j < NodeChildList->GetEntriesFast(); j++) {
            curNodeChild = static_cast<TGeoNode*>(NodeChildList->At(j));
            ExpandNodeForGDML(curNodeChild);
        }
    }
}

#else

void FairModule::ConstructGDMLGeometry(__attribute__((unused)) TGeoMatrix* posrot)
{
    LOG(error) << "Could not construct magnet geometry from gdml file.";
    LOG(error) << "The used ROOT version does not support gdml.";
    LOG(error) << "Please recompile ROOT with gdml support.";
    LOG(fatal) << "Stop execution at this point.";
}

void FairModule::ExpandNodeForGDML(__attribute__((unused)) TGeoNode* curNode) {}

#endif

void FairModule::ReAssignMediaId()
{
    // Initialise pointer to GeoBuilder
    FairGeoBuilder* geoBuilder = FairGeoLoader::Instance()->getGeoBuilder();
    // Get list of TGeo media
    TList* media = gGeoManager->GetListOfMedia();
    // Loop over new media which are not in GeoBase and shift the ID
    TGeoMedium* med;
    //    TGeoMedium* med2;
    for (Int_t i = geoBuilder->GetNMedia(); i < media->GetEntries(); i++) {
        med = static_cast<TGeoMedium*>(media->At(i));
        med->SetId(i + 1);
    }
    // Change GeoBase medium index
    geoBuilder->SetNMedia(media->GetEntries());

    // Revove dublicated materials
    TList* materials = gGeoManager->GetListOfMaterials();
    TIter next1(materials);
    // map for existing materials
    std::map<TString, Bool_t> mapMatName;
    TGeoMaterial* mat;
    while ((mat = static_cast<TGeoMaterial*>(next1()))) {
        // If material exist - delete dublicated. If not - set the flag
        if (mapMatName[mat->GetName()]) {
            materials->Remove(mat);
        } else {
            mapMatName[mat->GetName()] = kTRUE;
        }
    }
}

void FairModule::ConstructASCIIGeometry()
{
    LOG(warn) << "The method ConstructASCIIGeometry has to be implemented in the detector class which inherits from "
                 "FairModule";
}

//__________________________________________________________________________
Bool_t FairModule::IsSensitive(const std::string& name)
{
    //    LOG(warn) << "Implement IsSensitive in the detector class which inherits from FairModule";
    //    LOG(warn) << "For now calling the obsolete function CheckIfSensitive";
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    return CheckIfSensitive(name);
#pragma GCC diagnostic pop
}

void FairModule::ExpandNode(TGeoNode* fN)
{
    // FairGeoLoader* geoLoad = FairGeoLoader::Instance();
    // FairGeoInterface* geoFace = geoLoad->getGeoInterface();
    // FairGeoMedia* Media =  geoFace->getMedia();
    // FairGeoBuilder* geobuild=geoLoad->getGeoBuilder();
    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 = static_cast<TGeoNode*>(NodeList->At(Nod));
        TGeoMatrix* M = fNode->GetMatrix();
        // M->SetDefaultName();
        SetDefaultMatrixName(M);
        if (fNode->GetNdaughters() > 0) {
            ExpandNode(fNode);
        }
        TGeoVolume* v = fNode->GetVolume();
        AssignMediumAtImport(v);
        if (!gGeoManager->FindVolumeFast(v->GetName())) {
            LOG(debug2) << "Register Volume " << v->GetName();
            v->RegisterYourself();
        }
        if ((this->InheritsFrom("FairDetector")) && IsSensitive(v->GetName())) {
            LOG(debug2) << "Sensitive Volume " << v->GetName();
            AddSensitiveVolume(v);
        }
        v->GetShape()->AfterStreamer();
    }
}

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) {
        // In newer ROOT version also a TGeoVolumeAssembly has a material and medium.
        // This medium is called dummy and is automatically set when the geometry is constructed.
        // Since this material and medium is neither in the TGeoManager (at this point) nor in our
        // ASCII file we have to create it the same way it is done in TGeoVolume::CreateDummyMedium()
        // In the end the new medium and material has to be added to the TGeomanager, because this is
        // not done automatically when using the default constructor. For all other constructors the
        // newly created medium or material is added to the TGeomanger.
        // Create the medium and material only the first time.
        TString medName = static_cast<TString>(med1->GetName());
        if (medName.EqualTo("dummy") && nullptr == gGeoManager->GetMedium(medName)) {

            TGeoMaterial* dummyMaterial = new TGeoMaterial();
            dummyMaterial->SetName("dummy");

            TGeoMedium* dummyMedium = new TGeoMedium();
            dummyMedium->SetName("dummy");
            dummyMedium->SetMaterial(dummyMaterial);

            gGeoManager->GetListOfMedia()->Add(dummyMedium);
            gGeoManager->AddMaterial(dummyMaterial);
        }

        TGeoMaterial* mat1 = v->GetMaterial();
        TGeoMaterial* newMat = gGeoManager->GetMaterial(mat1->GetName());
        if (newMat == 0) {
            FairGeoMedium* FairMedium = Media->getMedium(mat1->GetName());
            if (!FairMedium) {
                LOG(fatal) << "Material " << mat1->GetName() << "is not defined in ASCII file nor in Root file.";
                //     FairMedium=new FairGeoMedium(mat1->GetName());
                //      Media->addMedium(FairMedium);
            } else {

                Int_t nmed = geobuild->createMedium(FairMedium);
                v->SetMedium(gGeoManager->GetMedium(nmed));
                gGeoManager->SetAllIndex();
            }
        } else {
            TGeoMedium* med2 = gGeoManager->GetMedium(mat1->GetName());
            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
            LOG(fatal) << "The volume " << v->GetName()
                       << "has no medium information and not an Assembly so we have to quit";
        }
    }
}

FairModule* FairModule::CloneModule() const
{
    Fatal("CloneModule", "Has to be overriden in multi-threading applications.");
    return 0;
}

ClassImp(FairModule);