vistpc.cxx

Go to the documentation of this file.

#include "vistpc.h"

// C++ STL

#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>

#include <cmath>

#include <map>

// ROOT

#include <TSystem.h>
#include <TColor.h>
#include <TList.h>
#include <TStyle.h>
#include <TRandom.h>
#include <TMath.h>

#include <TEveBoxSet.h>
#include <TEveBrowser.h>
#include <TEveCompound.h>
#include <TEveGeoNode.h>
#include <TEveManager.h>
#include <TEveRGBAPalette.h>
#include <TEvePointSet.h>
#include <TEveQuadSet.h>
#include <TEveStraightLineSet.h>
#include <TEveTrack.h>
#include <TEveTrackPropagator.h>
#include <TEveVSDStructs.h>

#include <TGeometry.h>
#include <TGeoManager.h>
#include <TGeoMaterial.h>
#include <TGeoMatrix.h>
#include <TGeoTube.h>
#include <TGeoXtru.h>

#include <TGListBox.h>
#include <TGListTree.h>
#include <TGMenu.h>

// GEAR

#include <gearimpl/Util.h>
#include <gearxml/GearXML.h>
#include <gear/GearMgr.h>
#include <gear/GEAR.h>
#include <gear/TPCParameters.h>

// LCIO

#include "lcio.h"

#include "EVENT/LCIO.h"
#include "IO/LCReader.h"
#include "IO/LCWriter.h"
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCCollectionVec.h"
#include "EVENT/Track.h"
#include "EVENT/TrackerData.h"
#include "EVENT/TrackerHit.h"
#include "EVENT/TrackerPulse.h"
#include "Exceptions.h"

using namespace lcio;
using namespace gear;

//==============================================================================
// Constructor / Destructor
//------------------------------------------------------------------------------


VisTPC::VisTPC(gear::GearMgr * aGearMgr) :
   debug(1),
   readoutFrequency(40.),
   driftVelocity(79.),
   maxPoints(0)
{
   // to avoid code duplication this code is in a helper function
   updateGearMgr( aGearMgr );

   TStyle* style = new TStyle();
   TGeoManager* geom = new TGeoManager("DetectorGeometry", "detector geometry");

   //--- define some materials
   materials["Vacuum"] = new TGeoMaterial("Vacuum", 0,0,0);
   materials["Al"]     = new TGeoMaterial("Al", 26.98,13,2.7);

   //--- define some media
   media["Vacuum"]     = new TGeoMedium("Vacuum",    1, materials["Vacuum"]);
   media["Al"]         = new TGeoMedium("Aluminium", 2, materials["Al"]);

   // default options
   options["module.parallelBorders"] = 0;
   options["module.drawPadsByDefault"] = kFALSE;
   options["module.drawOutlineByDefault"] = kTRUE; 
}

//==============================================================================
// Methods (Alphabetical order)
//------------------------------------------------------------------------------

void VisTPC::drawCircularPolygon(TEveStraightLineSet* ls, float r, int nDivisions, float z) {
   drawCircularPolygon(ls, r, 0, 2 * TMath::Pi(), nDivisions, z);
}

void VisTPC::drawCircularPolygon(TEveStraightLineSet* ls, float r, float phi_min, float phi_max, int nDivisions, float z) {

   // Checks of the arguments

   if (nDivisions == 0) {
      return;
   }

   if (r == 0) {
      return;
   }

   // Implementation 1
   /*
   float phi  = phi_min;
   float dPhi = (phi_max - phi_min) / nDivisions;

   while (phi + dPhi <= phi_max) {
      ls->AddLine(r*cos(phi),     r*sin(phi),     z,
                  r*cos(phi+dPhi), r*sin(phi+dPhi), z);

      phi += dPhi;
   }
   */

   // Implementation 2

   float dPhi = (phi_max - phi_min) / nDivisions;
   int   j    = 0;

   while (j < nDivisions) {
      ls->AddLine(r*cos(dPhi*j + phi_min),     r*sin(dPhi*j + phi_min),     z,
                  r*cos(dPhi*(j+1) + phi_min), r*sin(dPhi*(j+1) + phi_min), z);

      ++j;
   }
}

void VisTPC::drawEvent(EVENT::LCEvent* event, TList &selectedEntries, TEveElement* parent, const char* name) {
   std::vector<std::string> entries;

   TIter next(&selectedEntries);
   while (TGTextLBEntry* entry = (TGTextLBEntry*) next()) {
      entries.push_back(entry->GetTitle());
   }

   drawEvent(event, entries, parent, name);
}

void VisTPC::drawEvent(EVENT::LCEvent* event, std::vector<std::string> &selectedEntries, TEveElement* parent, const char* name) {
   TEveElementList *eveNode      = new TEveElementList(name);
   TEveElement     *eveChild     = NULL;
   TGListTreeItem  *listTreeItem = 0;

   LCCollectionVec *collection = 0;

   //

   eveNode->SetRnrSelfChildren(kFALSE, kFALSE);

   //

   for (
      std::vector<std::string>::iterator it = selectedEntries.begin();
      it != selectedEntries.end();
      ++it
   ) {

      try {

         //std::cout << "Trying collection \"" << *it
         //          << "\"." << std::endl;

         collection = dynamic_cast<LCCollectionVec*>(event->getCollection(*it));

         //std::cout << "Going to draw " << *it << " "
         //          << "of type "
         //          << collection->getTypeName()
         //          << "("
         //          << gBBQ->getCollectionType(collection->getTypeName())
         //          << ")."
         //          << std::endl;

      } catch (DataNotAvailableException &) {

         std::cout << "Collection \"" << *it
                   << "\" not available" << std::endl;

         continue;

      }

      switch (
         getCollectionType(collection->getTypeName())
      ) {
         case kUNDEFINEDTYPE:
            break;

         case kHIT:
            eveChild = drawHits(collection);
            eveChild->SetElementName(it->c_str());

            if (parent) {
               parent->AddElement(eveChild);
            } else {
               gEve->AddElement(eveChild);
            }

            break;

         case kPULSE:
            eveChild = drawPulses(collection);
            eveChild->SetElementName(it->c_str());
            eveChild->SetRnrSelf(kFALSE);

            if (parent) {
               parent->AddElement(eveChild);
            } else {
               gEve->AddElement(eveChild);
            }

            break;

         case kTRACK:
            eveChild = drawTracks(collection, parent);
            eveChild->SetElementName(it->c_str());

            // Rename "All hits" to "<collection name> - All hits"
            // "All hits" contains all the hits of all the tracks in the
            // current collection. This is done inside drawTracks().

            try {

               // This can be replaced by (UNTESTED!)
               // listTreeItem = ListTree->FindItemByObj(ListTree->GetFirstItem(), eveChild)->GetNextSibling

               listTreeItem = gEve->GetListTree()->FindItemByPathname("/Event/All hits");

               if (listTreeItem) {
                  eveChild = static_cast<TEveElement*>( listTreeItem->GetUserData() );

                  it->append(" - All hits");
                  eveChild->SetElementName(it->c_str());
               }

            } catch (...) {
               //
            }

            break;

         default:
            break;
      }

   }
}

void VisTPC::drawHit(EVENT::TrackerHit* hit, TEvePointSet* parent, int wantTooltip) {

   if (!parent)
      throw "VisTPC::drawHit() : parent can't be null.";

   TEvePointSet      *point      = 0;

   const double      *pos   = hit->getPosition();
   const FloatVec    &covariance = hit->getCovMatrix();
  
   double position[3];

   switch (tpcParameters->getCoordinateType())
   {
     case PadRowLayout2D::CARTESIAN:
       position[0] = pos[0];
       position[1] = pos[1];
       position[2] = pos[2];
       break;
 
     case PadRowLayout2D::POLAR: 
       position[0] = pos[0]*cos(pos[1]);
       position[1] = pos[0]*sin(pos[1]);
       position[2] = pos[2];
       break;
      default:
         throw "Unknown coordinate type.";
   }

   std::stringstream tooltip;

   if (wantTooltip) {

      tooltip
         << "Position [mm]: ("
            << position[0] << ", "
            << position[1] << ", "
            << position[2] << ")"
            << std::endl
         << "Covariance: ("
            << covariance[0] << ", "
            << covariance[1] << ", "
            << covariance[2] << ")"
            << std::endl
         << "Energy [GeV]: " << hit->getEDep()    << std::endl
         << "Time [ns]: "    << hit->getTime()    << std::endl
         << "Type: "         << hit->getType()    << std::endl
         << "Quality: "      << hit->getQuality() << std::endl;

      point = new TEvePointSet();
      point->SetOwnIds(kTRUE);
      point->SetMarkerSize(0.2);
      point->SetMarkerStyle(2);
      point->SetMarkerColor(kMagenta);
      point->SetElementTitle(tooltip.str().c_str());

      parent->AddElement(point);

      parent = point;

   }

   //

   parent->SetNextPoint(position[0], position[1], position[2]);

   // Add one projected to the end cap

   parent->SetNextPoint(position[0], position[1], maxDriftLength + 1);
}

