SwimGeo.cxx

Go to the documentation of this file.

// $Id: SwimGeo.cxx,v 1.20 2007/11/11 04:37:27 rhatcher Exp $
//
// Roy's code with three additional functions: GetMaterial,
// DistToNextPlane
//
// seun@huhepl.harvard.edu
#include "Swimmer/SwimGeo.h"
#include <iostream>
#include <string>

#include "Swimmer/SwimPlaneInterface.h"
#include "Swimmer/SwimParticle.h"
#include "Swimmer/SwimStepData.h"

#include "Plex/PlexPlaneId.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "UgliGeometry/UgliSteelPlnHandle.h"
#include "MessageService/MsgService.h"

#include "TMath.h"

CVSID("$Id: SwimGeo.cxx,v 1.20 2007/11/11 04:37:27 rhatcher Exp $");
using std::vector;
using std::pair;
static const double kInvSqrt2 = TMath::Sqrt(0.5);
//......................................................................

SwimGeo* SwimGeo::fInstance = 0;

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

SwimGeo *SwimGeo::Instance(const VldContext& vldc)
{
  // Return a compatible instance.

  // If there is an existing one and it is compatible use it
  if (fInstance && (fInstance->fVldRange).IsCompatible(vldc)) return fInstance;

  // What we have is unacceptable
  if (fInstance) {
    MSG("Swim",Msg::kDebug) 
      << "SwimGeo kill version for range " << fInstance->fVldRange << endl;
    delete fInstance;
  }
  fInstance = new SwimGeo(vldc);

  return fInstance;
}

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

SwimGeo::SwimGeo(const VldContext& vldc) :
  fNInterfaces(0),
  fNextPlaneZ(0.0*Munits::m)
{
  // Determine z positions of front/back of steel and scint planes.
  // Ignore planes not in U/V orientation ie. CalDet cosmic counters
  // and FarDet veto shield. But be careful that uninstrumented (bare) 
  // steel planes don't have a view and in general we want those.

  MSG("Swim",Msg::kDebug) << "SwimGeo created for " << vldc << endl;
  fVld=vldc;
  SwimPlaneInterface *spi=0;
  UgliGeomHandle ugh(vldc);
  fVldRange = ugh.GetVldRange();

  std::vector<UgliPlnHandle> uph_vec = ugh.GetPlnHandleVector();

  for (unsigned int i=0; i < uph_vec.size(); ++i) {
    UgliPlnHandle uph = uph_vec[i];
    PlexPlaneId plnid = uph.GetPlexPlaneId();
    PlaneView::PlaneView_t view = plnid.GetPlaneView();
    if (PlaneView::kA == view && PlaneView::kB == view) continue; // cosmics
    if (view > PlaneView::kUnknown) continue;  // veto shield, but not uninstrumented
    SwimGeo::SwimMaterial_t plnmat = 
      (plnid.IsSteel()) ? SwimGeo::kSteel : SwimGeo::kScint;
    Double_t z0 = uph.GetZ0();
    Double_t dz = uph.GetHalfThickness();
    spi = new SwimPlaneInterface(z0-dz,SwimGeo::kAir,plnmat);
    fInterface.Add(spi);
    spi = new SwimPlaneInterface(z0+dz,plnmat,SwimGeo::kAir);
    fInterface.Add(spi);
               
    MSG("Swim",Msg::kDebug) 
      << " add " << plnid.AsString("c") << " " << AsString(plnmat) 
      << " z [" << z0-dz << "," << z0+dz << "]" << endl;
  }

  fInterface.Sort();
  fNInterfaces=fInterface.GetLast()+1;
  if(fNInterfaces>5000)                
    MSG("Swim",Msg::kFatal) << " Interface list not long enough in SwimGeo.h " << endl;
  for(Int_t i=0;i<fNInterfaces;i++){
    fInterfaceList[i]= dynamic_cast<SwimPlaneInterface *>(fInterface.At(i));
  }
  /*
  // Print material and z position
  double prevz = 0.0;
  double inc;
  for (Int_t i=0; i<=fInterface.GetLast(); ++i) {
    SwimPlaneInterface *spi = 
      dynamic_cast<SwimPlaneInterface*>(fInterface.At(i));
    for (int j = 0; j < 10; ++j)
    {
        inc = (spi->fz - prevz)/10;
        cout << (prevz + inc*j) << " " <<  spi->imatbef << endl;
    }
    cout << spi->fz << " " << spi->imatbef << endl;
    prevz = spi->fz;
  }
  */

  // Initialize coil & collar information //

  //far detector
  fFarCoilRadius = 12*Munits::cm;
  fFarOuterCollarRadius = 17*Munits::cm;
  fFarInnerCollarRadius = 13*Munits::cm;
  fFarEndSM1Z = 14.5*Munits::m;
  fFarBegSM2Z = 16*Munits::m; 
  
  // near detector information
  fNearCoilLength = 24*Munits::cm;
  fNearOuterCollarLength = 34*Munits::cm;
  fNearInnerCollarLength = 26*Munits::cm;
  fNearCoil = InitNearCoil();
  fNearInnerCollar = InitNearInnerCollar();
  fNearOuterCollar = InitNearOuterCollar();

  MSG("Swim",Msg::kDebug)<<"Far Coil Radius = "<<fFarCoilRadius<<endl;
  MSG("Swim",Msg::kDebug)<<"Far Inner Coil Radius = "<<fFarInnerCollarRadius<<endl;
  MSG("Swim",Msg::kDebug)<<"Far Outer Collar Radius = "<<fFarOuterCollarRadius<<endl;
  
  if (MsgService::Instance()->IsActive("Swim",Msg::kDebug)) Print();
}

