ShieldProj.cxx

Go to the documentation of this file.

#include <cmath>
#include <iostream>
#include "TCanvas.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TVector3.h"
#include "TH1.h"
#include "TH2.h"
#include "TSystem.h"
#include "TChain.h"
#include "MessageService/MsgService.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRShieldStrip.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Plex/PlexHandle.h"
#include "Plex/PlexSEIdAltL.h"
#include "Plex/PlexStripEndId.h"
#include "Plex/PlexPlaneId.h"
#include <fstream>
#include "Validity/VldTimeStamp.h"
#include "CandShield/ShieldGeom.h"
#include "ShieldProj.h" 

using std::endl;
using std::cout;
using std::ifstream;
using std::ofstream;

// ******** STATIC DATA MEMBERS **********

CVSID("$Id: ShieldProj.cxx,v 1.5 2007/03/01 17:26:07 rhatcher Exp $");

// **************** CONSTRUCTOR & DESTRUCTOR *********************

ShieldProj::ShieldProj(Double_t vtx0_in, Double_t vtx1_in, Double_t vtx2_in, Double_t vtxCos0_in, Double_t vtxCos1_in, Double_t vtxCos2_in, Int_t plane_in, Int_t strip0_in, ShieldGeom* sg){
  vtx0 = vtx0_in;
  vtx1 = vtx1_in;
  vtx2 = vtx2_in;
  vtxCos0 = vtxCos0_in;
  vtxCos1 = vtxCos1_in;
  vtxCos2 = vtxCos2_in;
  plane = plane_in;
  strip0 = strip0_in;
  
  //Finding intersection and distance
  
  if(sg->PlankExists(plane,strip0)==true){
    planeLineDisFinder(sg->GetPlank_X(plane,strip0),sg->GetPlank_Y(plane,strip0),sg->GetPlank_Z(plane,strip0),sg->GetPlankStrip0_X(plane,strip0),sg->GetPlankStrip0_Y(plane,strip0),sg->GetPlankStrip0_Z(plane,strip0),vtx0,vtx1,vtx2,vtx0+vtxCos0,vtx1+vtxCos1,vtx2+vtxCos2,inter0,inter1,inter2);
    
    if(specialPlank(plane,strip0)==true && inter1<sg->GetPlank_Y(plane,strip0)){ //<--- for special case of "curved" planks (for instance plane 566, plank 3)  
      planeLineDisFinder(sg->GetPlank_X(plane,strip0),sg->GetPlank_Y(plane,strip0),sg->GetPlank_Z(plane,strip0),sg->GetPlank_X(plane-1,15),sg->GetPlank_Y(plane-1,15),sg->GetPlank_Z(plane-1,15),vtx0,vtx1,vtx2,vtx0+vtxCos0,vtx1+vtxCos1,vtx2+vtxCos2,inter0,inter1,inter2);
    }
    
    dpv=distanceFinder(sg->GetPlank_X(plane,strip0),sg->GetPlank_Y(plane,strip0),sg->GetPlank_Z(plane,strip0),inter0,inter1,inter2);
  }
  
  else{
    dpv=-1000;
    inter0=-1000;
    inter1=-1000;
    inter2=-1000;
    MSG("CandShield",Msg::kWarning) << "Plank given to ShieldProj does not exist in shield!" << endl;
  }

  //Projection in plank?
  HitInPlank=false;
  if(dpv<(sg->GetStripsInPlank(plane,strip0))*0.0205 && fabs(inter2-(sg->GetPlank_Z(plane,strip0)))<4.00){
    HitInPlank=true;
  } 
  
}//end of constructor

ShieldProj::~ShieldProj() {
  
}

// ************** PRIVATE METHODS **********************************
// -----------------------------------------------------------------------
// method that tells you wether a plank is "special" or not (special planks
// are the ones that are "twisted", or "curved" (half of them is inclined
// and the other half isn't)

Bool_t ShieldProj::specialPlank(Int_t pl,Int_t pk){
  
  Bool_t verdict;
  
  Int_t specialList[8];
  specialList[0]=566;
  specialList[1]=575;
  specialList[2]=630;
  specialList[3]=639;
  specialList[4]=694;
  specialList[5]=703;
  specialList[6]=758;
  specialList[7]=767;

  verdict=false;
  for(int pi=0;pi<8;pi++){
    if(pl==specialList[pi] && pk == 3){
      verdict = true;
    }
  }
  return verdict;
  
}

//------------------------------------------------------------------------
//Method that calculates the intersection between a plane and a line

