PisaBBC.C

Go to the documentation of this file.

#include "PisaBBC.h"

#include <algorithm>
#include <cassert>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <vector>

#include "PHCompositeNode.h"
#include "PHIODataNode.h"
#include "PHObjectVectorTCA.h"
#include "phool.h"
#include "PisaBBCHitv1.h"
#include "PisaHelper.h"
#include "PisaPhnxPar.h"

#include "TGeoManager.h"
#include "TGeoMedium.h"
#include "TGeoVolume.h"
#include "TLorentzVector.h"
#include "TVirtualMC.h"
#include "TMCProcess.h"
#include "TArrayI.h"

namespace 
{
  double max(double a, double b, double c)
  {
    return std::max(a,std::max(b,c));
  }

  double max(double a, double b, double c, double d)
  {
    return std::max(a,max(b,c,d));
  }
 
  double min(double a, double b, double c)
  {
    return std::min(a,std::min(b,c));
  }

  double min(double a, double b, double c, double d)
  {
    return std::min(a,min(b,c,d));
  }

  std::string status()
  {
    std::string rv;

    if ( gMC->IsTrackInside() )
      {
        rv += " INSIDE";
      }
    if ( gMC->IsTrackEntering() )
      {
        rv += " ENTERING"; 
      }
    if ( gMC->IsTrackExiting() )
      {
        rv += " EXITING";
      }
    if ( gMC->IsTrackOut() )
      {
        rv += " OUTSIDE";
      }
    if ( gMC->IsTrackAlive() )
      {
        rv += " ALIVE";
      }
    if ( gMC->IsTrackDisappeared() )
      {
        rv += " DISAPPEARED";
      }
    if ( gMC->IsTrackStop() )
      {
        rv += " STOP";
      }
    if ( gMC->IsNewTrack() )
      {
        rv += " NEWTRACK";
      }

    TArrayI array;
    int n = gMC->StepProcesses(array);

    for ( int i = 0; i < n; ++i )
      {
        rv += "+";
        rv += TMCProcessName[array[i]];
      }
    return rv;
  }
}

//_____________________________________________________________________________
PisaBBC::PisaBBC(bool enabled)
  : PisaDetector()
{
  enable(enabled);
}

//_____________________________________________________________________________
PisaBBC::~PisaBBC()
{
}

//_____________________________________________________________________________
void
PisaBBC::beginEvent(PHCompositeNode* topNode)
{
  PisaBBCHitVector* hits = 
    PisaHelper::getClass<PisaBBCHitVector>(topNode, "BBCHits","DST");

  if (!hits)
    {
      fHits = new PHObjectVectorTCA<PisaBBCHit,PisaBBCHitv1>;
      PHIODataNode<PisaBBCHitVector>* node =
        new PHIODataNode<PisaBBCHitVector>(fHits,"BBCHits","PHObject");
      PisaHelper::getPisaNode(topNode,"DST")->addNode(node);
    }
  else
    {
      fHits = hits;
    }
  
}