TEveElement* VisTPC::drawHits(IMPL::LCCollectionVec* hits, int wantTooltip) {

   TEvePointSet* ps  = new TEvePointSet();
   TrackerHit*   th;
   const double* pos = NULL;

   //===========================================================================

   ps->SetOwnIds(kTRUE);
   ps->SetMarkerColor(kMagenta);
   ps->SetMarkerSize(0.2);
   ps->SetMarkerStyle(2);

   int j = 0;
   for (
      LCCollectionVec::iterator inputHitsIter = hits->begin();
      inputHitsIter < hits->end();
      inputHitsIter++
   ) {
      if (maxPoints && j++ > maxPoints) {
         break;
      }

      th  = dynamic_cast<TrackerHit *>( *inputHitsIter );

      drawHit(th, ps, wantTooltip);
   }

   return ps;
}

TEveElement* VisTPC::drawHits(const EVENT::TrackerHitVec& hits, int wantTooltip) {

   TEvePointSet* ps  = new TEvePointSet();
   TrackerHit*   th;
   const double* pos = NULL;

   //===========================================================================

   ps->SetOwnIds(kTRUE);
   ps->SetMarkerColor(kMagenta);
   ps->SetMarkerSize(0.2);
   ps->SetMarkerStyle(2);

   int j = 0;
   for (
      TrackerHitVec::const_iterator inputHitsIter = hits.begin();
      inputHitsIter != hits.end();
      inputHitsIter++
   ) {
      if (maxPoints && j++ > maxPoints) {
         break;
      }

      th  = *inputHitsIter;

      drawHit(th, ps, wantTooltip);
   }

   return ps;
}

void VisTPC::drawModule(TPCModule * aModule, TEveElement* parent) {

   const std::vector<double>& extent = aModule->getLocalPadLayout().getPlaneExtent();
   double                     border = aModule->getBorderWidth();

   // Cartesian: [x_min, x_max, y_min, y_max]
   // Polar    : [r_min, r_max, phi_min, phi_max]

   float x1_min     = extent[0];
   float x1_max     = extent[1];
   float x2_min     = extent[2];
   float x2_max     = extent[3];

   float r_min, phi_min, r_max, phi_max = 0.;

   float pos_z      = 0;
   int   nDivisions = 100;

   TEveStraightLineSet* module_outline = new TEveStraightLineSet("Outline");
   std::stringstream    tooltip;

   //===========================================================================
   // Draw outline

   module_outline->SetLineColor(kBlue);
   module_outline->SetRnrSelf(options["module.drawOutlineByDefault"]);

   switch (aModule->getLocalPadLayout().getCoordinateType()) {
      case PadRowLayout2D::CARTESIAN:

         x1_min -= border;
         x1_max += border;
         x2_min -= border;
         x2_max += border;

         module_outline->AddLine(x1_min, x2_min, pos_z,
                                 x1_min, x2_max, pos_z);

         module_outline->AddLine(x1_max, x2_min, pos_z,
                                 x1_max, x2_max, pos_z);

         module_outline->AddLine(x1_min, x2_min, pos_z,
                                 x1_max, x2_min, pos_z);

         module_outline->AddLine(x1_min, x2_max, pos_z,
                                 x1_max, x2_max, pos_z);
         break;
      case PadRowLayout2D::POLAR:

         x1_min -= border;
         x1_max += border;

         // When not a full circle

         if (fabs((x2_max - x2_min) - 2*TMath::Pi()) > 1e-5) {

            if (!options["module.parallelBorders"]) {

               // Draw border at constant phi

               x2_min -= border / x1_min;
               x2_max += border / x1_min;

               drawCircularPolygon(module_outline, x1_min, x2_min, x2_max, nDivisions, pos_z);
               drawCircularPolygon(module_outline, x1_max, x2_min, x2_max, nDivisions, pos_z);

               module_outline->AddLine(x1_min*cos(x2_min), x1_min*sin(x2_min), pos_z,
                                       x1_max*cos(x2_min), x1_max*sin(x2_min), pos_z);

               module_outline->AddLine(x1_min*cos(x2_max), x1_min*sin(x2_max), pos_z,
                                       x1_max*cos(x2_max), x1_max*sin(x2_max), pos_z);

            } else {

               //=====================================
               // Draw border parallel to plane extent

               // Two arcs

               drawCircularPolygon(module_outline, x1_min, x2_min, x2_max, nDivisions, pos_z);
               drawCircularPolygon(module_outline, x1_max, x2_min, x2_max, nDivisions, pos_z);

               // Parallel border at right side with two lines perpendicular to
               // the border connecting to the arcs.

               phi_min = x2_min - border / x1_min;
               r_min   = x1_min / cos(border / x1_min);

               phi_max = x2_min - border / x1_max;
               r_max   = x1_max / cos(border / x1_max);

               module_outline->AddLine(r_min*cos(phi_min), r_min*sin(phi_min), pos_z,
                                       r_max*cos(phi_max), r_max*sin(phi_max), pos_z);

               module_outline->AddLine(x1_min*cos(x2_min), x1_min*sin(x2_min), pos_z,
                                       r_min*cos(phi_min), r_min*sin(phi_min), pos_z);

               module_outline->AddLine(x1_max*cos(x2_min), x1_min*sin(x2_min), pos_z,
                                       r_max*cos(phi_max), r_max*sin(phi_max), pos_z);


               // Parallel border at left side with two lines perpendicular to
               // the border connecting to the arcs.

               phi_min = x2_max + border / x1_min;
               r_min   = x1_min / cos(border / x1_min);

               phi_max = x2_max + border / x1_max;
               r_max   = x1_max / cos(border / x1_max);

               module_outline->AddLine(r_min*cos(phi_min), r_min*sin(phi_min), pos_z,
                                       r_max*cos(phi_max), r_max*sin(phi_max), pos_z);

               module_outline->AddLine(x1_min*cos(x2_max), x1_min*sin(x2_max), pos_z,
                                       r_min*cos(phi_min), r_min*sin(phi_min), pos_z);

               module_outline->AddLine(x1_max*cos(x2_max), x1_max*sin(x2_max), pos_z,
                                       r_max*cos(phi_max), r_max*sin(phi_max), pos_z);

            }

         } else {

            drawCircularPolygon(module_outline, x1_min, x2_min, x2_max, nDivisions, pos_z);
            drawCircularPolygon(module_outline, x1_max, x2_min, x2_max, nDivisions, pos_z);

         }

         break;

      default:
         throw "Unknown coordinate type.";
   }

   toGlobal(aModule, module_outline);

   if (parent) {

      parent->AddElement(module_outline);

      tooltip << parent->GetElementName() << " outline";
      module_outline->SetElementTitle(tooltip.str().c_str());

   } else {
      gEve->AddGlobalElement(module_outline);
   }

   //===========================================================================
   // Pads

   drawPads(aModule, parent);
}

void VisTPC::drawModules(TEveElement* parent) {

   // 1. Loop over all modules.
   // 2. Draw each module.
   // 3. Add tooltip for each module.

   //===========================================================================
   // Variables

   const std::vector<TPCModule *> & tpcModules = tpcParameters->getModules();
   TPCModule*                       tpcModule;
   TEveElementList*                 elements   = NULL;
   int                              i          = 0;

   //===========================================================================
   // Debug

   if (debug) {
      std::cout << tpcParameters->getNModules() << " Module(s)" << std::endl;
   }

   //===========================================================================
   // 1. Loop over all modules.

   for (
      std::vector<TPCModule *>::const_iterator iterator = tpcModules.begin();
      iterator != tpcModules.end();
      ++iterator
   ) {
      tpcModule = *iterator;

      elements = new TEveElementList(Form("Module %i", tpcModule->getModuleID()));
      elements->SetPickable(kFALSE);

      // 2. Draw each module.

      drawModule(tpcModule, elements);

      //

      if (parent) {
         parent->AddElement(elements);
      } else {
         gEve->AddGlobalElement(elements);
      }
   }

   //

   gEve->Redraw3D();

   // Finalize

   tpcModule = NULL;
   elements  = NULL;
}