//......................................................................
void SwimGeo::Print(Option_t * /*option*/ ) const
{

  MsgStream& msg = MSGSTREAM("Swim",Msg::kInfo);

  double prevz = 0.0;
  for (int i=0; i <= fInterface.GetLast(); ++i) {
    SwimPlaneInterface *spi =
      dynamic_cast<SwimPlaneInterface*>(fInterface.At(i));
    if ( ! spi ) continue;
    spi->Print("0");
    double thisz = spi->GetZ();
    msg << " dz = " << thisz - prevz << endl;
    prevz = thisz;
  }
}

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

SwimPlaneInterface* SwimGeo::GetInterface(Double_t z, Int_t idir,
                                          SwimGeo::SwimMaterial_t material) const
{
  Int_t iLow=0;
  Int_t iHi=fNInterfaces;
  Int_t iMid=0;
  while( (iHi-iLow)>1){
    iMid=(iLow+iHi) >> 1;
    SwimPlaneInterface *spi = fInterfaceList[iMid];

    if(z>=spi->GetZ())
      iLow=iMid;
    else
      iHi=iMid;
  }
  Int_t ibeg=0;
  Int_t iend=0;
  Int_t idirection=0;
  if (idir==0) {
    ibeg = iHi;
    iend = fNInterfaces;
    idirection = 1;
  }
  else {
    ibeg = iLow;
    iend = -1;
    idirection = -1;
  }
  for (Int_t i=ibeg; i!=iend; i+=idirection) {
    SwimPlaneInterface *spi = fInterfaceList[i];
    Double_t spiZ=spi->GetZ();
    if (spi->GetMaterialBefore() == material || 
        spi->GetMaterialAfter()  == material    ) {
      if (idir==0  && spiZ > z ) {
        return spi;
      }
      else if(idir==1 && spiZ < z){
        return spi;
      }
    }
  }
  MSG("SwimGeo",Msg::kFatal) 
    << "Could not find " << AsString(material) << " interface " 
    << z << " " << idir << " " << ibeg 
    << " " << iend << endl;
  
  return 0;
}

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

SwimGeo::SwimMaterial_t SwimGeo::GetSwimMaterial(const TVector3 xyz, bool zDir)
{

  // This function should not be used now that the coil model
  // has been included in the swimmer. Use the other one

  // Get the material for position xyz
  // If xyz lies on the boundary, get the material after the boundary,
  // which depends on dir

  Int_t iLow=0;
  Int_t iHi=fNInterfaces;
  Int_t iMid=0;
  Double_t xyzZ=xyz.Z();
  while( (iHi-iLow)>1){
    iMid=(iLow+iHi)>>1;
    SwimPlaneInterface *spi = fInterfaceList[iMid];
    if(xyzZ>=spi->GetZ())
      iLow=iMid;
    else
      iHi=iMid;
  }

  Int_t ibeg=0;
  Int_t iend=0;
  Int_t idirection=0;
  if (zDir==1) {
    ibeg = iHi;
    iend = fNInterfaces;
    idirection = 1;
  }
  else {
    ibeg = iLow;
    iend = -1;
    idirection = -1;
  }
 

  for (int i=ibeg; i!=iend; i+=idirection) {
    SwimPlaneInterface *spi = fInterfaceList[i];
    Double_t spiZ=spi->GetZ();
    if (zDir && spiZ > xyz.Z() ) {
      fNextPlaneZ = spi->GetZ();
      switch (spi->GetMaterialBefore() ) {
      case kAir:   return kAir;
      case kScint: return kScint;
      case kSteel: return kSteel;
      default: break;// do nothing
      }
    }
    else if (!zDir && spi->GetZ() < xyz.Z() ) {
      fNextPlaneZ = spi->GetZ();
      switch (spi->GetMaterialAfter() ) {
      case kAir:   return kAir;
      case kScint: return kScint;
      case kSteel: return kSteel;
      default: break; // do nothing
      }
    }
  }
  // Outside the detector
  fNextPlaneZ = 0.0;
  return kError;
}

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

