MuonCounterViewer.cc

Go to the documentation of this file.

#include "MuonCounterViewer.h"

#include <evt/Event.h>
#include <det/Detector.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mdet/Module.h>
#include <mdet/Scintillator.h>
#include <sdet/Station.h>
#include <utl/TimeStamp.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ErrorLogger.h>
#include <utl/TransformationMatrix.h>
#include <utl/Transformation.h>
#include <utl/MathConstants.h>
#include <utl/Line.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>

#include <sstream>
#include <memory>
#include <ctime>
#include <fstream>

#include <TGeoManager.h>
#include <TGeoMaterial.h>
#include <TGeoMedium.h> 
#include <TGeoMatrix.h>
#include <TGeoVolume.h> 
#include <TVirtualGeoPainter.h>
 
using namespace std;
using namespace utl;
using namespace fwk;


namespace MuonCounterViewerAG {


  void
  ScintBuilder::Box(const mdet::Scintillator& s, double w, double l, double h)
  {
    // A new scint name.
    std::ostringstream iStrScint;
    iStrScint << "scint" 
              << s.GetId() << "_" 
              << s.GetModule().GetId() << "_" 
              << s.GetModule().GetCounter().GetId();
    // Keep the reference.
    fVolume = fManager->MakeBox(iStrScint.str().c_str(), fMedium, w / utl::m / 2 , l / utl::m / 2, h / utl::m / 2);
  }


  VModule::ResultFlag
  MuonCounterViewer::Init()
  {
    det::Detector::GetInstance().Update(GetCurrentSystemTime());
    const Branch topB = CentralConfig::GetInstance()->GetTopBranch("MuonCounterViewer");
    topB.GetChild("CounterId").GetData(counterId);
    return eSuccess;
  }

