AtNuReco/ClusterCam.cxx

Go to the documentation of this file.

// Package: AtNuReco
//
// ClusterCam.cxx
//
// chapman@hep.phy.cam.ac.uk
// blake@hep.phy.cam.ac.uk
#include "AtNuReco/ClusterCam.h"
#include "AtNuReco/HitCam.h"

ClassImp(ClusterCam)

  //CVSID("$Id: ClusterCam.cxx,v 1.2 2006/04/25 15:59:09 rhatcher Exp $");


ClusterCam::ClusterCam(HitCam* hit) :
  fPlane(-1),fBegStrip(-1), fEndStrip(-1), 
  fBegTime(0.), fEndTime(0.),
  fBegTPos(0.), fEndTPos(0.),
  fZPos(0.), fTPos(0.), fCharge(0.),
  fTrkFlag(0), fShwFlag(0),
  fTrkPlnFlag(0), fShwPlnFlag(0),
  fPlaneView(-1),fDigits(0),fNDFlag(1), 
  StripWidth(4.108e-2)
{
  this->AddHit(hit);
}


ClusterCam::~ClusterCam()
{
  HitsInCluster.clear();
}


void ClusterCam::AddHit(HitCam* hit)
{
  if(HitsInCluster.size()==0) {
    HitsInCluster.push_back(hit);

    fPlane=hit->GetPlane();
    fBegStrip=hit->GetStrip();
    fEndStrip=hit->GetStrip();
    fBegTPos=hit->GetTPos();
    fEndTPos=hit->GetTPos();
    fBegTime=hit->GetTime();
    fEndTime=hit->GetTime();
    fZPos=hit->GetZPos();
    fPlaneView=hit->GetPlaneView();
  }

  else {
    if(this->ContainsHit(hit)==true) {return;}
    HitsInCluster.push_back(hit);   

    if(hit->GetStrip()<fBegStrip) fBegStrip=hit->GetStrip();
    if(hit->GetStrip()>fEndStrip) fEndStrip=hit->GetStrip();
    if(hit->GetTPos()<fBegTPos) fBegTPos=hit->GetTPos();
    if(hit->GetTPos()>fEndTPos) fEndTPos=hit->GetTPos();
    if(hit->GetTime()<fBegTime) fBegTime=hit->GetTime();
    if(hit->GetTime()>fEndTime) fEndTime=hit->GetTime();
  }
  fDigits += hit->GetCandStripHandle()->GetNDaughters();
  fTPos = (fTPos*fCharge+hit->GetTPos()*hit->GetCharge())/(fCharge+hit->GetCharge());
  fCharge += hit->GetCharge();

  return;

}


bool ClusterCam::ContainsHit(HitCam* hit)
{
  for(unsigned int i=0; i<HitsInCluster.size(); ++i) {
    if(hit==HitsInCluster[i]) {return true;}
  }

  return false;
}


int ClusterCam::IsHitAssoc(HitCam* hit) const
{
  double TimeWindow = 9999.9;
  if( hit->GetPlane()==fPlane 
      && hit->GetTPos()>(fBegTPos-2*StripWidth)  && hit->GetTPos()<(fEndTPos+2*StripWidth)
      && hit->GetTime()>(fBegTime-TimeWindow) && hit->GetTime()<(fEndTime+TimeWindow) ) {
    return 1; 
  }
  else {return 0;}
}


int ClusterCam::IsShwAssoc(ClusterCam* clust) const
{
  double TimeWindow = 99.9; int ShwAssocNum = 0;

  if( (fEndTime-clust->GetBegTime())>-TimeWindow || (clust->GetEndTime()-fBegTime)>-TimeWindow ) {

    if( abs(clust->GetPlane()-fPlane)<5 
        && (clust->GetEndTPos()-fBegTPos)>-6*StripWidth && (fEndTPos-clust->GetBegTPos())>-6*StripWidth ) {
      
      if( ( abs(clust->GetPlane()-fPlane)<3 
            && (clust->GetEndTPos()-fBegTPos)>-StripWidth && (fEndTPos-clust->GetBegTPos())>-StripWidth )

          || ( clust->GetPlane()==fPlane
               && (clust->GetEndTPos()-fBegTPos)>-3*StripWidth && (fEndTPos-clust->GetBegTPos())>-3*StripWidth ) )
        
        {ShwAssocNum=2;}
      
      else {ShwAssocNum=1;}
    }

  }

  if(ShwAssocNum==2 && this->GetEntries()<3 && clust->GetEntries()<2) {ShwAssocNum=1;}

  return ShwAssocNum;
}