void ShieldProj::planeLineDisFinder(Double_t p10,Double_t p11,Double_t p12,Double_t p30,Double_t p31,Double_t p32,Double_t p40, Double_t p41, Double_t p42,Double_t p50,Double_t p51,Double_t p52,Double_t& inter0,Double_t& inter1,Double_t& inter2){ 
  
  Double_t p20=p10;
  Double_t p21=p11;
  Double_t p22=p12+4;
  Double_t det_up;
  Double_t det_do;  
  Double_t t;
  
  det_up=p20*p31*p42-p20*p41*p32-p21*p30*p42+p21*p40*p32+p22*p30*p41-p22*p40*p31-p10*p31*p42+p10*p41*p32+p10*p21*p42-p10*p21*p32-p10*p22*p41+p10*p22*p31+p11*p30*p42-p11*p40*p32-p11*p20*p42+p11*p20*p32+p11*p22*p40-p11*p22*p30-p12*p30*p41+p12*p40*p31+p12*p20*p41-p12*p20*p31-p12*p21*p40+p12*p21*p30;
  
  det_do=-p22*p30*p41-p12*p40*p31+p20*p41*p32-p22*p31*p50+p10*p22*p41+p22*p30*p51+p21*p32*p50+p22*p40*p31+p10*p31*p42+p20*p31*p52+p10*p32*p51-p21*p30*p52-p10*p31*p52+p21*p30*p42+p12*p30*p41-p12*p20*p41-p11*p30*p42+p11*p20*p42-p11*p22*p40+p12*p21*p40-p10*p41*p32-p20*p31*p42-p20*p32*p51-p10*p21*p42-p21*p40*p32+p11*p40*p32+p10*p21*p52-p10*p22*p51+p11*p30*p52-p11*p32*p50-p11*p20*p52+p11*p22*p50-p12*p30*p51+p12*p31*p50+p12*p20*p51-p12*p21*p50;
  
  if(det_do!=0){
    t=-det_up/det_do;
    inter0=p40+(p50-p40)*t;
    inter1=p41+(p51-p41)*t;
    inter2=p42+(p52-p42)*t;
  }
  else{
  inter0=-100;
  inter1=-100;
  inter2=-100;
  MSG("CandShield",Msg::kWarning) << "Determinant in PlaneLineDisFinder was 0. Either bad luck (a lot of) or something fishy!" << endl;
  }
}

// ------------------------------------------------------------------------
// method that calculates the distance between a line and a point

Double_t ShieldProj::distanceFinder(Double_t n0,Double_t n1,Double_t n2,Double_t p0,Double_t p1,Double_t p2){
  
  Double_t x1[3];
  Double_t x2[3];
  Double_t xp[3];
  Double_t v1[3];
  Double_t v2[3];
  Double_t v3[3];
  Double_t magu;
  Double_t magd;

  x1[0]=n0;
  x1[1]=n1;
  x1[2]=n2;
  
  x2[0]=n0;
  x2[1]=n1;
  x2[2]=n2+4;

  xp[0]=p0;
  xp[1]=p1;
  xp[2]=p2;

  for(int ii=0;ii<3;ii++){
    v1[ii]=x2[ii]-x1[ii];
    v2[ii]=x1[ii]-xp[ii];   
  }

  v3[0]=v1[1]*v2[2]-v1[2]*v2[1];
  v3[1]=v1[2]*v2[0]-v1[0]*v2[2];
  v3[2]=v1[0]*v2[1]-v1[1]*v2[0];
  
  magu=sqrt(pow(v3[0],2)+pow(v3[1],2)+pow(v3[2],2));
  magd=sqrt(pow(v1[0],2)+pow(v1[1],2)+pow(v1[2],2));
  if(magd!=0){  
    return magu/magd;
  }
  else{
    return 100;
  }
}

//***************** PUBLIC METHODS ********************************

// -----------------------------------------------------------------
Int_t ShieldProj::GetProjPlane() const {

    return plane;

}

// -----------------------------------------------------------------
Int_t ShieldProj::GetProjStrip0() const {

  return strip0;

}

// ----------------------------------------------------------------
Double_t ShieldProj::GetProjInter_X() const{

  return inter0;

}

// ----------------------------------------------------------------
Double_t ShieldProj::GetProjInter_Y() const{

  return inter1;

}

// ----------------------------------------------------------------
Double_t ShieldProj::GetProjInter_Z() const{

  return inter2;

}

// ---------------------------------------------------------------
Double_t ShieldProj::GetProjDis() const{

  return dpv;

}

// ---------------------------------------------------------------
Bool_t ShieldProj::ProjHitPlank() const{

  return HitInPlank;

}

---