void VisTPC::drawNearestPad(const gear::TPCModule& tpcModule, TEveElement* parent, int padIndex, const gear::Vector2D& position, int i) {
   // Initialize variables

   TEveElementList*     node      = new TEveElementList(Form("%i", i));
   TEvePointSet*        point     = new TEvePointSet(Form("Test hit %i", i));
   TEvePointSet*        padcenter = 0;
   TEveQuadSet*         pad       = new TEveQuadSet(Form("Pad %i", i));
   TEveStraightLineSet* line      = new TEveStraightLineSet(Form("Line %i", i));

   PadMetric padMetric;

   float     a_pos, b_pos;
   float     z_pos = maxDriftLength;

   Vector2D             padPos;
   const PadRowLayout2D &layout = tpcModule.getLocalPadLayout();

   std::stringstream tooltip;

   // Initialize new graphical elements

   point->SetOwnIds(kTRUE);
   point->SetMarkerColor(kRed);
   point->SetMarkerSize(0.2);
   point->SetMarkerStyle(2);

   pad->Reset(TEveQuadSet::kQT_FreeQuad, kTRUE, 32);
   pad->RefitPlex();
   pad->SetOwnIds(kTRUE);

   tooltip.str("");

   padcenter = new TEvePointSet(Form("padcenter %i", i));
   padcenter->SetOwnIds(kTRUE);
   padcenter->SetMarkerColor(kGreen);
   padcenter->SetMarkerSize(0.2);
   padcenter->SetMarkerStyle(2);

   line->SetLineColor(kYellow);
   //line->SetLineColor(kMagenta);

   // Get pad using random point in local coordinate

   a_pos = position[0];
   b_pos = position[1];

   tooltip << "Test hit: (" << a_pos << ", " << b_pos << ")" << std::endl;

   padPos   = layout.getPadCenter(padIndex);

   tooltip << "Padcenter: (" << padPos[0] << ", " << padPos[1] << ")" << std::endl;

   tooltip << "Difference: ("
           << a_pos - padPos[0] << ", "
           << b_pos - padPos[1] << ")"
           << std::endl;

   // Draw point in local coordinate

   if (layout.getCoordinateType() == PadRowLayout2D::POLAR) {
      float x, y;

      x = a_pos*cos(b_pos);
      y = a_pos*sin(b_pos);

      a_pos = x;
      b_pos = y;
   }

   tooltip << "Test hit (cartesian): (" << a_pos << ", " << b_pos << ")" << std::endl;

   point->SetNextPoint(a_pos, b_pos, 1);

   // Draw pad in local coordinate

   getLocalPadMetric(&tpcModule, padIndex, padMetric);

   padcenter->SetNextPoint(padMetric.center_x, padMetric.center_y, 1);

   tooltip << "padMetric.center(cartesian): (" << padMetric.center_x << ", " << padMetric.center_y << ")" << std::endl;

   drawPadQuad(padMetric, 0, pad);
   pad->QuadColor(kBlack);

   // Connect the pad and hit with a line

   line->AddLine(
      a_pos, b_pos, 1,
      padMetric.center_x, padMetric.center_y, 1
   );

   // Transform everything into global coordinate system

   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(point));
   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(pad));
   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(padcenter));
   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(line));

   //

   node->SetElementTitle(Form("%i\n%s", i, tooltip.str().c_str()));
   point->SetElementTitle(Form("Test hit %i\n%s", i, tooltip.str().c_str()));
   pad->SetElementTitle(Form("Pad %i\n%s", i, tooltip.str().c_str()));
   padcenter->SetElementTitle(Form("Pad center %i\n%s", i, tooltip.str().c_str()));
   line->SetElementTitle(Form("Line %i\n%s", i, tooltip.str().c_str()));

   node->AddElement(point);
   node->AddElement(pad);
   node->AddElement(padcenter);
   node->AddElement(line);

   if (parent)
      parent->AddElement(node);
   else
      gEve->AddElement(node);
}

void VisTPC::drawPadBox(const PadMetric& pad, float z_front, float z_back, TEveBoxSet* bs) {
   float boxCorners[24];

   // top back left
   boxCorners[0]  = pad.corners_cartesian[6];
   boxCorners[1]  = pad.corners_cartesian[7];
   boxCorners[2]  = z_back;

   // top back right
   boxCorners[3]  = pad.corners_cartesian[4];
   boxCorners[4]  = pad.corners_cartesian[5];
   boxCorners[5]  = z_back;

   // top front right
   boxCorners[6]  = pad.corners_cartesian[4];
   boxCorners[7]  = pad.corners_cartesian[5];
   boxCorners[8]  = z_front;

   // top front left
   boxCorners[9]  = pad.corners_cartesian[6];
   boxCorners[10] = pad.corners_cartesian[7];
   boxCorners[11] = z_front;

   // bottom back left
   boxCorners[12] = pad.corners_cartesian[0];
   boxCorners[13] = pad.corners_cartesian[1];
   boxCorners[14] = z_back;

   // bottom back right
   boxCorners[15] = pad.corners_cartesian[2];
   boxCorners[16] = pad.corners_cartesian[3];
   boxCorners[17] = z_back;

   // bottom front right
   boxCorners[18] = pad.corners_cartesian[2];
   boxCorners[19] = pad.corners_cartesian[3];
   boxCorners[20] = z_front;

   // bottom front left
   boxCorners[21] = pad.corners_cartesian[0];
   boxCorners[22] = pad.corners_cartesian[1];
   boxCorners[23] = z_front;

   bs->AddBox(boxCorners);
}

void VisTPC::drawPadQuad(const PadMetric& pad, float z, TEveQuadSet* qs) {
   float quadCorners[12];

   // bottom left
   quadCorners[0]  = pad.corners_cartesian[0];
   quadCorners[1]  = pad.corners_cartesian[1];
   quadCorners[2]  = z;

   // bottom right
   quadCorners[3]  = pad.corners_cartesian[2];
   quadCorners[4]  = pad.corners_cartesian[3];
   quadCorners[5]  = z;

   // top right
   quadCorners[6]  = pad.corners_cartesian[4];
   quadCorners[7]  = pad.corners_cartesian[5];
   quadCorners[8]  = z;

   // top left
   quadCorners[9]  = pad.corners_cartesian[6];
   quadCorners[10] = pad.corners_cartesian[7];
   quadCorners[11] = z;

   qs->AddQuad(quadCorners);
}