SwimGeo::SwimMaterial_t SwimGeo::GetSwimMaterial(const TVector3 xyz, bool zDir, SwimStepData& stepdata)
{
  
  
  if(SetNextPlaneZ(xyz,zDir,stepdata)) {
    
    // Check if we're inside the coil or collar
    if(IsInsideCoil( xyz , fVld.GetDetector()) &&  
       fVld.GetDetector()==Detector::kNear )return kCoil_near;
    if(IsInsideCoil( xyz ,fVld.GetDetector()) &&  
       fVld.GetDetector()==Detector::kFar )return kCoil_far;
    if(IsInsideOuterCollar( xyz ,fVld.GetDetector()) ) {
      if( IsInsideInnerCollar(xyz,fVld.GetDetector()))  return kAir;
      else return kSteel;
    }
  
    // Get the material for position xyz
    // If xyz lies on the boundary, get the material after the boundary,
    // which depends on dir
    SwimPlaneInterface *spi = fInterfaceList[stepdata.GetSPI()];
    if(zDir) {
      switch (spi->GetMaterialBefore()) {
      case kAir:   return kAir;
      case kScint: return kScint;
      case kSteel: return kSteel;
      default:     break;
      }
    }
    else if(!zDir) {
      switch (spi->GetMaterialAfter()) {
      case kAir:   return kAir;
      case kScint: return kScint;
      case kSteel: return kSteel;
      default:     break;
      }
    }
    return kError;
  }
  
  else return kError;
    
}

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

Bool_t  SwimGeo::SetNextPlaneZ(const TVector3 xyz,bool zDir,SwimStepData& stepdata)
{
  // This function sets fNextPlaneZ,it returns false if out of detector
  int ibeg;
  int iend;
  int idirection;
  int i;
  //  cout << " in get mat " << stepdata.GetSPI() << " " << fNInterfaces << " " << zDir << endl; 
  // if stepdata.fSPI has not been set (is -1), then set it
  if ( stepdata.GetSPI() == -1 ) {
    
    Int_t iLow=0;
    Int_t iHi=fNInterfaces;
    Int_t iMid=0;
    Double_t z = xyz.Z();
    while( (iHi-iLow)>1){
      iMid=(iLow+iHi) >> 1;
      SwimPlaneInterface *spi = fInterfaceList[iMid];
      if(z>=spi->GetZ())
        iLow=iMid;
      else
        iHi=iMid;
    }
    
    if (zDir) {
      stepdata.SetSPI(iHi);
    } else {
      stepdata.SetSPI(iLow);
    }
  }            
 
  ibeg = stepdata.GetSPI();  
  if (zDir) {
    iend = fNInterfaces;
    idirection = 1;
  }
  else {
    iend = -1;
    idirection = -1;
  }
  //  cout << " lookup " << ibeg  << " " << iend << " " << idirection << endl;
  for (i=ibeg; i!=iend; i+=idirection) {
    SwimPlaneInterface *spi = fInterfaceList[i];
    Double_t spiZ=spi->GetZ();
    Double_t xyzZ=xyz.Z();
    //  cout  << " z " << spiZ << " " << xyzZ << endl;
    if (zDir && spiZ > xyzZ ) {
      stepdata.SetSPI(i);
      fNextPlaneZ = spiZ;
      return true;
    }
    else if (!zDir && spiZ < xyzZ ) {
      stepdata.SetSPI(i);
      fNextPlaneZ = spiZ;
      return true;
    }
  }
  
  
  // Outside the detector
  fNextPlaneZ = 0.0;
  return false;
  }
  
