SwimSwimmer.cxx

Go to the documentation of this file.

// $Id: SwimSwimmer.cxx,v 1.22 2007/11/11 04:37:27 rhatcher Exp $
//
// Swimming through a particle, forward or backward
//
// seun@huhepl.harvard.edu
#include "Swimmer/SwimSwimmer.h"
#include "Swimmer/SwimParticle.h"
#include "Swimmer/SwimCondition.h"
#include "Swimmer/SwimG4Stepper.h"
#include "Swimmer/SwimDefStepper.h"
#include "Swimmer/SwimGeo.h"
#include "Swimmer/SwimStepData.h"
#include "BField/BField.h"
#include "Validity/VldRange.h"
#include "MessageService/MsgService.h" // For MSG
#include <string>
#include <cassert>

CVSID("$Id: SwimSwimmer.cxx,v 1.22 2007/11/11 04:37:27 rhatcher Exp $");

//......................................................................
SwimSwimmer::SwimSwimmer(const VldContext& vldc) :
  fMagField(0),
  fStepMax(0.0254*Munits::m),
  fStepMin(1.0e-6*Munits::m),
  fAcc(1.0e-4),
  fNmaxStep(10000)
{
  // This is the "proper" constructor, it should be used for all reco purposes!!!!!!
  // Choice of BField determined from validity. 
  switch (vldc.GetDetector()) {
  case Detector::kFar:
    fMagField  = new BField(vldc);
    fOwnBfield = true;
    break;
  case Detector::kNear:
    fMagField = new BField(vldc);
    fOwnBfield = true;
    break;
  case Detector::kCalDet:
    fMagField = new BField(vldc,0,0);
    fOwnBfield = true;
    break;
  default:
    fOwnBfield = false;
    MSG("Swim",Msg::kError)
      << "SwimSwimmer does not support the "
      << Detector::AsString(vldc.GetDetector())
      << " detector " << endl;
    assert(0);
  }

  fSwimGeo = SwimGeo::Instance(vldc);
  fOwnSwimGeo = false;

  // DefaultStepper depends on fMagBield
  fStepper = this->CreateDefaultStepper();

  fStepData = new SwimStepData();
  fStepData->SetStepper(fStepper);
}

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

SwimSwimmer::SwimSwimmer(BField* magField) :
  fMagField(magField),
  fSwimGeo(0),
  fStepper(this->CreateDefaultStepper()),
  fStepMax(0.0254*Munits::m),   
  fStepMin(1.0e-6*Munits::m),  // Used by NumericalMethods
  fAcc(1.0e-4),                // Used by NumericalMethods
  fNmaxStep(10000),
  fOwnBfield(false),
  fOwnSwimGeo(false)
{
  // This is constructor for special testing ONLY!!!!
  // It allows to force  Swimmer to use some user defined BField. 
  // Don't use this for real reconstruction!!!!! 
  fStepData = new SwimStepData();
  fStepData->SetStepper(fStepper);
}

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

SwimSwimmer::SwimSwimmer(const VldContext& vldc, BField* magField) :
  fMagField(magField),
  fStepMax(0.0254*Munits::m),
  fStepMin(1.0e-6*Munits::m),
  fAcc(1.0e-4),
  fNmaxStep(10000),
  fOwnBfield(false)
{
  // This is constructor for special testing ONLY!!!!
  // It allows to force  Swimmer to use some user defined BField. 
  // Don't use this for real reconstruction!!!!! 
  fSwimGeo = SwimGeo::Instance(vldc);
  fOwnSwimGeo = false;

  // DefaultStepper depends on fMagBield
  fStepper = this->CreateDefaultStepper();
  
  fStepData = new SwimStepData();
  fStepData->SetStepper(fStepper);
}
//......................................................................

SwimSwimmer::~SwimSwimmer()
{
  this->DeleteBfield();
  this->DeleteSwimGeo();
  this->DeleteStepper();

  if (fStepData) {
    delete fStepData;
    fStepData = 0;
  }
}

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

SwimStepper* SwimSwimmer::CreateDefaultStepper()
{
//======================================================================
// Purpose: Create and return a default stepper object
//======================================================================
  return new SwimDefStepper(fMagField,fStepMin,fAcc);
}

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

void SwimSwimmer::DeleteBfield()
{
  if (fMagField && fOwnBfield) {
    delete fMagField;
    fMagField = 0;
  }
}

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

void SwimSwimmer::DeleteSwimGeo()
{
  if (fSwimGeo && fOwnSwimGeo) {
    delete fSwimGeo;
    fSwimGeo = 0;
  }
}

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

void SwimSwimmer::DeleteStepper()
{
  if (fStepper) {
    delete fStepper;
    fStepper = 0;
  }
}

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