void VisTPC::drawPads(TPCModule * aModule, TEveElement* parent) {

   // 1. Loop over all rows.
   // 2. Draw radial lines and a center marker of each pad.
   // 3. Draw circular lines of each row.

   //===========================================================================
   // Variables

   const PadRowLayout2D &layout = aModule->getLocalPadLayout();
   PadMetric            pad;
   std::stringstream    tooltip;

   // 1. Loop over all rows.

   int nRows = layout.getNRows();
       //nRows = nRows > 25 ? 25 : nRows; // CHANGE: Just for rapid prototyping. Showing full module is slow process

   // 2. Draw radial lines of each pad.

   const std::vector<int> *         padsInRow = NULL;
   std::vector<int>::const_iterator iterator;

   TEveStraightLineSet* ls = new TEveStraightLineSet("Pad borders");
   TEvePointSet*        ps = new TEvePointSet("Pad centers");

   float x_pos, y_pos = 0.;
   float z_pos        = 0.;

   // 3. Draw circular lines of each row.

   const std::vector<double>& extent = layout.getPlaneExtent();

   int    nPadsInRow       = 0;
   int    nDivisions       = 0;
   double row_extent_min   = 0.;
   double row_extent_max   = 0.;
   double padRowHeightDiff = 0.;
   double padGap           = 0.;

   // Debug output

   //std::cout << "Rows: " << nRows << std::endl;

   //===========================================================================

   ls->SetRnrSelf(options["module.drawPadsByDefault"]);
   ls->SetLineColor(kRed);
   ls->SetPickable(kFALSE);

   ps->SetRnrSelf(kFALSE);
   ps->SetOwnIds(kTRUE);
   ps->SetMarkerColor(kGreen);
   ps->SetMarkerSize(1);
   ps->SetMarkerStyle(1);
   ps->SetPickable(kFALSE);

   //

   switch (layout.getCoordinateType()) {
      case PadRowLayout2D::CARTESIAN:
         row_extent_min = extent[0];
         row_extent_max = extent[1];
         break;

      case PadRowLayout2D::POLAR:
         row_extent_min = extent[2];
         row_extent_max = extent[3];
         break;

      default:
         throw "Unknown coordinate type.";
   }

   // 1. Loop over all rows.

   for (int i = 0; i < nRows; i++) {
      padsInRow  = &layout.getPadsInRow(i);
      iterator   = padsInRow->begin();
      nPadsInRow = padsInRow->size();

      x_pos = 0;
      y_pos = 0;
      z_pos = 0;

      // Debug

      //std::cout <<"Cols[" << i << "]: " << nPadsInRow << std::endl;

      // 2. Draw radial lines of each pad.

      padGap = layout.getPadPitch(*(padsInRow->begin())) - layout.getPadWidth(*(padsInRow->begin()));

      while (iterator != padsInRow->end()) {
         getLocalPadMetric(aModule, *iterator, pad);

         //if (debug) {
         //   std::cout << "Pad corners: (" << pad.corners_cartesian[0] << ", " << pad.corners_cartesian[1] << ")" << std::endl
         //             << "             (" << pad.corners_cartesian[2] << ", " << pad.corners_cartesian[3] << ")" << std::endl
         //             << "             (" << pad.corners_cartesian[4] << ", " << pad.corners_cartesian[5] << ")" << std::endl
         //             << "             (" << pad.corners_cartesian[6] << ", " << pad.corners_cartesian[7] << ")" << std::endl;
         //}

         x_pos = pad.center_x;
         y_pos = pad.center_y;

         switch (layout.getCoordinateType()) {
            case PadRowLayout2D::CARTESIAN:
               if (iterator == padsInRow->begin() || padGap > 1e-6) {
                  ls->AddLine(pad.corners_cartesian[0], pad.corners_cartesian[1], z_pos,
                              pad.corners_cartesian[6], pad.corners_cartesian[7], z_pos);
               }

               ls->AddLine(pad.corners_cartesian[2], pad.corners_cartesian[3], z_pos,
                           pad.corners_cartesian[4], pad.corners_cartesian[5], z_pos);

               break;

            case PadRowLayout2D::POLAR:
               if (iterator == padsInRow->begin() || padGap > 1e-6) {
                  ls->AddLine(pad.corners_cartesian[0], pad.corners_cartesian[1], z_pos,
                              pad.corners_cartesian[2], pad.corners_cartesian[3], z_pos);
               }

               ls->AddLine(pad.corners_cartesian[6], pad.corners_cartesian[7], z_pos,
                           pad.corners_cartesian[4], pad.corners_cartesian[5], z_pos);

               break;

            default:
               throw "Unknown coordinate type.";
         }

         // 2. and the center marker

         ps->SetNextPoint(x_pos, y_pos, z_pos);

         //

         ++iterator;
      }

      // 3. Draw row lines of each row.

      padRowHeightDiff = layout.getRowHeight(i) - layout.getPadHeight(*(padsInRow->begin()));

      switch (layout.getCoordinateType()) {
         case PadRowLayout2D::CARTESIAN:
            ls->AddLine(row_extent_min, pad.corners_cartesian[1], z_pos,
                        row_extent_max, pad.corners_cartesian[1], z_pos);

            if (i == nRows - 1 || padRowHeightDiff > 1e-6) {
               ls->AddLine(row_extent_min, pad.corners_cartesian[7], z_pos,
                           row_extent_max, pad.corners_cartesian[7], z_pos);
            }

            break;

         case PadRowLayout2D::POLAR:
            nDivisions = ceil(nPadsInRow / 5.);

            drawCircularPolygon(ls, pad.corners_polar[0], row_extent_min, row_extent_max, nDivisions, z_pos);

            if (i == nRows - 1 || padRowHeightDiff > 1e-6) {
               drawCircularPolygon(ls, pad.corners_polar[4], row_extent_min, row_extent_max, nDivisions, z_pos);
            }

            break;

         default:
            throw "Unknown coordinate type.";
      }

   }

   toGlobal(aModule, ls);
   toGlobal(aModule, ps);

   // Add Line set to Eve

   if (parent) {

      parent->AddElement(ls);
      parent->AddElement(ps);

      tooltip << parent->GetElementName() << " pads";
      ls->SetElementTitle(tooltip.str().c_str());

      tooltip.str("");
      tooltip << parent->GetElementName() << " pads";
      ps->SetElementTitle(tooltip.str().c_str());

   } else {
      gEve->AddGlobalElement(ls);
      gEve->AddGlobalElement(ps);
   }
}

TEveElement* VisTPC::drawPulses(IMPL::LCCollectionVec* pulses) {

   // 1. Loop over all pulses.
   // 2. Calculate the 8 box corners of the pad projected into 3D space.
   // 3. Draw the box and marker.
   // 4. Draw global max charges.

   //===========================================================================
   // Variables

   TEveRGBAPalette* pal = new TEveRGBAPalette();

   TEveElementList* topElement = new TEveElementList();

   TEveElementList* el  = new TEveElementList("Pulses");
   TEveBoxSet*      bs  = new TEveBoxSet     ("Boxes");
   TEvePointSet*    ps  = new TEvePointSet   ("Markers");

   TEveElementList* el2 = new TEveElementList("Max charge deposit");
   TEveBoxSet*      bs2 = new TEveBoxSet     ("3D lego");
   TEveQuadSet*     qs2 = new TEveQuadSet    ("2D projection");
   TEvePointSet*    ps2 = new TEvePointSet   ("Markers");

   //

   TrackerPulse*            tp;
   TrackerData*             td;
   int                      moduleId        = 0;
   int                      padId           = 0;
   float                    time            = 0.;
   float                    localMaxCharge  = 0.;
   float                    overallCharge   = 0.;
   float                    startTime       = 0.;
   float                    endTime         = 0.;

   typedef std::map<int, float>        MapPadCharge;
   typedef std::map<int, MapPadCharge> MapModulePadCharge;

   MapModulePadCharge globalCharge;
   MapPadCharge*      moduleCharge;

   //

   const TPCModule* tpcModule = NULL;

   PadMetric pad;
   float     x_pos, y_pos = 0;
   float     z_pos        = maxDriftLength;

   //

   //===========================================================================

   // Setup the visualization stuffs

   gStyle->SetPalette(1, 0);
   pal->SetMinMax(0, 1023);
   pal->SetInterpolate(kTRUE);
   pal->SetOverColor(kWhite);
   pal->SetOverflowAction(TEveRGBAPalette::kLA_Mark);
   pal->SetUnderColor(kBlack);
   pal->SetUnderflowAction(TEveRGBAPalette::kLA_Mark);

   bs->Reset(TEveBoxSet::kBT_FreeBox, kFALSE, 64);
   //bs->SetRnrSelf(kFALSE);
   bs->SetMainAlpha(0.9);
   bs->SetPalette(pal);

   ps->SetRnrSelf(kFALSE);
   //ps->SetOwnIds(kTRUE);
   ps->SetMarkerColor(kYellow);
   ps->SetMarkerSize(0.2);
   ps->SetMarkerStyle(2);

   bs2->Reset(TEveBoxSet::kBT_FreeBox, kFALSE, 64);
   bs2->SetRnrSelf(kFALSE);
   //bs2->SetMainAlpha(0.5);
   bs2->SetPalette(pal);

   ps2->SetRnrSelf(kFALSE);
   //ps2->SetOwnIds(kTRUE);
   ps2->SetMarkerColor(kYellow);
   ps2->SetMarkerSize(0.2);
   ps2->SetMarkerStyle(2);

   qs2->Reset(TEveQuadSet::kQT_FreeQuad, kFALSE, 32);
   qs2->SetPalette(pal);

   // 1. Loop over all pulses.

   int j = 0;
   for (
      LCCollectionVec::iterator inputPulsIter = pulses->begin();
      inputPulsIter < pulses->end();
      inputPulsIter++
   ) {
      if (maxPoints && j++ > maxPoints) {
         break;
      }

      tp = dynamic_cast<TrackerPulse *>( *inputPulsIter );
      td = tp->getTrackerData();

      //

      moduleId       = tp->getCellID1();
      padId          = tp->getCellID0();

      time           = tp->getTime();
      localMaxCharge = 0;
      startTime      = td->getTime() / readoutFrequency;
      endTime        = startTime + td->getChargeValues().size() / readoutFrequency;
      tpcModule      = &tpcParameters->getModule(moduleId);

      // 2. Calculate the 8 box corners of the pad projected into 3D space.

      getGlobalPadMetric(tpcModule, padId, pad);

      x_pos = pad.center_x;
      y_pos = pad.center_y;

      // Find max charge

      for (
         FloatVec::const_iterator i = td->getChargeValues().begin();
         i != td->getChargeValues().end();
         ++i
      ) {
         if (*i > localMaxCharge) {
            localMaxCharge = *i;
         }
      }

      if (!globalCharge.count(moduleId)) {
         globalCharge[moduleId] = MapPadCharge();
      }

      moduleCharge = &globalCharge[moduleId];

      if (!moduleCharge->count(padId)) {
         (*moduleCharge)[padId] = 0;
      }

      if ((*moduleCharge)[padId] < localMaxCharge) {
         (*moduleCharge)[padId] = localMaxCharge;
      }

      // 3. Draw the marker and box.

      ps->SetNextPoint(x_pos, y_pos, z_pos - startTime* driftVelocity);

      drawPadBox(
         pad,
         z_pos - endTime * driftVelocity,
         z_pos - startTime * driftVelocity,
         bs
      );
      bs->DigitValue(localMaxCharge);
   }

   // 4. Draw global max charges.

   for (
      MapModulePadCharge::iterator j = globalCharge.begin();
      j != globalCharge.end();
      ++j
   ) {

      moduleId  = (*j).first;
      tpcModule = &tpcParameters->getModule(moduleId);

      for (
         MapPadCharge::iterator k = (*j).second.begin();
         k != (*j).second.end();
         ++k
      ) {

         padId          = (*k).first;
         localMaxCharge = (*k).second;

         getGlobalPadMetric(tpcModule, padId, pad);

         x_pos = pad.center_x;
         y_pos = pad.center_y;

         // 2D

         drawPadQuad(pad, z_pos, qs2);
         qs2->DigitValue(localMaxCharge);

         ps2->SetNextPoint(x_pos, y_pos, z_pos + 1);

         // 3D lego

         drawPadBox(pad, z_pos, z_pos + localMaxCharge / 10, bs2);
         bs2->DigitValue(localMaxCharge);

         //ps2->SetNextPoint(x_pos, y_pos, z_pos + 50 + localMaxCharge);
      }

   }

   //

   bs->RefitPlex();
   bs2->RefitPlex();
   qs2->RefitPlex();

   topElement->AddElement(el);
   el->AddElement(bs);
   el->AddElement(ps);

   topElement->AddElement(el2);
   el2->AddElement(bs2);
   el2->AddElement(ps2);
   el2->AddElement(qs2);

   return topElement;
}