  VModule::ResultFlag
  MuonCounterViewer::Run(evt::Event& /*event*/)
  {
    const mdet::MDetector& d = det::Detector::GetInstance().GetMDetector();
    typedef mdet::MDetector::CounterConstIterator Itc;
    typedef mdet::Counter::ModuleConstIterator Itm;
    typedef mdet::Module::ScintillatorConstIterator Its;
    auto_ptr<TGeoManager> geom(new TGeoManager("counter", "A muon counter."));
    TGeoMaterial* matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    TGeoMedium* vacuum      = new TGeoMedium("Vacuum",1, matVacuum);
    // For constructing Scintillator volumes
    ScintBuilder builder;
    // XXX Consider to represent the hierarchy also here...
    // Make a box for the top:
    TGeoVolume *top = geom->MakeBox("TOP", vacuum, 1, 1, 1);
    geom->SetTopVolume(top);

    for (Itc ic = d.AllCountersBegin(); ic != d.AllCountersEnd(); ++ic) {
      const mdet::Counter& c = *ic;
      
      if (ic->GetId() != counterId)
        continue;
    
      // Ground.
      const double groundHeight = 1 * utl::cm;
      TGeoVolume* ground = geom->MakeBox("ground", vacuum, 10 * utl::m , 10 * utl::m, groundHeight / 2);
      TGeoTranslation* transGround = new TGeoTranslation(0,0, 
        c.GetAssociatedTank().GetPosition().GetZ(c.GetLocalCoordinateSystem()) / utl::m);
      ground->SetLineColor(kOrange-1);
      top->AddNode(ground, 1, transGround);
      // Build tags.
      std::ostringstream iStrCounter;
      iStrCounter<< "counter" << c.GetId();
      TGeoVolume* cCenter = geom->MakeSphere(iStrCounter.str().c_str(),vacuum, 10*utl::cm,12*utl::cm);
      // No need to rotate, as where putting an sphere.
      TGeoTranslation* transCounter = new TGeoTranslation(c.GetPosition().GetX(c.GetLocalCoordinateSystem()) / utl::m, 
                                                          c.GetPosition().GetY(c.GetLocalCoordinateSystem()) / utl::m, 
                                                          c.GetPosition().GetZ(c.GetLocalCoordinateSystem()) / utl::m);
      cCenter->SetLineColor(kBlack);
      top->AddNode(cCenter, 1, transCounter);
      // Tank.
      TGeoVolume* tank = geom->MakeTube("tank", vacuum, 0, 
        c.GetAssociatedTank().GetRadius() / utl::m, 
        c.GetAssociatedTank().GetHeight() / utl::m / 2);
      TGeoTranslation* transTank = new TGeoTranslation(
        c.GetAssociatedTank().GetPosition().GetX(c.GetLocalCoordinateSystem()) / utl::m,
        c.GetAssociatedTank().GetPosition().GetY(c.GetLocalCoordinateSystem()) / utl::m,        
        c.GetAssociatedTank().GetHeight() /  utl::m / 2  + 
          c.GetAssociatedTank().GetPosition().GetZ(c.GetLocalCoordinateSystem()) / utl::m);
      tank->SetLineColor(41);
      top->AddNode(tank, 1, transTank);
      for (Itm im = c.ModulesBegin(); im != c.ModulesEnd(); ++im) {
        const mdet::Module& m = *im;
        std::ostringstream iStrModule;
        iStrModule << "module" << m.GetId() << "_" << c.GetId();
        TGeoVolume* mCenter = geom->MakeSphere(iStrModule.str().c_str(),vacuum, 10*utl::cm,12*utl::cm);
        // No need to rotate, as we're putting an sphere.
        TGeoTranslation* transModule = new TGeoTranslation(m.GetPosition().GetX(c.GetLocalCoordinateSystem()) / utl::m, 
                                                           m.GetPosition().GetY(c.GetLocalCoordinateSystem()) / utl::m, 
                                                           m.GetPosition().GetZ(c.GetLocalCoordinateSystem()) / utl::m);
        mCenter->SetLineColor(kViolet);
        top->AddNode(mCenter, 1, transModule);
        for (Its is = m.ScintillatorsBegin(); is != m.ScintillatorsEnd(); ++is) {
          const mdet::Scintillator& s = *is;
          // Scintillator.
          utl::Point sCenter(0, 0, 0, s.GetLocalCoordinateSystem());
          double xxx = sCenter.GetX(c.GetLocalCoordinateSystem());
          double yyy = sCenter.GetY(c.GetLocalCoordinateSystem());
          double zzz = sCenter.GetZ(c.GetLocalCoordinateSystem());
          // Find the rotation from scin to counter reference frame.
          utl::TransformationMatrix scintRot = s.GetLocalCoordinateSystem()
            ->GetTransformationTo(c.GetLocalCoordinateSystem()).Decompose().get<1>();
          // Construct the intermediate object to perform "active" transformation.
          utl::TransformationPolicy scintPol(scintRot,s.GetLocalCoordinateSystem());
          // Since there's no matrix project we build a local basis...
          utl::Point xUnit(1, 0, 0, s.GetLocalCoordinateSystem());
          utl::Point yUnit(0, 1, 0, s.GetLocalCoordinateSystem());
          utl::Point zUnit(0, 0, 1, s.GetLocalCoordinateSystem());
          // ... rotate the vectors to the counter system.
          xUnit = scintPol * xUnit;
          yUnit = scintPol * yUnit;
          zUnit = scintPol * zUnit;
          // Plain'ol vector for Root.
          double rotMatrix[9];
          // So that by querying the resulting compoenents (in the same system)
          // we can retrieve the matrix.
          rotMatrix [ 0 ] = xUnit.GetX(s.GetLocalCoordinateSystem());
          rotMatrix [ 1 ] = yUnit.GetX(s.GetLocalCoordinateSystem());
          rotMatrix [ 2 ] = zUnit.GetX(s.GetLocalCoordinateSystem());
          rotMatrix [ 3 ] = xUnit.GetY(s.GetLocalCoordinateSystem());
          rotMatrix [ 4 ] = yUnit.GetY(s.GetLocalCoordinateSystem());
          rotMatrix [ 5 ] = zUnit.GetY(s.GetLocalCoordinateSystem());
          rotMatrix [ 6 ] = xUnit.GetZ(s.GetLocalCoordinateSystem());
          rotMatrix [ 7 ] = yUnit.GetZ(s.GetLocalCoordinateSystem());
          rotMatrix [ 8 ] = zUnit.GetZ(s.GetLocalCoordinateSystem());
          // So build the "rootation".
          TGeoRotation* rot = new TGeoRotation();
          rot->SetMatrix(rotMatrix);
          /*
           * "The rotation won't be owned" tells the doxy from the constructor of 
           * TGeoCombiTrans: so register it.
           * In other parts of the docs (see TGeoMatrix): "unless explicitly used 
           * for positioning nodes (TGeoVolume::AddNode()) all
           * matrices deletion have to be managed by users. Matrices passed to geometry
           * have to be created by using new() operator and their deletion is done by
           * TGeoManager class".
           *
           * The use of an stack variable for rot was used, but it led to a seg-fault:
           * it isn't owned but may be pointed: this seems to be the case, see
           * void TGeoCombiTrans::SetRotation(const TGeoRotation *rot)
           * which is the method called by the constructor.
           */
          rot->RegisterYourself();
          // Compose the translation with the rotation.
          // Build the box with its dimensions.
          TGeoVolume* v = s.VisitShape(builder(geom.get(), vacuum)).GetVolume();
          // Some colors for the Scint, due to its module tag.
          switch(m.GetAreaKind()) {
          case mdet::Module::eLarge:
            v->SetLineColor(kRed);
            break;
          case mdet::Module::eSmall:
            v->SetLineColor(kBlue);
          }
          top->AddNode(v, 1, new TGeoCombiTrans(xxx / utl::m, yyy / utl::m, zzz / utl::m, rot));
        }
      }
      //Only one counter, no more.
      //break;
    }

    /*
     * With this line the "triggering" initial Draw is added automatically
     * to the script.
     */ 
    geom->GetGeomPainter()->SetTopVolume(top);
    /*
     * Give the script as output.
     * It seems that it isn't neccesary to
     * "close" the geometry.
     * See the following commented call to this function for more info.
     *
     * Enclose the dumped code in a method. 
     * I just tried in root-5.24 and note that the output of SavePrimitive method
     * uses some variales without prior declaration. If the script is used directly,
     * that gives problems (just because of the lack of declaration) but if the script
     * is run inside a methos, it is just fine.
     */

     std::ostringstream iStrFile;
     iStrFile << "CounterViewer_" << counterId << ".C";
     ofstream fout(iStrFile.str().c_str());
     fout << "void CounterViewer_" << counterId << "() {" << endl;
     top->SavePrimitive(fout);
     fout << "}" << endl;

    /*
     * This call, in certain oportunities, gave a problem during Voxelizing: Wooops!!! appears in 
     * standard-output:
     * /lib/tls/i686/cmov/libc.so.6[0xb0ce53f4]
     * /usr/lib/libstdc++.so.6(_ZdlPv+0x21)[0xb0e98031]
     * /usr/lib/libstdc++.so.6(_ZdaPv+0x1d)[0xb0e9808d]
     * /usr/local/auger/External/ROOT/pro/lib/libGeom.so(_ZN15TGeoVoxelFinder7SortAllEPKc+0x6f1)[0xb0fede31]
     * /usr/local/auger/External/ROOT/pro/lib/libGeom.so(_ZN15TGeoVoxelFinder8VoxelizeEPKc+0x3a)[0xb0fe7b5a]
     * /usr/local/auger/External/ROOT/pro/lib/libGeom.so(_ZN10TGeoVolume8VoxelizeEPKc+0x188)[0xb0fdf5e8]
     * /usr/local/auger/External/ROOT/pro/lib/libGeom.so(_ZN11TGeoManager8VoxelizeEPKc+0xb1)[0xb0f849e1]
     * /usr/local/auger/External/ROOT/pro/lib/libGeom.so(_ZN11TGeoManager13CloseGeometryEPKc+0x1cc)[0xb0f84c1c]
     * ./userAugerOffline(_ZN19MuonCounterViewerAG17MuonCounterViewer3RunERN3evt5EventE+0x2315)[0x8059295]
     *
     * This problem couldn't be isolated and tested with Root 5.20.
     *
     * geom->CloseGeometry();
     */
    /*
     * This call works but afterwards the program gives seg-fault in root 5.16. Not in 5.20. 
     * Both of course working in my box with Ubuntu 7.10.
     *
     * #0  0x08be481c in ?? ()
     * #1  0xb2b627da in TCollection::GarbageCollect () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #2  0xb2b6aefb in TObjArray::Delete () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #3  0xb0f74832 in TGeoManager::~TGeoManager () from /usr/local/auger/External/ROOT/pro/lib/libGeom.so
     * #4  0xb2b627da in TCollection::GarbageCollect () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #5  0xb2b67d66 in TList::Delete () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #6  0xb2b20d23 in TROOT::~TROOT () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #7  0xb2b20e70 in __tcf_0 () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #8  0xb0c91101 in __cxa_finalize () from /lib/tls/i686/cmov/libc.so.6
     * #9  0xb2ae0693 in __do_global_dtors_aux () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #10 0xb2f11d7c in _fini () from /usr/local/auger/External/ROOT/pro/lib/libCore.so
     * #11 0xb7f82153 in ?? () from /lib/ld-linux.so.2
     * #12 0xb7f90260 in _rtld_global () from /lib/ld-linux.so.2
     * #13 0x00000000 in ?? ()
     *
     * top->SavePrimitive(cout);
     */
    return eSuccess;
  }


  VModule::ResultFlag
  MuonCounterViewer::Finish()
  {
    cout << "\n\n******************\n";
    INFO("ROOT macro script created");

    std::ostringstream iStrMessage;
    iStrMessage << "From command line, execute: root CounterViewer_" << counterId << ".C";
    INFO(iStrMessage.str().c_str());
    cout << "n******************\n\n";

    return eSuccess;
  }

}