//_____________________________________________________________________________
void
PisaBBC::createGeometry(PHCompositeNode* topNode)
{

  // This routine defines the geometry for the Beam-Beam
  // Counter. Structure of BBC is place in the mother volume "HALL".
  //
  //     HALL --- BBCM ---                         ! BBC mother volume
  //                     /-- BBCE                  ! Electron absorber
  //                     /-- BBCF                  ! Front panel
  //                     /-- BBCS                  ! Structure tube
  //                     /-- BBCB                  ! Back panel
  //                     /-- BBCC                  ! Thin barrel cover
  //                     /-- BBCD --               ! Detector volume
  //                               /-- BBCA        ! Attachment plate
  //                               /-- BBCQ        ! Quartz block
  //                               /-- BBCP        ! PMT
  //                               /-- BBCR        ! Breeder module
  //                               /-- BBCH        ! Metal shield 
  //

  const int MAXPMT = 100; // Maximum # of PMTs to be installed.
  const int MAXREM =  10; // Maximum # of PMTs to be removed.

  PisaPhnxPar* parameters = 
    PisaHelper::getClass<PisaPhnxPar>(topNode, "PHNXPAR");
  assert(parameters != 0);

  const PisaPhnxSubPar* bbc_par = parameters->get("bbc_par");
  assert(bbc_par != 0);

  TGeoManager* gm = PisaHelper::getGeoManager();

  double BBBCOVER[3];  // BarrelCover with Rmin, Rmax, thick 
  double BBDETECT[10]; // Detector element for PGON
  double BBMOTHER[3];  // BBC Mother Volume with Rmin, Rmax, thick
  double BBSHIELD[10]; // Metal shield PGON

  // Construct needed parameters from the phnxpar parameters.

  BBDETECT[0] = bbc_par->getDouble("BBQUARTZ",0);
  BBDETECT[1] = bbc_par->getDouble("BBQUARTZ",1);
  BBDETECT[2] = bbc_par->getDouble("BBQUARTZ",2);
  BBDETECT[3] = bbc_par->getDouble("BBQUARTZ",3);
  BBDETECT[4] = 
    + bbc_par->getDouble("BBATTACH",4) 
    + bbc_par->getDouble("BBQUARTZ",4) 
    - bbc_par->getDouble("BBPMTSIZ",2) 
    - bbc_par->getDouble("BBBREEDE",2);
  BBDETECT[5] = bbc_par->getDouble("BBQUARTZ",5);
  BBDETECT[6] = max(bbc_par->getDouble("BBATTACH",6),
                    bbc_par->getDouble("BBQUARTZ",6),
                    bbc_par->getDouble("BBPMTSIZ",1));
  BBDETECT[7] = 
    + bbc_par->getDouble("BBATTACH",7) 
    + bbc_par->getDouble("BBQUARTZ",7) 
    + bbc_par->getDouble("BBPMTSIZ",2)  
    + bbc_par->getDouble("BBBREEDE",2);
  BBDETECT[8] = bbc_par->getDouble("BBQUARTZ",8);
  BBDETECT[9] = max(bbc_par->getDouble("BBATTACH",9),
                    bbc_par->getDouble("BBQUARTZ",9),
                    bbc_par->getDouble("BBPMTSIZ",1));

  BBBCOVER[0] = max(bbc_par->getDouble("BBABSORB",1),
                    bbc_par->getDouble("BBFRONTB",1),
                    bbc_par->getDouble("BBSTRUCT",1),
                    bbc_par->getDouble("BBBACKBD",1));
  BBBCOVER[1] = BBBCOVER[0] + bbc_par->getDouble("BBCOVERT");
  BBBCOVER[2] = bbc_par->getDouble("BBSTRUCT",2);
   
  BBMOTHER[0] = min(bbc_par->getDouble("BBABSORB",0),
                    bbc_par->getDouble("BBFRONTB",0),
                    bbc_par->getDouble("BBSTRUCT",0),
                    bbc_par->getDouble("BBBACKBD",0));
  BBMOTHER[1] = BBBCOVER[1];
  BBMOTHER[2] = bbc_par->getDouble("BBSTRUCT",2);

  BBSHIELD[0] = bbc_par->getDouble("BBQUARTZ",0);
  BBSHIELD[1] = bbc_par->getDouble("BBQUARTZ",1);
  BBSHIELD[2] = bbc_par->getDouble("BBQUARTZ",2);
  BBSHIELD[3] = bbc_par->getDouble("BBQUARTZ",3);
  BBSHIELD[4] = 
    + bbc_par->getDouble("BBQUARTZ",4)
    - bbc_par->getDouble("BBPMTSIZ",2)
    - bbc_par->getDouble("BBBREEDE",2);
  BBSHIELD[5] = max(bbc_par->getDouble("BBATTACH",6),
                    bbc_par->getDouble("BBQUARTZ",6),
                    bbc_par->getDouble("BBPMTSIZ",1))
    - bbc_par->getDouble("BBSHITHI");
                    
  BBSHIELD[6] = max(bbc_par->getDouble("BBATTACH",6),
                    bbc_par->getDouble("BBQUARTZ",6),
                    bbc_par->getDouble("BBPMTSIZ",1));
  BBSHIELD[7] =
    + bbc_par->getDouble("BBQUARTZ",7)
    + bbc_par->getDouble("BBPMTSIZ",2)
    + bbc_par->getDouble("BBBREEDE",2); 
  BBSHIELD[8] =  max(bbc_par->getDouble("BBATTACH",9),
                     bbc_par->getDouble("BBQUARTZ",9),
                     bbc_par->getDouble("BBPMTSIZ",1))
    - bbc_par->getDouble("BBSHITHI");
  BBSHIELD[9] = max(bbc_par->getDouble("BBATTACH",9),
                     bbc_par->getDouble("BBQUARTZ",9),
                     bbc_par->getDouble("BBPMTSIZ",1));

  double YSTP = bbc_par->getDouble("BBQUARTZ",6) 
    + bbc_par->getDouble("BBSPACIN")*0.5;
  double XSTP = YSTP*sqrt(3.0);
 
  int MROW = static_cast<int>(floor(bbc_par->getDouble("BBFRONTB",1)/YSTP))+1;
  int MCOL = static_cast<int>(floor(bbc_par->getDouble("BBFRONTB",1)/XSTP))+1;

  double RMAX = bbc_par->getDouble("BBFRONTB",1)
    - BBDETECT[6]*2.0/sqrt(3.0);
  double RMIN = bbc_par->getDouble("BBFRONTB",0)
    + BBDETECT[6]*2.0/sqrt(3.0);

  double RRMA = RMAX*RMAX;
  double RRMI = RMIN*RMIN;

  if ( bbc_par->getInt("BBNUMREM") > MAXREM )
    {
      std::cerr << PHWHERE << "BBC: Number of detector to be removed"
                << " exceeded MAXREM=" << MAXREM
                << std::endl;
      exit(1);
    }

  if ( YSTP <= BBDETECT[6] ) 
    {
      std::cerr << PHWHERE << "BBC : Spacing conflict" << std::endl;
      exit(1);
    }

  // Create a BBC mother volume, BBCM

  TGeoVolume* bbcm = gm->Volume("BBCM","TUBE",
                                bbc_par->getInt("BBMEDMOT"),BBMOTHER,3);

  // Place the mother volume BBCM in the HALL volume
  // BBZPOSIT gives the z position of a quartz radiator front face.

  TGeoVolume* hall = gm->GetVolume("HALL");
  assert(hall != 0);

  double zpos = 
    bbc_par->getDouble("BBZPOSIT",0) 
    + BBMOTHER[2]
    - 2.0*bbc_par->getDouble("BBABSORB",2)
    - 2.0*bbc_par->getDouble("BBFRONTB",2)
    - ( bbc_par->getDouble("BBATTACH",7)
        - bbc_par->getDouble("BBABSORB",2));

  hall->AddNode(bbcm,1,new TGeoTranslation(0.0,0.0,zpos));

  zpos = bbc_par->getDouble("BBZPOSIT",1) 
    - BBMOTHER[2]
    + 2.0*bbc_par->getDouble("BBABSORB",2)
    + 2.0*bbc_par->getDouble("BBFRONTB",2)
    + ( bbc_par->getDouble("BBATTACH",7)
        - bbc_par->getDouble("BBABSORB",2));

  TGeoTranslation t1(0.0,0.0,zpos);
  TGeoRotation r1;

  r1.SetAngles(90.,0.,90.,90.,180.,0.); // mirror symmetry

  TGeoCombiTrans* tr1 = new TGeoCombiTrans(t1,r1);

  hall->AddNode(bbcm,2,tr1);

  // Create BBCE, BBCF, BBCS, BBCB, BBCC, BBCD volumes
  // and place them into BBCM.

  TGeoVolume* bbce = gm->MakeTube("BBCE",
                                  gm->GetMedium(bbc_par->getInt("BBMEDABS")),
                                  bbc_par->getDouble("BBABSORB",0),
                                  bbc_par->getDouble("BBABSORB",1),
                                  bbc_par->getDouble("BBABSORB",2));
  bbcm->AddNode(bbce,1,
                new TGeoTranslation(0.0,0.0,-BBMOTHER[2]+
                                    bbc_par->getDouble("BBABSORB",2)));

  TGeoVolume* bbcf = gm->MakeTube("BBCF",
                                  gm->GetMedium(bbc_par->getInt("BBMEDFRO")),
                                  bbc_par->getDouble("BBFRONTB",0),
                                  bbc_par->getDouble("BBFRONTB",1),
                                  bbc_par->getDouble("BBFRONTB",2));
  bbcm->AddNode(bbcf,1,
                new TGeoTranslation(0.0,0.0,-BBMOTHER[2]+
                                    2.0*bbc_par->getDouble("BBABSORB",2)+
                                    bbc_par->getDouble("BBFRONTB",2)));

  TGeoVolume* bbcs = gm->MakeTube("BBCS",
                                  gm->GetMedium(bbc_par->getInt("BBMEDSTR")),
                                  bbc_par->getDouble("BBSTRUCT",0),
                                  bbc_par->getDouble("BBSTRUCT",1),
                                  bbc_par->getDouble("BBSTRUCT",2));  
  bbcm->AddNode(bbcs,1);

  TGeoVolume* bbcb = gm->MakeTube("BBCB",
                                  gm->GetMedium(bbc_par->getInt("BBMEDBAC")),
                                  bbc_par->getDouble("BBBACKBD",0),
                                  bbc_par->getDouble("BBBACKBD",1),
                                  bbc_par->getDouble("BBBACKBD",2)); 

  bbcm->AddNode(bbcb,1,
                new TGeoTranslation(0.0,0.0,BBMOTHER[2]-
                                    bbc_par->getDouble("BBBACKBD",2)));


  TGeoVolume* bbcc = gm->Volume("BBCC","TUBE",bbc_par->getInt("BBMEDCOV"),
                                BBBCOVER,3);

  bbcm->AddNode(bbcc,1);

  // Now the envelops of the detectors (vacuum)
  
  TGeoVolume* bbcd = gm->Volume("BBCD","PGON",bbc_par->getInt("BBMEDCOV"),
                                BBDETECT,10);
  
  //  bbcm->AddNode(bbcd,1);

  zpos = -BBMOTHER[2] 
    + 2.0*bbc_par->getDouble("BBFRONTB",2)
    + 2.0*bbc_par->getDouble("BBABSORB",2)
    + ( BBDETECT[7]-BBDETECT[4] ) * 0.5;

  int NROW = MROW;
  int NCOL = MCOL;
  int NDET = 0;

  // Find detector element location inside the BBC boundaries.

  double XYPS[2][MAXPMT];
  bool LXYP[MAXPMT];

  while ( NROW >= -MROW ) // Start from top left and go right
    {
      int LROW = abs(NROW/2*2-NROW); // 0 for even, 1 for odd
      double YPOS = NROW*YSTP;
      NCOL = MCOL;
      while ( NCOL >= -MCOL )
        {
          int LCOL = abs(NCOL/2*2-NCOL); // 0 for even, 1 for odd
          if ( LROW == LCOL ) // since elements are staggered
            {
              double XPOS = NCOL*XSTP;
              double RRAD = XPOS*XPOS+YPOS*YPOS;
              if ( (RRAD<=RRMA) && (RRAD>=RRMI) ) // in range ?
                {
                  if ( NDET < MAXPMT )
                    {
                      XYPS[0][NDET] = XPOS;
                      XYPS[1][NDET] = YPOS;
                      LXYP[NDET] = true;
                    }
                  else
                    {
                      std::cerr << PHWHERE << "BBC : max # of PMT exceeded"
                                << std::endl;
                      exit(1);
                    }
                  ++NDET;
                } // in range
            } // even-odd
          --NCOL;
        } // column loop
      --NROW;
    } // row loop

  double RREM=0.0;

  for ( int IP1 = 0; IP1 < bbc_par->getInt("BBNUMREM"); ++IP1 )
    {
      int RD = 1;
      double XX = XYPS[0][NDET] - bbc_par->getDouble("BBPOSREM",IP1*2);
      double YY = XYPS[1][NDET] - bbc_par->getDouble("BBPOSREM",IP1*2+1);
      RREM = XX*XX + YY*YY;
      for ( int IP2 = 1; IP2 < NDET; ++IP2 )
        {
          double XX = XYPS[0][IP2] - bbc_par->getDouble("BBPOSREM",IP1*2);
          double YY = XYPS[1][IP2] - bbc_par->getDouble("BBPOSREM",IP1*2+1);
          double RR = XX*XX + YY*YY;
          if ( RR <= RREM )
            {
              RD = IP2;
              RREM = RR;
            }
        }
      LXYP[RD] = false;
    }

  for ( int IP1 = 0; IP1 < NDET; ++IP1 )
    {
      if ( LXYP[IP1] )
        {
          bbcm->AddNode(bbcd,IP1+1,
                        new TGeoTranslation(XYPS[0][IP1],XYPS[1][IP1],zpos));
        }
    }

  // Now fills the detector empty envelops with usefull stuff =
  // attachment, quartz, PMT, breeder and metal shield.

  double BBATTACH[10];
  for ( size_t i = 0; i < 10; ++i )
    {
      BBATTACH[i] = bbc_par->getDouble("BBATTACH",i);
    }

  TGeoVolume* bbca = gm->Volume("BBCA","PGON",bbc_par->getInt("BBMEDATT"),
                                BBATTACH,10);

  zpos = -(BBDETECT[7]-BBDETECT[4])*0.5
    + (BBATTACH[7]-BBATTACH[4])*0.5;

  bbcd->AddNode(bbca,1,new TGeoTranslation(0.0,0.0,zpos));

  double BBQUARTZ[10];
  for ( size_t i = 0; i < 10; ++i )
    {
      BBQUARTZ[i] = bbc_par->getDouble("BBQUARTZ",i);
    }

  TGeoVolume* bbcq = gm->Volume("BBCQ","PGON",bbc_par->getInt("BBMEDQUA"),
                                BBQUARTZ,10);

  zpos = -(BBDETECT[7]-BBDETECT[4])*0.5
    + (BBQUARTZ[7]-BBQUARTZ[4])*0.5
    + (BBATTACH[7]-BBATTACH[4]);

  bbcd->AddNode(bbcq,1,new TGeoTranslation(0.0,0.0,zpos));

  TGeoVolume* bbcp = gm->MakeTube("BBCP",
                                  gm->GetMedium(bbc_par->getInt("BBMEDPMT")),
                                  bbc_par->getDouble("BBPMTSIZ",0),
                                  bbc_par->getDouble("BBPMTSIZ",1),
                                  bbc_par->getDouble("BBPMTSIZ",2));

  zpos = -(BBDETECT[7]-BBDETECT[4])*0.5
    + (BBQUARTZ[7]-BBQUARTZ[4])
    + (BBATTACH[7]-BBATTACH[4])
    + bbc_par->getDouble("BBPMTSIZ",2);

  bbcd->AddNode(bbcp,1,new TGeoTranslation(0.0,0.0,zpos));
 
  TGeoVolume* bbcr = gm->MakeTube("BBCR",
                                  gm->GetMedium(bbc_par->getInt("BBMEDBRE")),
                                  bbc_par->getDouble("BBBREEDE",0),
                                  bbc_par->getDouble("BBBREEDE",1),
                                  bbc_par->getDouble("BBBREEDE",2));
  
  zpos = -(BBDETECT[7]-BBDETECT[4])*0.5
    + (BBQUARTZ[7]-BBQUARTZ[4])
    + (BBATTACH[7]-BBATTACH[4])
    + bbc_par->getDouble("BBPMTSIZ",2)*2.0
    + bbc_par->getDouble("BBBREEDE",2);

  bbcd->AddNode(bbcr,1,new TGeoTranslation(0.0,0.0,zpos));

  TGeoVolume* bbch = gm->Volume("BBCH","PGON",bbc_par->getInt("BBMEDSHI"),
                                BBSHIELD,10);

  zpos = -(BBDETECT[7]-BBDETECT[4])*0.5
    + (BBATTACH[7]-BBATTACH[4])
    + (BBSHIELD[7]-BBSHIELD[4])*0.5;

  bbcd->AddNode(bbch,1,new TGeoTranslation(0.0,0.0,zpos));  
}