void VisTPC::drawTPC() {

   volumes["top"] = gGeoManager->MakeBox("Detector", media["Vacuum"],
                1.5*outerFieldCageRadius,
                1.5*outerFieldCageRadius,
                3*maxDriftLength);

   // Make a tube as TPC.
   volumes["TPC"] = gGeoManager->MakeTube("TPC", media["Al"], innerFieldCageRadius,
                 outerFieldCageRadius, maxDriftLength);

   volumes["TPC"]->SetTransparency(75);
   volumes["top"]->AddNode(volumes["TPC"], 0, new TGeoTranslation(0,0,0));

   gGeoManager->SetTopVolume(volumes["top"]);

   TGeoNode*       node = gGeoManager->GetTopNode();
   TEveGeoTopNode* en   = new TEveGeoTopNode(gGeoManager, node);

   en->SetVisLevel(1);
   en->GetNode()->GetVolume()->SetVisibility(kFALSE);

   gEve->AddGlobalElement(en);

   drawModules(en);

   //---------------------------------------------------------------------------
   // Close the geometry

   //gGeoManager->CloseGeometry();

   gEve->Redraw3D();
}

TEveElement* VisTPC::drawTracks(LCCollectionVec* tracks, TEveElement* parent) {

   //---------------------------------------------------------------------------
   // Contents
   //---------------------------------------------------------------------------
   // 1. Declarations & initializations
   // 2. Init stuffs
   // 3. Draw tracks
   // 4. Add tracks to parent or global
   //---------------------------------------------------------------------------

   //---------------------------------------------------------------------------
   // 1. Declarations & initializations
   //---------------------------------------------------------------------------

   const float         cm             = 1.;
   const float         mm             = 0.1 * cm;

   Track               *lcio_track    = NULL;
   const TrackerHitVec *lcio_hits     = NULL;

   //TEveElementList     *eveNodeEvent  = new TEveElementList(name);

   TEveTrackList       *track_list    = new TEveTrackList();
   TEveTrackPropagator *propagator    = track_list->GetPropagator();
   TEveRecTrack        *eve_track_rec = NULL;
   TEveTrack           *eve_track     = NULL;
   TEveStraightLineSet *eve_line      = NULL;
   TEveElement         *eve_hits      = NULL;
   TEveElementList     *eve_all_hits  = new TEveElementList("All hits");

   Vector3D            b_field(0,0,0);

   float               p_transversal  = 0.;

   const float         *reference_point               = NULL;
   float               start_point[3];
   float               d0, phi0, z0, tanLambda, omega = 0.;
   int                 charge                         = 0;

   int j = 0;

   //   const std::vector<double>& extent = tpcParameters->getPlaneExtent();
   //   float                      rmin   = 0.;
   // float                      rmax   = 0.;

   // 30 = TrackFitterBase::FITFAILEDBIT
   const int FITFAILEDBIT = 30;
   bool      fitFailed    = 0;

   //---------------------------------------------------------------------------
   // 2. Init stuffs
   //---------------------------------------------------------------------------

   try{
      b_field = gearMgr->getBField().at(Vector3D(0, 0, 0));
   }
   catch (UnknownParameterException & e) {
     std::cout << "WARNING:" << e.what() << std::endl;
     std::cout << "B field is assumed to be 0" << std::endl;
   }

   // Set properties of the propagator used for reconstructing the tracks.

   propagator->SetMagFieldObj(
      new TEveMagFieldConst(
         b_field.x(),
         b_field.y(),
         b_field.z()
      )
   );

   propagator->SetMaxR(outerFieldCageRadius);
   propagator->SetMaxZ(maxDriftLength);
   propagator->SetMaxOrbs(50);

   //---------------------------------------------------------------------------
   // 3. Draw tracks
   //---------------------------------------------------------------------------

   for (
      LCCollectionVec::iterator tracksIter = tracks->begin();
      tracksIter != tracks->end();
      ++tracksIter
   ) {
      if (maxPoints && j++ > maxPoints) {
         break;
      }

      lcio_track = dynamic_cast<Track*>( *tracksIter );

      //

      reference_point = lcio_track->getReferencePoint();
      d0              = lcio_track->getD0();
      omega           = lcio_track->getOmega();
      phi0            = lcio_track->getPhi();
      tanLambda       = lcio_track->getTanLambda();
      z0              = lcio_track->getZ0();

      start_point[0]  = reference_point[0] - d0 * sin(phi0);
      start_point[1]  = reference_point[1] + d0 * cos(phi0);
      start_point[2]  = reference_point[2] + z0;

      fitFailed       = (lcio_track->getType() >> FITFAILEDBIT) & 0x1 ? 1 : 0;
      //fitFailed       = lcio_track->getType() & (1<<FITFAILEDBIT) ;

      // Do some necessary calculations

      if (omega) {
         charge = b_field.z() / omega >= 0. ? 1 : -1;
      } else {
         charge = 0;
      }

      // MouseOver tooltip

      std::stringstream tooltip;

      tooltip << std::left << "d0 = "        << d0        << std::endl
              << std::left << "omega = "     << omega     << std::endl
              << std::left << "phi0 = "      << phi0      << std::endl
              << std::left << "tanLambda = " << tanLambda << std::endl
              << std::left << "starting point = ("
                           << start_point[0] << ", "
                           << start_point[1] << ", "
                           << start_point[2] << ")"       << std::endl
              << std::left << "fitFailed = " << fitFailed << std::endl;

      // First draw the hits, otherwise it would require a second if-else
      // construction

      lcio_hits = &lcio_track->getTrackerHits();

      eve_hits = drawHits(*lcio_hits);
      ((TEvePointSet*) eve_hits)->SetMarkerColor(!fitFailed ? kRed : kGreen);
      ((TEvePointSet*) eve_hits)->SetRnrSelfChildren(kFALSE, kFALSE);

      eve_all_hits->AddElement(eve_hits);

      // Make a distinction between a magnetic field (tracks)
      // and an absent magnetic field (straight line) because
      // TEveTrack apparently doesn't work with absent magnetic
      // field.
      //
      // Without a B field there is no p_transversal -> No track.
      // With omega=0 p_transversal is infinite (divide-by-zero).
      //
      // TODO: Investigate the possibility of using TEveTrackPropagator for
      // both curved and straight tracks. It has code to do that, just need to
      // find out how to implement it here. It would simplify the code here.

      if (fabs(b_field.z()) > 0. && omega) {

         // in case the particle does not come from the IP just starting
         // the particle at the PCA is not enough. Draw a second arm into
         // the direction where it came from (negative momentum with same start point)
         // This is only done if there is no inner field cage, otherwise the
         // particle will (most probably) come from the vertex.

         for (
            int direction = 0; 
            direction < ( innerFieldCageRadius==0?2:1 );
            direction++
         ) {

            //-------------------------------------
            // Initial vertex location of the track

            eve_track_rec = new TEveRecTrack();
            eve_track_rec->fV.Set(
               start_point[0],
               start_point[1],
               start_point[2]
            );

            //-----------------
            // Initial momentum

            // 1. fV without *mm, but fP with *mm.. not really consistent.
            // 2. When two tracks are needed the drawn, the first track gets the
            // correct momentum and the second track opposite momentum and
            // charge as if it's an anti particle. This will create two tracks
            // that goes smooth in each other such that it looks like one
            // track.

            p_transversal = 3. * pow(10, -2) * fabs(b_field.z() / omega);
            
            eve_track_rec->fP.Set(
               static_cast<float>( (direction?-1:1) * p_transversal * cos(phi0) ) * mm,
               static_cast<float>( (direction?-1:1) * p_transversal * sin(phi0) ) * mm,
               static_cast<float>( (direction?-1:1) * p_transversal * tanLambda ) * mm
            );
            eve_track_rec->fSign = (direction?1:-1) * charge;

            //

            eve_track = new TEveTrack(eve_track_rec, propagator);
            eve_track->SetElementTitle(tooltip.str().c_str());
            eve_track->SetLineColor(fitFailed ? kRed : kGreen);
            eve_track->SetLineWidth(1);

            //

            track_list->AddElement(eve_track);
            eve_track->MakeTrack();

            eve_track->AddElement(eve_hits);

         }

      } else {

    if ( eve_line == NULL ) {
      eve_line = new TEveStraightLineSet();
    }

         // Calculate how much we have to multiply our vector to hit the
         // detector's boundaries.

         float scaling_min_xy = sqrt(pow(start_point[0] + cos(phi0) * innerFieldCageRadius, 2) 
                 + pow(start_point[1] + sin(phi0) * innerFieldCageRadius, 2));
         float scaling_min_z  = (maxDriftLength - start_point[2]) / tanLambda;
         float scaling_min    = tanLambda * scaling_min_xy > maxDriftLength ? scaling_min_z : scaling_min_xy;

         float scaling_max_xy = sqrt(pow(start_point[0] + cos(phi0) * outerFieldCageRadius, 2) 
                 + pow(start_point[1] + sin(phi0) * outerFieldCageRadius, 2));
         float scaling_max_z  = (maxDriftLength - start_point[2]) / tanLambda;
         float scaling_max    = tanLambda * scaling_max_xy > maxDriftLength ? scaling_max_z : scaling_max_xy;

         eve_line->AddLine(
            start_point[0] + cos(phi0) * scaling_min,
            start_point[1] + sin(phi0) * scaling_min,
            start_point[2] + tanLambda * scaling_min,
            start_point[0] + cos(phi0) * scaling_max,
            start_point[1] + sin(phi0) * scaling_max,
            start_point[2] + tanLambda * scaling_max
         );

         eve_line->SetElementTitle(tooltip.str().c_str());
         eve_line->SetLineColor(fitFailed ? kRed : kGreen);
         eve_line->SetLineWidth(1);

         //track_list->AddElement(eve_line);

         eve_line->AddElement(eve_hits);
      
      } // end if

   } // end for

   //---------------------------------------------------------------------------
   // 4. Add tracks to parent or global
   //---------------------------------------------------------------------------

   if (parent) {
      if (track_list)   parent->AddElement(track_list);
      if (eve_line)     parent->AddElement(eve_line);
      if (eve_all_hits) parent->AddElement(eve_all_hits);
   } else {
      if (track_list)   gEve->AddElement(track_list);
      if (eve_line)     gEve->AddElement(eve_line);
      if (eve_all_hits) gEve->AddElement(eve_all_hits);
   }

   return track_list;
}

