GeoSwimmer.cxx

Go to the documentation of this file.

// GeoSwimmer.cxx
//
// January 20, 2006  M.Ishitsuka First version for GeoSwimmer
//                   ref. SwimSwimmer.cxx
//
// $Id: GeoSwimmer.cxx,v 1.15 2008/08/01 04:23:12 ishi Exp $
//
// Swimming through a particle, forward or backward
//
// seun@huhepl.harvard.edu
#include <string>
#include <cassert>
#include <cmath>

#include <TMCProcess.h>  // MC process codes, e.g. kPPrimary
#include <TVirtualMC.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TInterpreter.h>
#include <TGeoManager.h>
#include <TRandom.h>

#include "GeoSwimmer/GeoSwimZCondition.h"
#include "GeoSwimmer/GeoSwimApplication.h"
#include "GeoSwimmer/GeoSwimParticle.h"
#include "GeoSwimmer/GeoSwimStack.h"
#include "Util/LoadMinosPDG.h"
#include "MCApplication/MCApplication.h"
#include "GeoSwimmer/GeoSwimmer.h"
#include "Validity/VldContext.h"
#include "MessageService/MsgService.h" // For MSG
#include "Conventions/Munits.h"
#include "Conventions/Detector.h"

#include "Swimmer/SwimParticle.h"

#if defined(MACOSX) && (__GNUC__ <= 3)
  extern "C" int isnan(double);
  extern "C" int isinf(double);
#endif

CVSID("$Id: GeoSwimmer.cxx,v 1.15 2008/08/01 04:23:12 ishi Exp $");

ClassImp(GeoSwimmer)
  
  GeoSwimmer* GeoSwimmer::fgInstance = 0;

// Constructor
GeoSwimmer::GeoSwimmer() :
  fMCConfig("PTSim_g3Config.C"), fHallSizeSet(false)
{
  fMCAppUser = new GeoSwimApplication("GeoSwimApplication","GeoSwimmer Application");
  fUseGeoSwimmer = false;
  fgInstance = this;

  LoadMinosPDG(); // initialize MINOS specific particles
}

// Destructor
GeoSwimmer::~GeoSwimmer()
{
  if ( fMCAppUser ) delete fMCAppUser; fMCAppUser = 0;
}

//......................................................................

void GeoSwimmer::Initialize(const VldContext& vldc)
{
  switch (vldc.GetDetector()) {
  case Detector::kFar:
    break;
  case Detector::kNear:
    break;
  case Detector::kCalDet:
    break;
  default:
    MSG("GeoSwim",Msg::kError)
      << "GeoSwimmer does not support the "
      << Detector::AsString(vldc.GetDetector())
      << " detector " << endl;
    assert(0);
  }

  MSG("GeoSwim",Msg::kDebug) << "Initialize." << endl;
  if ( !fMCAppUser ) {
    MSG("GeoSwim",Msg::kWarning) << "No MC application instance." << endl;
    abort();
  }

  MCApplication& mcapp = const_cast<MCApplication&>(MCApplication::Instance());
  if ( !(mcapp.GetMC() ) ) {
    mcapp.InitMC(fMCConfig.c_str(),vldc);
  }
  mcapp.InitSnarl(vldc);
  
  // Tell the MCApplication Instance about the current User application
  mcapp.SetUserApplication(fMCAppUser);

  if ( !fHallSizeSet ) { 
    UgliGeomHandle ugh(vldc);
    fHallMin = ugh.GetHallExtentMin();
    fHallMax = ugh.GetHallExtentMax();
    fMCAppUser->SetHallSize(fHallMin,fHallMax); 
    fHallSizeSet = true;
  }
}

//......................................................................

bool GeoSwimmer::SwimForward(GeoSwimParticle& particle, double zFinal)
{
  MSG("GeoSwim",Msg::kDebug) << "GeoSwimmer::SwimForward\n";

  GeoSwimZCondition zc(zFinal);

  return Swim(particle, zc, 1);
}

//......................................................................

bool GeoSwimmer::SwimBackward(GeoSwimParticle& particle, double zFinal)
{
  MSG("GeoSwim",Msg::kDebug) << "GeoSwimmer::SwimBackward\n";

  GeoSwimZCondition zc(zFinal);

  return Swim(particle, zc, -1);
}

//......................................................................