//......................................................................

double SwimGeo::DistToNextPlane(const TVector3 xyz, const TVector3 dxyz)
{
// GetMaterial has to be called prior this function to set fNextPlaneZ
//
// Return the distance to the next interface traveling from position xyz
// along direction dxyz
//
// Inputs:
//  xyz - position vector in X,Y,Z coordinates
//  dxyz - normalized direction vector
//
// Returns:
//  distance to next interface

  return TMath::Abs((fNextPlaneZ-xyz.Z())/dxyz.Z());

}
                                                             

vector< pair<Double_t,Double_t> > SwimGeo::InitNearCoil(){

  // coil model for near detector
  vector< pair<Double_t,Double_t> > v;
  Double_t a=fNearInnerCollarLength;
  Double_t b=fNearCoilLength;
   
  // Add points from beam's eye view
  // Assume lower point of coil and inner collar are the same
  
  // Right point
  Double_t rx=-b*kInvSqrt2;
  Double_t ry=-(a-b)*kInvSqrt2;
  // Left point
  Double_t lx=b*kInvSqrt2;
  Double_t ly=-(a-b)*kInvSqrt2;
  // Upper point
  Double_t ux=0;
  Double_t uy=(2*b-a)*kInvSqrt2;
  // Lower point
  Double_t bx=0;
  Double_t by=-a*kInvSqrt2;
  
  v.push_back(pair<Double_t,Double_t>(lx,ly));
  v.push_back(pair<Double_t,Double_t>(ux,uy));   
  v.push_back(pair<Double_t,Double_t>(rx,ry));   
  v.push_back(pair<Double_t,Double_t>(bx,by));   
  v.push_back(pair<Double_t,Double_t>(lx,ly));
  
  MSG("Swim",Msg::kDebug)<<"In Beams Eye View "<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Coil Top(X,Y) = "<<ux<<","<<uy<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Coil Left(X,Y) = "<<lx<<","<<ly<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Coil Bottom(X,Y) = "<<bx<<","<<by<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Coil Right(X,Y) = "<<rx<<","<<ry<<endl;
  return v;
}

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

vector< pair<Double_t,Double_t> > SwimGeo::InitNearOuterCollar(){
  // outer collar model for near detector
  vector< pair<Double_t,Double_t> > v;
  Double_t s=fNearOuterCollarLength*kInvSqrt2;
  
  // Add points from beam's eye view
  // Assume center of collar is center of coordinate system
  v.push_back(pair<Double_t,Double_t>(s,0)); // left
  v.push_back(pair<Double_t,Double_t>(0,s)); // upper   
  v.push_back(pair<Double_t,Double_t>(-s,0));// right  
  v.push_back(pair<Double_t,Double_t>(0,-s));// lower   
  v.push_back(pair<Double_t,Double_t>(s,0));
 
  MSG("Swim",Msg::kDebug)<<"Near Collar Top(X,Y) = "<<0<<","<<s<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Collar Left(X,Y) = "<<s<<","<<0<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Collar Bottom(X,Y) = "<<0<<","<<-s<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Collar Right(X,Y) = "<<-s<<","<<0<<endl;
  
  return v;
}

vector< pair<Double_t,Double_t> > SwimGeo::InitNearInnerCollar(){
  // inner collar model for near detector  
  vector< pair<Double_t,Double_t> > v;
  Double_t s=fNearInnerCollarLength*kInvSqrt2;
  
  // Add points from beam's eye view
  // Assume center of collar is center of coordinate system
  v.push_back(pair<Double_t,Double_t>(s,0)); // left
  v.push_back(pair<Double_t,Double_t>(0,s)); // upper   
  v.push_back(pair<Double_t,Double_t>(-s,0));// right  
  v.push_back(pair<Double_t,Double_t>(0,-s));// lower   
  v.push_back(pair<Double_t,Double_t>(s,0));
    
  MSG("Swim",Msg::kDebug)<<"Near Air Top(X,Y) = "<<0<<","<<s<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Air Left(X,Y) = "<<s<<","<<0<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Air Bottom(X,Y) = "<<0<<","<<-s<<endl;
  MSG("Swim",Msg::kDebug)<<"Near Air Right(X,Y) = "<<-s<<","<<0<<endl;
  
return v;
}

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

 Bool_t SwimGeo::IsInsideCoil(TVector3 xyz, Detector::Detector_t dt)
{
  // Check to see if point is inside coil region
  Double_t x = xyz.X();
  Double_t y = xyz.Y();
  Double_t z = xyz.Z();
  
  Double_t r=TMath::Sqrt(x*x+y*y);
  if ( dt == Detector::kNear ) {
 
    if( IsInside(x,y,fNearCoil) ) return true;
    
  }
  else if ( dt == Detector::kFar ) {
    
    if( r < fFarCoilRadius && z > fFarBegSM2Z && z < fFarEndSM1Z ) return true;
    
  }
  
  return false;
  
}

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