void VisTPC::drawTwoConnectingPads(const gear::TPCModule& tpcModule, TEveElement* parent, int pad1, int pad2, LeftRight_enum leftright) {

   const PadRowLayout2D &layout    = tpcModule.getLocalPadLayout();

   TEveElementList*     eNode      = 0;
   TEveQuadSet*         ePad       = new TEveQuadSet(Form("Pad %i", pad1));
   TEveQuadSet*         eNeighbour = new TEveQuadSet(Form("Neighbour %i - %i", pad1, pad2));
   TEveStraightLineSet* eLine      = new TEveStraightLineSet(Form("Line %i", pad1));

   if (leftright == kLEFT) {
      eNode = new TEveElementList(Form("Left %i", pad1));
   } else {
      eNode = new TEveElementList(Form("Right %i", pad1));
   }

   PadMetric padMetric;
   Vector2D  padPos       = layout.getPadCenter(pad1);
   Vector2D  neighbourPos = layout.getPadCenter(pad2);

   //

   ePad->Reset(TEveQuadSet::kQT_FreeQuad, kTRUE, 32);
   ePad->RefitPlex();
   ePad->SetOwnIds(kTRUE);

   eNeighbour->Reset(TEveQuadSet::kQT_FreeQuad, kTRUE, 32);
   eNeighbour->RefitPlex();
   eNeighbour->SetOwnIds(kTRUE);

   // Draw test pad

   getLocalPadMetric(&tpcModule, pad1, padMetric);
   drawPadQuad(padMetric, 0, ePad);
   ePad->QuadColor(kYellow);

   // Draw neighbour pad

   getLocalPadMetric(&tpcModule, pad2, padMetric);
   drawPadQuad(padMetric, 0, eNeighbour);
   eNeighbour->QuadColor(kRed);

   // Draw connecting line

   if (layout.getCoordinateType() == PadRowLayout2D::POLAR) {
      float x, y;

      x = padPos[0]*cos(padPos[1]);
      y = padPos[0]*sin(padPos[1]);

      padPos[0] = x;
      padPos[1] = y;

      x = neighbourPos[0]*cos(neighbourPos[1]);
      y = neighbourPos[0]*sin(neighbourPos[1]);

      neighbourPos[0] = x;
      neighbourPos[1] = y;
   }

   eLine->AddLine(
      padPos[0], padPos[1], 1,
      neighbourPos[0], neighbourPos[1], 1
   );
   eLine->SetLineColor(kCyan);

   //

   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(ePad));
   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(eNeighbour));
   toGlobal(&tpcModule, dynamic_cast<TEveElement*>(eLine));

   //

   eNode->AddElement(ePad);
   eNode->AddElement(eNeighbour);
   eNode->AddElement(eLine);

   if (parent)
      parent->AddElement(eNode);
   else
      gEve->AddElement(eNode);
}

const VisTPC::LCCollectionType_enum VisTPC::getCollectionType(const std::string& collectionType) {
        if (!strcmp(collectionType.c_str(), EVENT::LCIO::TRACKERHIT))        {return kHIT;}
   else if (!strcmp(collectionType.c_str(), EVENT::LCIO::TRACKERPULSE))      {return kPULSE;}
   else if (!strcmp(collectionType.c_str(), EVENT::LCIO::TRACK))             {return kTRACK;}

   return kUNDEFINEDTYPE;
}

void VisTPC::getLocalPadMetric(const TPCModule * aModule, const int padIndex, PadMetric& pad) {
   const PadRowLayout2D* layout = &aModule->getLocalPadLayout();
   Vector2D              center = layout->getPadCenter(padIndex);

   // TODO
   // Add caching

   //

   pad.height = layout->getPadHeight(padIndex);
   pad.width  = layout->getPadWidth(padIndex);
  
   float x1_min = 0.;
   float x1_max = 0.;
   float x2_min = 0.;
   float x2_max = 0.;

   switch (layout->getCoordinateType()) {
      case PadRowLayout2D::CARTESIAN:

         pad.coordinate_system = PadRowLayout2D::CARTESIAN;

         x1_min = center[0] - pad.width / 2.;
         x1_max = center[0] + pad.width / 2.;
         x2_min = center[1] - pad.height / 2.;
         x2_max = center[1] + pad.height / 2.;

         pad.center_x = center[0];
         pad.center_y = center[1];

         pad.corners_cartesian[0] = x1_min;
         pad.corners_cartesian[1] = x2_min;
         pad.corners_cartesian[2] = x1_max;
         pad.corners_cartesian[3] = x2_min;
         pad.corners_cartesian[4] = x1_max;
         pad.corners_cartesian[5] = x2_max;
         pad.corners_cartesian[6] = x1_min;
         pad.corners_cartesian[7] = x2_max;

         // Polar coordinates are not needed in drawings when
         // we work in cartesian coordinate system.

         pad.corners_polar[0] = 0;
         pad.corners_polar[1] = 0;
         pad.corners_polar[2] = 0;
         pad.corners_polar[3] = 0;
         pad.corners_polar[4] = 0;
         pad.corners_polar[5] = 0;
         pad.corners_polar[6] = 0;
         pad.corners_polar[7] = 0;

         break;

      case PadRowLayout2D::POLAR:

         x1_min = center[0] - pad.height / 2.;
         x1_max = center[0] + pad.height / 2.;
         x2_min = center[1] - pad.width / 2.;
         x2_max = center[1] + pad.width / 2.;

         pad.coordinate_system = PadRowLayout2D::POLAR;

         pad.center_x = center[0] * cos(center[1]);
         pad.center_y = center[0] * sin(center[1]);

         pad.center_r = center[0];
         pad.center_phi = center[1];

         pad.corners_polar[0] = x1_min;
         pad.corners_polar[1] = x2_min;
         pad.corners_polar[2] = x1_max;
         pad.corners_polar[3] = x2_min;
         pad.corners_polar[4] = x1_max;
         pad.corners_polar[5] = x2_max;
         pad.corners_polar[6] = x1_min;
         pad.corners_polar[7] = x2_max;

         pad.corners_cartesian[0] = x1_min*cos(x2_min);
         pad.corners_cartesian[1] = x1_min*sin(x2_min);
         pad.corners_cartesian[2] = x1_max*cos(x2_min);
         pad.corners_cartesian[3] = x1_max*sin(x2_min);
         pad.corners_cartesian[4] = x1_max*cos(x2_max);
         pad.corners_cartesian[5] = x1_max*sin(x2_max);
         pad.corners_cartesian[6] = x1_min*cos(x2_max);
         pad.corners_cartesian[7] = x1_min*sin(x2_max);

         break;
      default:
         throw "Unknown coordinate type.";
   }
}