int ClusterCam::IsTrkAssoc(ClusterCam* clustm, ClusterCam* clustp, int PlaneGap) const
{
  double TimeWindow = 99.9; int TrkAssocNum = 0;
  double min=0.2; double max=0.8;
  double NDScale=1;


  // Configure for correct detector instrumentation
  if(PlaneGap==0) {
    if(this->GetNDFlag()==2) {PlaneGap=10;}
    else {PlaneGap=2;}
  }

  if(PlaneGap==10) {NDScale=2;}



  // Check timing proximity
  if(( (fEndTime-clustm->GetBegTime())>-TimeWindow && (clustm->GetEndTime()-fBegTime)>-TimeWindow)
     && ( (fEndTime-clustp->GetBegTime())>-TimeWindow && (clustp->GetEndTime()-fBegTime)>-TimeWindow) ) {


    // If more than two planes away, scale back width of cluster
    // and then treat as if only two planes away
    double mBegTPos=clustm->GetBegTPos();
    double mEndTPos=clustm->GetEndTPos();
    
    if((fPlane-clustm->GetPlane())>PlaneGap) {
      double mScale = double(PlaneGap)/double(fPlane-clustm->GetPlane());
    
      mBegTPos=fBegTPos+mScale*(mBegTPos-fBegTPos);
      mEndTPos=fEndTPos+mScale*(mEndTPos-fEndTPos);
    }
    
    double pBegTPos=clustp->GetBegTPos();
    double pEndTPos=clustp->GetEndTPos();
    
    if((clustp->GetPlane()-fPlane)>PlaneGap) {
      double pScale = double(PlaneGap)/double(clustp->GetPlane()-fPlane);
      
      pBegTPos=fBegTPos+pScale*(pBegTPos-fBegTPos);
      pEndTPos=fEndTPos+pScale*(pEndTPos-fEndTPos);
    } 
    
    
    // Scale tolerance for matching clusters for cases where there are gaps
    double k0 = 0.5*(clustp->GetPlane()-clustm->GetPlane()-(2*PlaneGap));
    // Make sure we do the same thing for the ND +/- 10 plane sections
    if(PlaneGap==10) {k0/=5;}

    min = min + 0.1*k0; 
    max = max - 0.1*k0; 
        

    // Determine associations
    // For tracks with +ve dtpos/dz
    if( (fEndTPos-clustm->GetBegTPos())>-(NDScale*1.1*StripWidth) && (clustp->GetEndTPos()-fBegTPos)>-(NDScale*1.1*StripWidth) ) {

      // Clusters don't overlap
      if( (fBegTPos-clustm->GetEndTPos())>-(NDScale*0.1*StripWidth) || (clustp->GetBegTPos()-fEndTPos)>-(NDScale*0.1*StripWidth) ) {
        if( TMath::Abs( (clustp->GetBegTPos()-fEndTPos)-(fBegTPos-clustm->GetEndTPos()) )<(NDScale*2.1*StripWidth)
            || ( ((min*mEndTPos)+(max*pEndTPos))>(fBegTPos-(NDScale*0.5*StripWidth)) && ((max*mBegTPos)+(min*pBegTPos))<(fEndTPos+(NDScale*0.5*StripWidth)) ) )
          {TrkAssocNum=2;} 
      }
      
      // Overlapping clusters
      if( (clustm->GetEndTPos()-fBegTPos)>-(NDScale*1.1*StripWidth) && (fEndTPos-clustp->GetBegTPos())>-(NDScale*1.1*StripWidth) ) {
        if(TrkAssocNum<1) TrkAssocNum=1;
      }
    }
    
    
    // For tracks with -ve dtpos/dz
    if( (fBegTPos-clustm->GetEndTPos())<(NDScale*1.1*StripWidth) && (clustp->GetBegTPos()-fEndTPos)<(NDScale*1.1*StripWidth) ) {
      
      // Clusters don't overlap
      if( (fEndTPos-clustm->GetBegTPos())<(NDScale*0.1*StripWidth) || (clustp->GetEndTPos()-fBegTPos)<(NDScale*0.1*StripWidth) ) {
        if( TMath::Abs( (clustp->GetEndTPos()-fBegTPos)-(fEndTPos-clustm->GetBegTPos()) )<(NDScale*2.1*StripWidth)
            || ( ((min*pEndTPos)+(max*mEndTPos))>(fBegTPos-(NDScale*0.5*StripWidth)) && ((max*pBegTPos)+(min*mBegTPos))<(fEndTPos+(NDScale*0.5*StripWidth)) ) )
          {TrkAssocNum=2;}
      }
      
      // Overlapping clusters
      if( (clustm->GetBegTPos()-fEndTPos)<(NDScale*1.1*StripWidth) && (fBegTPos-clustp->GetEndTPos())<(NDScale*1.1*StripWidth) ) {
        if(TrkAssocNum<1) TrkAssocNum=1;
      }
    }
    
  }
  
  return TrkAssocNum;
}


HitCam* ClusterCam::GetHit(unsigned int i) const
{
  if(i<HitsInCluster.size()) {return HitsInCluster[i];}
  else {return 0;}
}


int ClusterCam::IsDiffuseShwAssoc(ClusterCam* clr) const
{
  double win = 99.9;
  int assoc = 0;
  if(fEndTime-clr->GetBegTime()>-win || clr->GetEndTime()-fBegTime>-win){
    if( clr->GetPlane()-fPlane<9 && clr->GetPlane()-fPlane>-9 && clr->GetEndStrip()-fBegStrip>-21 && fEndStrip-clr->GetBegStrip()>-21 )
      {
        if( ( clr->GetPlane()-fPlane<5 && clr->GetPlane()-fPlane>-5 && clr->GetEndStrip()-fBegStrip>-11 && fEndStrip-clr->GetBegStrip()>-11 ) ) assoc=2;
        else assoc=1;
      }
  }
  return assoc;
}

---