//_____________________________________________________________________________
void
PisaBBC::createMaterials(PHCompositeNode*)
{
  TGeoManager* gm = PisaHelper::getGeoManager();

  int nmat = 16;   // vacuum

  gm->Material("Vacuum",0.0,0.0,0.0,nmat,0.0,0.0);  

  nmat = 9; // aluminium

  gm->Material("Aluminium",26.980,13.000,2.700,nmat);

  nmat = 46; // iron

  gm->Material("Low Mag Iron",55.85,26.0,7.87,nmat,1.76,17.1);

  nmat = 201;

  double aquartz[] = { 28.09 ,    16.0 };
  double zquartz[] = { 14.00 ,     8.0 };
  double wquartz[] = {  0.467,   0.532 };

  double dens = 2.64;

  gm->Mixture("Quartz-BBQ",aquartz,zquartz,dens,2,wquartz,nmat);

  nmat = 202; // why a second material identical to the first one ?

  gm->Mixture("Quartz-PMT",aquartz,zquartz,dens,2,wquartz,nmat);

  nmat = 802;

  // Plastic (lucite) for cover plates (C5H8O2)
  double ap[] = { 12.0,  1.0, 16.0 };
  double zp[] = {  6.0,  1.0,  8.0 };
  double wp[] = {  0.6, 0.08, 0.32 };
  double dp = 1.18;

  gm->Mixture("Plastic",ap,zp,dp,3,wp,nmat);

  // Tracking media # 32 - Vacuum + magnetic field
  int isvol = 0;   // not sensitive
  int ifield = 1;  // magnetic field
  double fieldm = 5.0;  // max field
  double tmaxfd = 5.0;    // maximum angle due to field (one step) in degrees
  double dmaxms = 0.2;    // max disp. due to mulsct. in one step (cm)
  double deemax =  0.2;   // max fractional energy loss in one step
  double epsil =   0.1;  // tracking precision (cm)
  double stmin =   5.00;  // min step due to e loss or mulsct. (cm)

  gm->Medium("Air + Field",32,nmat,isvol,ifield,fieldm,tmaxfd,dmaxms,deemax,epsil,stmin);

  // Tracking media # 26 

  nmat = 9 ; // aluminium
  isvol = 0   ; // insensitive
  ifield = 1  ; // magnetic field
  fieldm = 5.0  ; // max field
  tmaxfd = 1.0   ; // maximum angle due to field (one step) in degrees
  dmaxms = 0.2   ; // max disp. due to mulsct. in one step (cm)
  deemax = 0.5   ; // max fractional energy loss in one step
  epsil = 0.01 ; // tracking precision (cm)
  stmin = 0.01 ; // min step due to e loss or mulsct. (cm)
  
  gm->Medium("Aluminum Frame",26,nmat,isvol,ifield,fieldm,tmaxfd,dmaxms,deemax,epsil,stmin);

  // Tracking medium #201 -- Quartz for BBQ (sensitive detector)
  nmat = 201;
  isvol = 1; // sensitive
  ifield = 1;
  fieldm = 20.0;
  tmaxfd = 45.0;
  dmaxms = 0.2;
  deemax = 0.1;
  epsil = 0.01;
  stmin = 0.1;
  
  gm->Medium("Quartz-BBQ",201,nmat,isvol,ifield,fieldm,tmaxfd,dmaxms,deemax,epsil,stmin);

  nmat = 202;

  gm->Medium("Quartz-PMT",202,nmat,isvol,ifield,fieldm,tmaxfd,dmaxms,deemax,epsil,stmin);

  nmat = 802;
  gm->Medium("Plastic",802,nmat,0,0,0.,10.,1.,0.1,0.01,0.05);

  // Tracking media # 46 - Mu Magnet Yoke (Iron)
  nmat = 46   ; // iron (low threshold)
  isvol = 0   ; // not sensitive
  ifield = 0  ; // no magnetic field
  fieldm = 0.0  ; // max field
  tmaxfd = 0.3    ; // maximum angle due to field (one step) in degrees
  dmaxms = 2.0   ; // max disp. due to mulsct. in one step (cm)
  deemax = 0.2   ; // max fractional energy loss in one step
  epsil = 0.1  ; // tracking precision (cm)
  stmin = 0.1 ; // min step due to e loss or mulsct. (cm)
  gm->Medium("MU non-mag-low",46,nmat,isvol,ifield,fieldm,tmaxfd,dmaxms,deemax,epsil,stmin);

}