void VisTPC::getGlobalPadMetric(const TPCModule * aModule, const int padIndex, PadMetric& pad) {
   double         angle  = aModule->getAngle();
   const Vector2D offset = aModule->getOffset();
   Vector2D       center = aModule->getPadCenter(padIndex);


   //Vector2D       global;

   // TODO
   // Add caching

   //

   getLocalPadMetric(aModule, padIndex, pad);

   //

   // To global coordinate 

   switch (pad.coordinate_system) {
      case PadRowLayout2D::CARTESIAN:
         pad.center_x = center[0];
         pad.center_y = center[1];
         break;

      case PadRowLayout2D::POLAR:
         pad.center_x = center[0] * cos(center[1]);
         pad.center_y = center[0] * sin(center[1]);
         break;
      default:
         throw "Unknown coordinate type.";
   }

   // Rotation first..

   PadMetric tempPad = pad; 
   double sinangle   = sin(angle);
   double cosangle   = cos(angle);
   tempPad.corners_cartesian[0] = pad.corners_cartesian[0] * cosangle - pad.corners_cartesian[1] * sinangle;
   tempPad.corners_cartesian[1] = pad.corners_cartesian[0] * sinangle + pad.corners_cartesian[1] * cosangle;
   tempPad.corners_cartesian[2] = pad.corners_cartesian[2] * cosangle - pad.corners_cartesian[3] * sinangle;
   tempPad.corners_cartesian[3] = pad.corners_cartesian[2] * sinangle + pad.corners_cartesian[3] * cosangle;
   tempPad.corners_cartesian[4] = pad.corners_cartesian[4] * cosangle - pad.corners_cartesian[5] * sinangle;
   tempPad.corners_cartesian[5] = pad.corners_cartesian[4] * sinangle + pad.corners_cartesian[5] * cosangle;
   tempPad.corners_cartesian[6] = pad.corners_cartesian[6] * cosangle - pad.corners_cartesian[7] * sinangle;
   tempPad.corners_cartesian[7] = pad.corners_cartesian[6] * sinangle + pad.corners_cartesian[7] * cosangle;
   
   pad = tempPad;

   if (pad.coordinate_system == PadRowLayout2D::POLAR) {
      pad.corners_polar [1] += angle;
      pad.corners_polar [3] += angle;
      pad.corners_polar [5] += angle;
      pad.corners_polar [7] += angle;
   }

   // and then translation.

   pad.corners_cartesian[0] += offset[0]; pad.corners_cartesian[1] += offset[1];
   pad.corners_cartesian[2] += offset[0]; pad.corners_cartesian[3] += offset[1];
   pad.corners_cartesian[4] += offset[0]; pad.corners_cartesian[5] += offset[1];
   pad.corners_cartesian[6] += offset[0]; pad.corners_cartesian[7] += offset[1];

   if (pad.coordinate_system == PadRowLayout2D::POLAR) {
      // TODO for pad.corners_polar[]
   }
}

void VisTPC::testLeftRightNeighbour(int moduleID) {
   Vector2D             padPos;
   const TPCModule      &tpcModule   = tpcParameters->getModule(moduleID);
   const PadRowLayout2D &layout      = tpcModule.getLocalPadLayout();

   std::vector<int>     pads;

   //

   for (
      int nRow = 0;
      nRow < layout.getNRows();
      ++nRow
   ) {

      const std::vector<int> &padsInRow = layout.getPadsInRow(nRow);

      for (
         std::vector<int>::const_iterator it = padsInRow.begin();
         it != padsInRow.end();
         ++it
      ) {

         pads.push_back(*it);

      }

   }

   testLeftRightNeighbourWithPads(moduleID, pads, kLEFT);
   testLeftRightNeighbourWithPads(moduleID, pads, kRIGHT);
}

void VisTPC::testLeftRightNeighbourWithPads(int moduleID, const std::vector<int>& pads, LeftRight_enum leftright) {
   // Initialize variables

   const TPCModule      &tpcModule   = tpcParameters->getModule(moduleID);
   const PadRowLayout2D &layout      = tpcModule.getLocalPadLayout();

   int                  pad       = 0;
   int                  neighbour = 0;
   Vector2D             padPos;
   Vector2D             neighbourPos;

   double               distance;
   double               threshold;

   // Initialize TEve variables

   TEveElementList*     eParent = 0;

   if (leftright == kLEFT) {
      eParent = new TEveElementList("Test getLeftNeighbour");
   } else {
      eParent = new TEveElementList("Test getRightNeighbour");
   }

   //

   for (
      std::vector<int>::const_iterator it = pads.begin();
      it != pads.end();
      ++it
   ) {

      try {

         pad       = *it;

         if (leftright == kLEFT)
            neighbour = layout.getLeftNeighbour(pad);
         else
            neighbour = layout.getRightNeighbour(pad);

         padPos       = layout.getPadCenter(pad);
         neighbourPos = layout.getPadCenter(neighbour);

         // calculate distance between pad and neighbour

         switch (layout.getCoordinateType()) {
            case PadRowLayout2D::CARTESIAN:
               distance = sqrt(
                  pow(padPos[0] - neighbourPos[0], 2)
                  + pow(padPos[1] - neighbourPos[1], 2)
               );

               threshold = 1.5 * layout.getPadWidth(pad);

               break;

            case PadRowLayout2D::POLAR:
               distance = sqrt(
                  pow(padPos[0], 2)
                  + pow(neighbourPos[0], 2)
                  - 2*padPos[0]*neighbourPos[0]*cos(padPos[1] - neighbourPos[1])
               );

               threshold = 1.5 * padPos[0] * layout.getPadWidth(pad);

               break;

            default:
               distance = 1000;
               break;
         }

         if (distance > threshold)
            throw pad;

      } catch (gear::Exception e) {

         int rowNumber = layout.getRowNumber(*it);

         int firstPad = layout.getPadIndex(rowNumber, 0);
         int lastPad  = firstPad + layout.getPadsInRow(rowNumber).size() - 1;

         // Last pad should give an exception when getLeftNeighbour

         if (
            (leftright == kLEFT)
            && (*it == lastPad)
         ) {
            //std::cout << "left: it's ok to have no neighbour." << std::endl;
         }

         // First pad should give an exception when getRightNeighbour

         else if (
            (leftright == kRIGHT)
            && (*it == firstPad)
         ) {
            //std::cout << "right: it's ok to have no neighbour." << std::endl;
         }

         // The middle pads shouldn't give exception, so let's visualize.

         else {
            drawTwoConnectingPads(tpcModule, eParent, pad, neighbour, leftright);
         }

      } catch (int pad) {

         drawTwoConnectingPads(tpcModule, eParent, pad, neighbour, leftright);

      } catch (...) {

         std::cout << "exception on pad " << pad << std::endl;

      }

   }

   gEve->AddElement(eParent);
   gEve->Redraw3D();
}

