JDrawDetector3D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <limits>

#include "TROOT.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TView.h"
#include "TGeometry.h"
#include "TGeoManager.h"
#include "TGeoMatrix.h"
#include "TGeoMaterial.h"
#include "TGeoMedium.h"
#include "TGeoVolume.h"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"


/**
 * \file
 * Auxiliary program to draw the detector in 3D.
 * \author mdejong
 */
int main(int argc, char**argv)
{
  using namespace std;

  string         detectorFile;
  int            debug;

  try {

    JParser<> zap("Auxiliary program to draw the detector in 3D.");

    zap['a'] = make_field(detectorFile);
    zap['d'] = make_field(debug)             = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPP;
  
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  if (!detector.empty()) {

    JPosition3D pos(0.0, 0.0, 0.0);

    for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
      pos += *i;
    }
    
    pos /= detector.size();

    NOTICE("Detector center: " << pos.getX() << ' ' << pos.getY() << ' ' << pos.getZ() << endl);

    double zmin = numeric_limits<double>::max();
    double zmax = numeric_limits<double>::lowest();

    double radius = 0.0;

    for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

      const double x = i->getX() - pos.getX();
      const double y = i->getY() - pos.getY();
      const double z = i->getZ() - pos.getZ();

      const double R = sqrt(x*x + y*y);

      if (z < zmin) zmin = z;
      if (z > zmax) zmax = z;

      if (R > radius) radius = R;
    }

    NOTICE("can: " << zmin << ' ' << zmax << ' ' << radius << endl);
  }


  double xmax = 0.0;
  double ymax = 0.0;
  double zmax = 0.0;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    if (i->getX() > xmax) xmax = i->getX();
    if (i->getY() > ymax) ymax = i->getY();
    if (i->getZ() > zmax) zmax = i->getZ();
  }

  xmax *= 2.5;
  ymax *= 2.5;
  zmax *= 2.5;


  TApplication* tapple = new TApplication("user", NULL, NULL);
  TCanvas*      canvas = new TCanvas("a", detectorFile.c_str(), 1200, 1200);

  canvas->SetFillColor(0);

  TGeoManager*  geom     = new TGeoManager ("geometry", "");
  TGeoMaterial* material = new TGeoMaterial("vacuum", 0, 0, 0);
  TGeoMedium*   medium   = new TGeoMedium  ("vacuum", 1, material);
  TGeoVolume*   top      = geom->MakeBox   ("Top",    medium, xmax, ymax, zmax);
  TGeoVolume*   shape    = geom->MakeSphere("Module", medium, 0.4, 0.5); 

  shape->SetLineColor(kBlue);

  int N = 0;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    top->AddNode(shape, N++, new TGeoTranslation(i->getX(), i->getY(), i->getZ()));
  }

  geom->SetTopVolume(top);
  geom->CloseGeometry();
  //geom->SetVisLevel(4);

  top->Draw();

  canvas->GetView()->ShowAxis();
  canvas->Update();
  tapple->Run();
}