//_____________________________________________________________________________
void
PisaBBC::finishEvent(PHCompositeNode*)
{
  std::cout << PHWHERE << " : nhits=" << fHits->size() << std::endl;
  
}

//_____________________________________________________________________________
void
PisaBBC::stepManager()
{
  static unsigned int nevt = 0;

  ++nevt;

  if ( gMC->GetMedium() != 5 )
    { 
      // This should of course be done by the PisaServer...
      // using e.g. a map tracking media -> PisaDetector
      // to know at each step which stepManager to call...
      return;
    }

  if ( gMC->TrackCharge() == 0.0 )
    {
      // Do not deal with e.g. photons.
      return;
    }

  int copy;
  int volid = gMC->CurrentVolID(copy);
  
  int track = gMC->GetStack()->GetCurrentTrackNumber();
  int pid = gMC->TrackPid();
  
  std::cout << " Track # " << std::setw(4) << track
            << " " << gMC->ParticleName(pid)
            << " VolID=" << volid << "(copy=" << copy << ") "
            << " VolName=" << gMC->CurrentVolName()
            << " Medium=" << gMC->GetMedium()
            << status()
            << " Length=" << gMC->TrackLength()
            << std::endl;


  if ( gMC->IsTrackEntering() )
    {
      //      PisaBBCHit* last = fHits->back();
  //     int id = 0;
//       if (!last)
//      {
//        id = last->id()+1;
//      }
      PisaBBCHit* hit = fHits->add(fHits->size());
      //      hit->id(id);
      hit->track(gMC->GetStack()->GetCurrentTrackNumber());
      hit->pid(gMC->TrackPid());
      TLorentzVector pos;
      TLorentzVector mom;
      gMC->TrackPosition(pos);
      gMC->TrackMomentum(mom);
      hit->position(pos.X(),pos.Y(),pos.Z());
      hit->momentum(mom.X(),mom.Y(),mom.Z());
      hit->del(0.0);
      hit->tof(gMC->TrackTime()*1e9);
      hit->len(0.0);
      std::cout << PHWHERE << nevt << " Created hit " << (*hit) << std::endl;
    }

  if ( gMC->IsTrackInside() )
    {
      PisaBBCHit* last = fHits->back();
      assert(last!=0);
      assert(last->track()==gMC->GetStack()->GetCurrentTrackNumber());
      assert(last->pid()==gMC->TrackPid());
      last->del(last->del()+gMC->Edep());
      std::cout << PHWHERE << nevt << " Updating hit " << std::endl;
      std::cout << (*last) << std::endl;
    }

  if ( !gMC->IsTrackAlive() || gMC->IsTrackExiting() || gMC->IsTrackOut() )
    {
      PisaBBCHit* last = fHits->back();
      assert(last!=0);
      assert(last->track()==gMC->GetStack()->GetCurrentTrackNumber());
      assert(last->pid()==gMC->TrackPid());
      last->len(gMC->TrackLength()-last->len());
      std::cout << PHWHERE << nevt << " Really closing hit" << (*last) << std::endl;
    }

 //  if ( gMC->IsTrackEntering() )
//     {
//       TLorentzVector track;
//       gMC->TrackPosition(track);
//       TLorentzVector mom;
//       gMC->TrackMomentum(mom);
//       hit->position(track.X(),track.Y(),track.Z());
//       hit->momentum(mom.X(),mom.Y(),mom.Z());
//       hit->len(gMC->TrackLength());
//       hit->del(0.0);
//       hit->tof(gMC->TrackTime()*1E9);
//       hit->etot(gMC->Etot());
//       hit->track(gMC->GetStack()->GetCurrentTrackNumber());
//       hit->pid(gMC->TrackPid());
//       std::cout << "VolID=" << volid << "(copy=" << copy << ") "
//              << " VolName=" << gMC->CurrentVolName()
//              << " Medium=" << gMC->GetMedium()
//              << " IsExiting=" << gMC->IsTrackExiting()
//              << " IsInside=" << gMC->IsTrackInside()
//              << " IsAlive=" << gMC->IsTrackAlive()
//              << " IsStop=" << gMC->IsTrackStop()
//              << " IsDisappeared=" << gMC->IsTrackDisappeared()
//              << std::endl;
//       hit->identify();
//     }

//   if ( gMC->IsTrackInside() )
//     {
//       std::cout << "Length=" 
//              << gMC->TrackLength()
//              << " VolID=" << volid << "(copy=" << copy << ") "
//              << " VolName=" << gMC->CurrentVolName()
//              << " Medium=" << gMC->GetMedium()
//              << " IsExiting=" << gMC->IsTrackExiting()
//              << " IsInside=" << gMC->IsTrackInside()
//              << " IsAlive=" << gMC->IsTrackAlive()
//              << " IsStop=" << gMC->IsTrackStop()
//              << " IsDisappeared=" << gMC->IsTrackDisappeared()
//              << std::endl;
      
//     }
}