void VisTPC::testNearestPad(int moduleID) {
   // Initialize variables

   Vector2D                  padPos;
   const TPCModule           &tpcModule   = tpcParameters->getModule(moduleID);
   const PadRowLayout2D      &layout      = tpcModule.getLocalPadLayout();

   const std::vector<double> &planeExtent = layout.getPlaneExtent();
   std::vector<double>       testExtent;

   std::vector<Vector2D>     points;

   // Initialize random generator

   if (!gRandom)
      gRandom = new TRandom(0);

   TRandom& r= *gRandom;
   r.SetSeed(0);

   // Extend the plane extent

   testExtent.resize(4);
   testExtent[0] = planeExtent[0] - (planeExtent[1] - planeExtent[0]) / 10;
   testExtent[1] = planeExtent[1] + (planeExtent[1] - planeExtent[0]) / 10;
   testExtent[2] = planeExtent[2] - (planeExtent[3] - planeExtent[2]) / 10;
   testExtent[3] = planeExtent[3] + (planeExtent[3] - planeExtent[2]) / 10;

   //if (tpcModule.getPadLayoutType() == PadRowLayout2D::POLAR) {
   //   if (testExtent[0] < 0.)
   //      testExtent[0] = 0.;
   //}

   // Constant in one coordinate

   for (int i = 0; i < 25; i++) {
      points.push_back(Vector2D(
         testExtent[0],
         r.Uniform(testExtent[2], testExtent[3])
      ));
   }

   for (int i = 0; i < 100; i++) {
      points.push_back(Vector2D(
         (planeExtent[0] + planeExtent[1]) / 2,
         r.Uniform(planeExtent[2], planeExtent[3])
      ));
   }

   for (int i = 0; i < 25; i++) {
      points.push_back(Vector2D(
         testExtent[1],
         r.Uniform(testExtent[2], testExtent[3])
      ));
   }

   // Constant in the other coordinate

   for (int i = 0; i < 25; i++) {
      points.push_back(Vector2D(
         r.Uniform(testExtent[0], testExtent[1]),
         testExtent[2]
      ));
   }

   for (int i = 0; i < 100; i++) {
      points.push_back(Vector2D(
         r.Uniform(planeExtent[0], planeExtent[1]),
         (planeExtent[2] + planeExtent[3]) / 2
      ));
   }

   for (int i = 0; i < 25; i++) {
      points.push_back(Vector2D(
         r.Uniform(testExtent[0], testExtent[1]),
         testExtent[3]
      ));
   }

   // Random points

   for (int i = 0; i < 100; i++) {
      points.push_back(Vector2D(
         r.Uniform(testExtent[0], testExtent[1]),
         r.Uniform(testExtent[2], testExtent[3])
      ));
   }

   try {
      testNearestPadWithPoints(moduleID, points);
   } catch (std::out_of_range e) {
      std::cout << e.what() << std::endl;
   }
}

void VisTPC::testNearestPadWithPoints(int moduleID, const std::vector<Vector2D>& points) {
   // Initialize variables

   TEveElementList* parent = new TEveElementList("Test getNearestPad (module)");

   const TPCModule           &tpcModule   = tpcParameters->getModule(moduleID);
   const PadRowLayout2D      &layout      = tpcModule.getLocalPadLayout();

   std::stringstream tooltip;

   int padIndex;

   //

   int i = 0;

   for (
      std::vector<Vector2D>::const_iterator it = points.begin();
      it != points.end();
      ++it
   ) {
      ++i;

      padIndex = layout.getNearestPad((*it)[0], (*it)[1]);

      drawNearestPad(tpcModule, (TEveElement*) parent, padIndex, *it, i);
   }

   gEve->AddElement(parent);

   gEve->Redraw3D();
}

void VisTPC::testNearestPadGlobal() {
   // Initialize variables

   const std::vector<double> &planeExtent = tpcParameters->getPlaneExtent();
   std::vector<double>       testExtent;

   std::vector<Vector2D>     points;

   // Initialize random generator

   if (!gRandom)
      gRandom = new TRandom(0);

   TRandom& r= *gRandom;
   r.SetSeed(0);

   // Extend the plane extent

   testExtent.resize(4);
   testExtent[0] = planeExtent[0] - (planeExtent[1] - planeExtent[0]) / 10;
   testExtent[1] = planeExtent[1] + (planeExtent[1] - planeExtent[0]) / 10;
   testExtent[2] = planeExtent[2] - (planeExtent[3] - planeExtent[2]) / 10;
   testExtent[3] = planeExtent[3] + (planeExtent[3] - planeExtent[2]) / 10;

   //if (tpcModule.getPadLayoutType() == PadRowLayout2D::POLAR) {
   //   if (testExtent[0] < 0.)
   //      testExtent[0] = 0.;
   //}

   // Random points

   for (int i = 0; i < 5000; i++) {
      points.push_back(Vector2D(
         r.Uniform(testExtent[0], testExtent[1]),
         r.Uniform(testExtent[2], testExtent[3])
      ));
   }

   try {
      testNearestPadGlobalWithPoints(points);
   } catch (std::out_of_range e) {
      std::cout << e.what() << std::endl;
   }
}

void VisTPC::testNearestPadGlobalWithPoints(const std::vector<Vector2D>& points) {
   // Initialize variables

   TEveElementList* parent = new TEveElementList("Test getNearestPad (global)");

   const TPCModule *tpcModule = 0;
   Vector2D        local;

   //

   int i = 0;

   for (
      std::vector<Vector2D>::const_iterator it = points.begin();
      it != points.end();
      ++it
   ) {
      ++i;

      GlobalPadIndex gpi = tpcParameters->getNearestPad((*it)[0], (*it)[1]);

      tpcModule = & tpcParameters->getModule(gpi.getModuleID());
      local = tpcModule->globalToLocal((*it)[0], (*it)[1]);

      drawNearestPad(*tpcModule, (TEveElement*) parent, gpi.getPadIndex(), local, i);
   }

   gEve->AddElement(parent);

   gEve->Redraw3D();
}

void VisTPC::toGlobal(const TPCModule * aModule, TEveElement* element) {
   TEveTrans*       transformation = NULL;

   const Vector2D*  offset         = &aModule->getOffset();
   double           angle          = aModule->getAngle();
   double           x_pos, y_pos   = 0.;
   float            z_pos          = maxDriftLength;

   switch (aModule->getTPCCoordinateType()) {
      case PadRowLayout2D::CARTESIAN:
         x_pos = (*offset)[0];
         y_pos = (*offset)[1];
       break;
      case PadRowLayout2D::POLAR:
         x_pos = (*offset)[0] * cos((*offset)[1]);
         y_pos = (*offset)[0] * sin((*offset)[1]);
       break;
       default:
         throw "Unknown coordinate type.";
   }

   transformation = element->PtrMainTrans();
   transformation->SetRotByAngles(angle, 0., 0.);
   transformation->Move3PF(x_pos, y_pos, z_pos);

   // Finalize

   transformation = NULL;
}

void  VisTPC::updateGearMgr( gear::GearMgr * aGearMgr )
{
  gearMgr = aGearMgr;

  // set the class to a safe state in case the GearMgr is 0 (default constructor)
  if ( gearMgr == 0 )
  {
    tpcParameters=0;
    innerFieldCageRadius=0;
    outerFieldCageRadius=0;
    maxDriftLength=0;
    return;
  }

  // the GearMgr is valid, set the variables
  tpcParameters = &(gearMgr->getTPCParameters());
  maxDriftLength = tpcParameters->getMaxDriftLength();
  
  const std::vector<double>& extent = tpcParameters->getPlaneExtent();
  switch (tpcParameters->getCoordinateType())
  {
    // for cartesian coordinates we take the edge that is farthest away from 0,0
    case PadRowLayout2D::CARTESIAN:
    {
       // only draw one cylinder, inner radius is 0
       innerFieldCageRadius = 0;
       // test all four corners of the plane extent for max radius
       double rMaxSquare = extent[0]*extent[0] + extent[2]*extent[2];// xMin^2 + yMin^2

       double rSquare = (extent[0]*extent[0] + extent[3]*extent[3]); //check xMin^2 + yMax^2
       if ( rSquare > rMaxSquare )  rMaxSquare = rSquare;

       rSquare = (extent[1]*extent[1] + extent[2]*extent[2]); //check xMax^2 + yMin^2
       if ( rSquare > rMaxSquare )  rMaxSquare = rSquare;

       rSquare = (extent[1]*extent[1] + extent[3]*extent[3]); //check xMax^2 + yMax^2
       if ( rSquare > rMaxSquare )  rMaxSquare = rSquare;
       
       outerFieldCageRadius = sqrt( rMaxSquare );
    }
     break;

     case PadRowLayout2D::POLAR:
       innerFieldCageRadius = extent[0];
       outerFieldCageRadius = extent[1];
       break;

     default:
       throw gear::Exception("Unknown coordinate type in TPCParameters");
   }
 
}