bool GeoSwimmer::Swim(GeoSwimParticle& particle, GeoSwimZCondition& zc, int dir)
{

  double              pThres = 0.05*Munits::GeV;
  double              fUnreasonableEnergy = 1.e+6*Munits::GeV;
  bool                satisfied = false;

  fMCAppUser->SetDirection(dir); 
  fMCAppUser->SetZCondition(zc); 

  // Check if particle has a momentum < 50MeV
  if (particle.GetMomentumModulus()<pThres) {
    satisfied = zc.Satisfied(particle);
    return satisfied;
  }

  Int_t toBeDone = 1;  // mark as primary
  Int_t parent = -1; // -1 => no parent track 
  TVector3 pol;
  Int_t trkId= -1; // filled by stack

  TVector3 momentum = particle.GetMomentum();
  Double_t energy = particle.GetEnergy();
  TVector3 position = particle.GetPosition();  

  if (position.X() > fHallMax.X() ||
      position.X() < fHallMin.X() ||
      position.Y() > fHallMax.Y() ||
      position.Y() < fHallMin.Y() ||
      position.Z() > fHallMax.Z() ||
      position.Z() < fHallMin.Z()) {
    MSG("GeoSwim",Msg::kDebug) 
      << "GeoSwimmer::Swim too large position (m) "
      << position.x() << " " << position.y() << " " << position.z() << "\n";
    return false;
  }

  if (isnan(momentum.X()) || isinf(momentum.X()) || 
      isnan(momentum.Y()) || isinf(momentum.Y()) || 
      isnan(momentum.Z()) || isinf(momentum.Z())) {
    MSG("GeoSwim",Msg::kError) 
      << "Momentum is not number\n";
    return false;
  }

  if (TMath::Abs(momentum.X()) > fUnreasonableEnergy ||
      TMath::Abs(momentum.Y()) > fUnreasonableEnergy ||
      TMath::Abs(momentum.Z()) > fUnreasonableEnergy) {
    MSG("GeoSwim",Msg::kDebug) 
      << "GeoSwimmer::Swim too large momrutum (GeV) "
      << momentum.X() << " " << momentum.Y() << " " << momentum.Z() << "\n";
    return false;
  }

  MSG("GeoSwim",Msg::kDebug)
    << "GeoSwimmer::Swim Begin Momentum "
    << momentum.x() << " " << momentum.y() << " " << momentum.z() << " "
    << " Position(m) " << position.x() << " " << position.y() << " "
    << position.z() << " " << dir << "\n";
  
  if(position.Mag()>10000.) return satisfied; // Temporary to avoid error.

  gMC -> GetRandom() -> SetSeed(0);

  GeoSwimStack* fStack = fMCAppUser -> GetStack();
  fStack -> Reset();  

  // Add particle to stack 
  fStack->PushTrack(toBeDone,parent, dir*particle.GetID(),
                    dir*momentum.x(), dir*momentum.y(), dir*momentum.z(), 
                    energy, 
                    position.x()/Munits::cm, position.y()/Munits::cm, position.z()/Munits::cm, 
                    1.0,
                    pol.X(), pol.Y(), pol.Z(),
                    kPPrimary,  // sets bit in UniqueId
                    trkId,1.,0);  
 
  // Running geant
  gMC -> ProcessRun(1);

  // Collect results
  satisfied = fMCAppUser->GetSatisfied(); 

  TLorentzVector pos4;
  gMC->TrackPosition(pos4);
  TVector3 pos3(pos4.X()*Munits::cm,pos4.Y()*Munits::cm,pos4.Z()*Munits::cm);
  particle.SetPosition(pos3);

  TLorentzVector mom4;
  gMC->TrackMomentum(mom4);
  TVector3 mom3(dir*mom4.X(),dir*mom4.Y(),dir*mom4.Z());
  particle.SetMomentum(mom3);

  particle.SetS(fMCAppUser->GetS()*Munits::cm);
  particle.SetRange(fMCAppUser->GetRange());

  MSG("GeoSwim",Msg::kDebug)
    << "GeoSwimmer::Swim End Momentum "
    << mom3.X() << " " << mom3.Y() << " " << mom3.Z() << " "
    << " Positiopn (m) " << pos3.X() << " " << pos3.Y() << " "
    << pos3.Z() << " " << satisfied << "\n";
  
  return satisfied;
}

// SwimParticle is alloewed as input (temporary)

bool GeoSwimmer::SwimForward(SwimParticle& swpt, double zFinal)
{
  MSG("GeoSwim",Msg::kDebug) << "GeoSwimmer::SwimForward\n";

  GeoSwimZCondition zc(zFinal);
  GeoSwimParticle particle(swpt.GetPosition(), swpt.GetMomentum(),
                           swpt.GetCharge(), (int)swpt.GetCharge()*-1*13);

  bool satisfied = Swim(particle, zc, 1);

  swpt.SetPosition(particle.GetPosition());
  swpt.SetMomentum(particle.GetMomentum());
  swpt.AddS(particle.GetS());
  swpt.AddRange(particle.GetRange());

  return satisfied;
}

//......................................................................

bool GeoSwimmer::SwimBackward(SwimParticle& swpt, double zFinal)
{
  MSG("GeoSwim",Msg::kDebug) << "GeoSwimmer::SwimBackward\n";

  GeoSwimZCondition zc(zFinal);
  GeoSwimParticle particle(swpt.GetPosition(), swpt.GetMomentum(),
                           swpt.GetCharge(), (int)swpt.GetCharge()*-1*13);

  bool satisfied = Swim(particle, zc, -1);

  swpt.SetPosition(particle.GetPosition());
  swpt.SetMomentum(particle.GetMomentum());
  swpt.AddS(particle.GetS());
  swpt.AddRange(particle.GetRange());

  return satisfied;
}

---