Bool_t SwimGeo::IsInsideOuterCollar(TVector3 xyz, Detector::Detector_t dt)
{
  Double_t x = xyz.X();
  Double_t y = xyz.Y();
  Double_t z = xyz.Z();
  
  Double_t r=TMath::Sqrt(x*x+y*y);
  
  // Check to see if point is inside collar region
  if ( dt == Detector::kNear ) {
    
    if( IsInside(x,y,fNearOuterCollar) ) return true;
    
  }
  else if ( dt == Detector::kFar ) {
   
    if( r < fFarOuterCollarRadius && z > fFarBegSM2Z && z < fFarEndSM1Z ) return true;
    
  }
  
  return false;
}

Bool_t SwimGeo::IsInsideInnerCollar(TVector3 xyz, Detector::Detector_t dt)
{
  Double_t x = xyz.X();
  Double_t y = xyz.Y();
  Double_t z = xyz.Z();
  
  Double_t r=TMath::Sqrt(x*x+y*y);
 
  if ( dt == Detector::kNear ) {
    
    if( IsInside(x,y,fNearInnerCollar) ) return true;
  }
  
  else if ( dt == Detector::kFar ) {
    
    if( r < fFarInnerCollarRadius && z > fFarBegSM2Z && z < fFarEndSM1Z ) return true;
  }
  return false;
}

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

bool SwimGeo::IsInside(Double_t xp, Double_t yp, 
                            const vector< pair<Double_t,Double_t> >& v) const
{
  // Is (xp,yp) inside the plane outline represented by v?
  // taken from TMath::IsInside and modified to use my data structure
  Double_t xint;
  Int_t inter = 0;
  Int_t np=v.size(); 
  for (Int_t i=0;i<np-1;i++) {
    if (v[i].second == v[i+1].second) continue;
    if (yp <= v[i].second && yp <= v[i+1].second) continue;
    if (v[i].second < yp && v[i+1].second < yp) continue;
    xint = v[i].first 
      + (yp-v[i].second)*(v[i+1].first-v[i].first)/(v[i+1].second-v[i].second);
    if (xp < xint) inter++;
  }
  if (inter%2) return true;
  return false;
}

double SwimGeo::GetSwimMaterialDensity(SwimGeo::SwimMaterial_t mat)
{
  switch (mat)  {
  case (SwimGeo::kAir):   return 0.001*Munits::g/(Munits::cm3);
  case (SwimGeo::kScint): return MinosMaterial::Density(MinosMaterial::ePolystyreneMinos)*Munits::g/(Munits::cm3);
  case (SwimGeo::kSteel): return MinosMaterial::Density(MinosMaterial::eIronMinos)*Munits::g/(Munits::cm3);
  case (SwimGeo::kAluminum): return 2.699*Munits::g/(Munits::cm3);
  case (SwimGeo::kCopper): return 8.96*Munits::g/(Munits::cm3);
  case (SwimGeo::kHydrogen): return 0.063*Munits::g/(Munits::cm3);
  case (SwimGeo::kOxygen): return 1.140*Munits::g/(Munits::cm3);
  case (SwimGeo::kWater): return  1.0*Munits::g/(Munits::cm3);
  case (SwimGeo::kCarbon): return  2.265*Munits::g/(Munits::cm3);
  case (SwimGeo::kNitrogen): return  0.808*Munits::g/(Munits::cm3);
  case (SwimGeo::kCoil_near):
    // This does not include a small fraction of
    // insulation or teflon, only Al and H20
    return 2.35731*Munits::g/(Munits::cm3);
    
  case (SwimGeo::kCoil_far): 
    // From Robert's GMINOS code in mix_farcoil.F 
    return 2.3488*Munits::g/(Munits::cm3);
  case (SwimGeo::kError): return 0.0*Munits::g/(Munits::cm3);
  default:                return 0.0*Munits::g/(Munits::cm3);
  };
}  //

---