bool SwimSwimmer::SetStepper(const char* name)
{
  // Delete any stepper we may have already created
  this->DeleteStepper();

  // Create a new stepper
  if (name == 0) fStepper = this->CreateDefaultStepper();

  string s(name);
  if (s == "")        fStepper = this->CreateDefaultStepper();
  if (s == "default") fStepper = this->CreateDefaultStepper();
  if (s == "G4") fStepper = new SwimG4Stepper(fMagField,fStepMin,fAcc);
  // Provide a template to show people how to add new steppers
  // if (s == "NewStepper") fStepper = new NewStepper();

  if (fStepper) return true;
  // else
  // Sorry, the stepper you asked for does not exist...
  MSG("Swim",Msg::kError) << "Unknown stepper '" << name << "'\n";
  return false;
}

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

bool SwimSwimmer::SwimForward(SwimParticle& particle, SwimCondition& c)
{
  // Arrow of time, true for dt>0, false for dt<0
  fStepData->SetIsForward(true);
  return Swim(particle, c);
}

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

bool SwimSwimmer::SwimBackward(SwimParticle& particle, SwimCondition& c)
{
  // Arrow of time, true for dt>0, false for dt<0
  fStepData->SetIsForward(false);
  return Swim(particle, c);
}

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

bool SwimSwimmer::Swim(SwimParticle& particle, SwimCondition& c)
{
  SwimGeo::SwimMaterial_t material;
  double              stepSize;
  double              distToNextPlane;
  double              pThres = 0.05*Munits::GeV;
  bool                satisfied = false;
  bool                zDirInitial = false, zDir;   // z-momentum direction

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

  // initial layer is undefined (gets set in GetSwimMaterial)
  fStepData->SetSPI(-1);
  
  for (int i=0; i<fNmaxStep; ++i) {
    // Check if condition is satisfied, if yes break out of loop
    satisfied = c.Satisfied(particle);
    if (satisfied) break;

    const TVector3 xyz = particle.GetPosition();
    const TVector3 direction = particle.GetDirection();
    Double_t momZ= particle.GetMomentum().Z();

   // particle traveling in position/negative z-direction
    if ((fStepData->GetIsForward() && momZ>0.0) ||
        (!(fStepData->GetIsForward()) && momZ<0.0))
      zDir = true;
    else
      zDir = false;

    if (i==0)
      zDirInitial = zDir;
    if (zDir!=zDirInitial)
      break;          // particle changed swimming direction. Exit for loop

    // Determine the material of the step
    //material = fSwimGeo->GetSwimMaterial(particle.GetPosition(), zDir);
    material = fSwimGeo->GetSwimMaterial(xyz, zDir, *fStepData);
    
    if (material==SwimGeo::kError)
      break;          // Outside the detector. Exit for loop

    fStepData->SetSwimMaterial(material);

    // Calculate the distance to the next plane
    distToNextPlane = fSwimGeo->DistToNextPlane(xyz,
                                                direction);

    // Far detector limit: 11.6 meters, which always works for near detector
    if (distToNextPlane==0.0 || distToNextPlane>=11.6*Munits::m)
      break;           // Outside the detector. Exit for loop

    //======================================================================
    // Stepsize Calculation
    //======================================================================
    // stepSize changes accordingly if ZCondition is provided
    // stepSize = 1cm if momentum condition is provided
    stepSize = c.StepSize(particle);
  
    // if particle is not traveling vertically
    if (particle.GetDirection().Z() != 0.0) {
      distToNextPlane = fSwimGeo->DistToNextPlane(xyz,
                                                  direction);
      
      // If ((non-zero) distance to the next plane < stepsize)
      if ((distToNextPlane<stepSize) && distToNextPlane!=0.0)
        stepSize = distToNextPlane;
      
      // Set Maximum stepsize
      if (material==SwimGeo::kAir && stepSize>1.0*Munits::m)
        stepSize = 1.0*Munits::m;
      
      if ((material!=SwimGeo::kAir) && (stepSize>fStepMax))
        stepSize = fStepMax;
    }
    else {
      if (material==SwimGeo::kAir) stepSize = 1.0*Munits::m;
      if (material!=SwimGeo::kAir) stepSize = fStepMax;
    } // particle travels vertically

    fStepData->SetStepSize(stepSize);

    //======================================================================
    // Step Action: StepOnce, dEdx
    //======================================================================
    fStepper->Action(particle,fStepData);

    // Check if particle has a momentum < 50MeV
    if (particle.GetMomentumModulus()<pThres) {
      satisfied = c.Satisfied(particle);
      break; // Exit for loop
    }
    
  }
    
  return satisfied;
}

---