TrackSegmentCamAtNu.cxx

Go to the documentation of this file.

// Package: AtNuReco
//
// TrackSegmentCamAtNu.cxx
//
// chapman@hep.phy.cam.ac.uk
// blake@hep.phy.cam.ac.uk
#include "AtNuReco/TrackSegmentCamAtNu.h"
#include "AtNuReco/ClusterCamAtNu.h"
#include "AtNuReco/HitCamAtNu.h"
#include "math.h"

ClassImp(TrackSegmentCamAtNu)
 
CVSID("$Id: TrackSegmentCamAtNu.cxx,v 1.1 2009/01/25 16:55:49 blake Exp $");


TrackSegmentCamAtNu::TrackSegmentCamAtNu(ClusterCamAtNu* clustm, ClusterCamAtNu* clust0, ClusterCamAtNu* clustp):
  fSeedSegment(0), fPartner(0), fUID(0), fBegPlane(999), fEndPlane(-999), fBegVtxZ(999.), 
  fEndVtxZ(-999.), fTrkFlag(0), fTmpTrkFlag(0), fPlaneView(-1), fNPlanes(0), 
  fBegTime(0.), fEndTime(0.), StripWidth(4.108e-2)
{
  fPlaneView=clust0->GetPlaneView();
  fBegTime=clust0->GetBegTime();
  fEndTime=clust0->GetEndTime();

  this->AddCluster(clustm); this->AddCluster(clust0); this->AddCluster(clustp);
}


TrackSegmentCamAtNu::~TrackSegmentCamAtNu()
{
  ClustersInSegment.clear();
  fBegAssociatedSegList.clear();
  fEndAssociatedSegList.clear();
  fBegPreferredSegList.clear();
  fEndPreferredSegList.clear();
  fBegMatchedSegList.clear();
  fEndMatchedSegList.clear();
}


void TrackSegmentCamAtNu::AddCluster(ClusterCamAtNu* clust)
{
  if(clust) {
    if(this->ContainsCluster(clust)==false) 
      {
        ClustersInSegment.push_back(clust);
        if( fBegPlane > clust->GetPlane()){ fBegPlane = clust->GetPlane(); fBegVtxZ =  clust->GetZPos(); }
        if( fEndPlane < clust->GetPlane()){ fEndPlane = clust->GetPlane(); fEndVtxZ =  clust->GetZPos(); }
        
        if( clust->GetBegTime()<fBegTime ) { fBegTime=clust->GetBegTime(); }
        if( clust->GetEndTime()>fEndTime ) { fEndTime=clust->GetEndTime(); }
      }
  }
  
  return;
}


bool TrackSegmentCamAtNu::ContainsCluster(ClusterCamAtNu* clust)
{
  for(unsigned int i=0; i<ClustersInSegment.size(); ++i) {
    if(clust==ClustersInSegment[i]) {return true;}
  }

  return false;
}


ClusterCamAtNu* TrackSegmentCamAtNu::GetCluster(unsigned int i)
{
  if(i<ClustersInSegment.size()) {return ClustersInSegment[i];}
  
  else {return 0;}
}


void TrackSegmentCamAtNu::AddSegment(TrackSegmentCamAtNu* segment)
{
  for(unsigned int i=0; i<segment->GetEntries(); ++i) {
    this->AddCluster(segment->GetCluster(i));    
  }

  return;
}


unsigned int TrackSegmentCamAtNu::GetEntries() const
{
  return ClustersInSegment.size();
}


int TrackSegmentCamAtNu::GetBegPlane() const
{
  return fBegPlane;
}


int TrackSegmentCamAtNu::GetEndPlane() const
{
  return fEndPlane;
}


bool TrackSegmentCamAtNu::IsAssoc(TrackSegmentCamAtNu* segment)
{
  unsigned int i;
  bool assoc=false;
  bool flag=false; 
  TrackSegmentCamAtNu* Seg1 = this;
  TrackSegmentCamAtNu* Seg2 = segment;

  //If the two clusters overlap at all
  if(Seg1->GetEndPlane()>=Seg2->GetBegPlane()){
    flag=true;
    //All clusters in Seg1 in planes overlapping Seg2 should also be in Seg2
    for(i=0;i<Seg1->GetEntries();i++){
      ClusterCamAtNu* clr = Seg1->GetCluster(i);
      if(clr->GetPlane()>=Seg2->GetBegPlane()){ 
        if(!(Seg2->ContainsCluster(clr))){
          flag=false; 
          break; //one mismatch is too many
        }
      }
    }
    if(flag)//only bother with 2nd check if 1st was ok
      {
        //All clusters in Seg2 in planes overlapping Seg1 should also be in Seg1
        for(i=0;i<Seg2->GetEntries();i++){
          ClusterCamAtNu* clr = Seg2->GetCluster(i);
          if(clr->GetPlane()<=Seg1->GetEndPlane()){ 
            if(!(Seg1->ContainsCluster(clr))){
              flag=false; 
              break; //one mismatch is too many 
            }
          }
        }
      }
    //If this is true and they overlap by more than one plane then they are associated
    if(Seg1->GetEndPlane()>Seg2->GetBegPlane()){ if(flag) assoc=true; }
  }
    //If the segments overlap by 1 or fewer planes need to do some more work...    
  if(Seg1->GetEndPlane()<=Seg2->GetBegPlane()){
    int idb=0,idb0=0;
    int bpln=Seg2->GetEndPlane()+1;
    //find the first (lowest plane no) two clusters in the second segment
    for(i=0;i<Seg2->GetEntries();i++){
      ClusterCamAtNu* clr = Seg2->GetCluster(i);
      if(clr->GetPlane()<bpln && clr->GetPlane()>Seg2->GetBegPlane()+2){
        bpln=clr->GetPlane(); idb=i;
      }
      else{ if(clr->GetPlane()==Seg2->GetBegPlane()) idb0=i; }
    }
    ClusterCamAtNu* clrb = Seg2->GetCluster(idb); ClusterCamAtNu* clrb0 = Seg2->GetCluster(idb0);
    double bBegTPos = clrb->GetBegTPos(); double bEndTPos = clrb->GetEndTPos();
    double b0BegTPos = clrb0->GetBegTPos(); double b0EndTPos = clrb0->GetEndTPos();

    int ide=0,ide0=0;
    int epln=Seg1->GetBegPlane()-1;
    //find the last (highest plane no) two clusters in the first segment
    for(i=0;i<Seg1->GetEntries();i++){
      ClusterCamAtNu* clr = Seg1->GetCluster(i);
      if(clr->GetPlane()>epln && clr->GetPlane()<Seg1->GetEndPlane()-2){
        epln=clr->GetPlane(); ide=i;
      }
      else{ if(clr->GetPlane()==Seg1->GetEndPlane()) ide0=i; }
    }
    ClusterCamAtNu* clre = Seg1->GetCluster(ide); ClusterCamAtNu* clre0 = Seg1->GetCluster(ide0);
    double eBegTPos = clre->GetBegTPos(); double eEndTPos = clre->GetEndTPos();
    double e0BegTPos = clre0->GetBegTPos(); double e0EndTPos = clre0->GetEndTPos();
    //Look at how these clusters overlap stripwise    
    if( ( bEndTPos-b0BegTPos>-2*StripWidth && b0EndTPos-bBegTPos>-2*StripWidth && eEndTPos-e0BegTPos>-2*StripWidth && e0EndTPos-eBegTPos>-2*StripWidth )
        || ( ( bBegTPos-b0BegTPos>-2*StripWidth && e0BegTPos-eBegTPos>-2*StripWidth ) || ( bEndTPos-b0EndTPos>-2*StripWidth && e0EndTPos-eEndTPos>-2*StripWidth ) )
        || ( ( bBegTPos-b0BegTPos<2*StripWidth && e0BegTPos-eBegTPos<2*StripWidth ) || ( bEndTPos-b0EndTPos<2*StripWidth && e0EndTPos-eEndTPos<2*StripWidth ) ) 
        ){
      //if the segments overlap by 1 plane then we can now say if they are associated or not
      if(Seg2->GetBegPlane()==Seg1->GetEndPlane()){ if(flag) assoc=true; }
      else{
        double z1,z2,t1,t2,dt1,dt2,dir1,dir2,dirtmp,win;
        z1=Seg1->GetEndZPos();  z2=Seg2->GetBegZPos();
        t1=Seg1->GetEndTPos();  t2=Seg2->GetBegTPos();
        dir1=Seg1->GetEndDir(); dt1=t2-t1-dir1*(z2-z1); 
        dir2=Seg2->GetBegDir(); dt2=t2-t1-dir2*(z2-z1); 
        win=0.1+0.35*(z2-z1);   dirtmp=(1+dir1*dir2)/(sqrt(1+dir1*dir1)*sqrt(1+dir2*dir2));
        
        if( ( (dt1>-win&&dt1<win)||(dt2>-win&&dt2<win) ) && ( dirtmp>0.65 ) ) assoc=true; 
      }
    }
  }
  return assoc;
}



double TrackSegmentCamAtNu::GetBegTPos()
{
  double tot=0.0,begt=0.0;
  //Loop over clusters in segment
  unsigned int nclusters = this->GetEntries();
  unsigned int nhits = 0;
  HitCamAtNu* hit=0;
  for( unsigned int i=0; i<nclusters; ++i)
    {//find the clusters on the first plane in the segment
      if(ClustersInSegment[i]->GetPlane()==fBegPlane)
        {
          nhits = ClustersInSegment[i]->GetEntries();

          //loop over hits in cluster
          for(unsigned int j=0;j<nhits;++j)
            {
              hit = ClustersInSegment[i]->GetHit(j);
              if(hit)
                {
                  begt+=hit->GetTPos();
                  tot+=1.0;
                }
            }
        }
    }
  if(tot>0) return (begt/tot); 
  else return 0;
}



double TrackSegmentCamAtNu::GetEndTPos()
{
  double tot=0.0,endt=0.0;
  //Loop over clusters in segment
  unsigned int nclusters = this->GetEntries();
  unsigned int nhits = 0;
  HitCamAtNu* hit=0;
  for( unsigned int i=0; i<nclusters; ++i)
    {//find the clusters on the first plane in the segment
      if(ClustersInSegment[i]->GetPlane()==fEndPlane)
        {
          nhits = ClustersInSegment[i]->GetEntries();

          //loop over hits in cluster
          for(unsigned int j=0;j<nhits;++j)
            {
              hit = ClustersInSegment[i]->GetHit(j);
              if(hit)
                {
                  endt+=hit->GetTPos();
                  tot+=1.0;
                }
            }
        }
    }
  if(tot>0) return (endt/tot); 
  else return 0;
}



double TrackSegmentCamAtNu::GetBegDir()
{
  double z=0.0,t=0.0;
  double sw=0.0,swx=0.0,swx2=0.0;
  double swy=0.0,swyx=0.0;
  unsigned int nclusters = this->GetEntries();
  unsigned int nhits = 0;
  HitCamAtNu* hit=0;
  for( unsigned int i=0; i<nclusters; ++i)
    {//find the clusters on the first plane in the segment
      if(ClustersInSegment[i]->GetPlane()<fBegPlane+10)
        {
          nhits = ClustersInSegment[i]->GetEntries();

          //loop over hits in cluster
          for(unsigned int j=0;j<nhits;++j)
            {
              hit = ClustersInSegment[i]->GetHit(j);
              if(hit)
                {
                  z=hit->GetZPos(); t=hit->GetTPos();
                  sw+=1.0; swx+=z; swx2+=z*z; 
                  swy+=t; swyx+=t*z;;
                }
            }
        }
    }
  if((swx*swx-sw*swx2)!=0) {return (swx*swy-sw*swyx)/(swx*swx-sw*swx2);}
  else return 0;
}


double TrackSegmentCamAtNu::GetEndDir()
{
  double z=0.0,t=0.0;
  double sw=0.0,swx=0.0,swx2=0.0;
  double swy=0.0,swyx=0.0;
  unsigned int nclusters = this->GetEntries();
  unsigned int nhits = 0;
  HitCamAtNu* hit=0;
  for( unsigned int i=0; i<nclusters; ++i)
    {//find the clusters on the first plane in the segment
      if(ClustersInSegment[i]->GetPlane()>fEndPlane-10)
        {
          nhits = ClustersInSegment[i]->GetEntries();

          //loop over hits in cluster
          for(unsigned int j=0;j<nhits;++j)
            {
              hit = ClustersInSegment[i]->GetHit(j);
              if(hit)
                {
                  z=hit->GetZPos(); t=hit->GetTPos();
                  sw+=1.0; swx+=z; swx2+=z*z; 
                  swy+=t; swyx+=t*z;;
                }
            }
        }
    }
  
  if((swx*swx-sw*swx2)!=0) {return (swx*swy-sw*swyx)/(swx*swx-sw*swx2);}
  else return 0;
}


double TrackSegmentCamAtNu::GetBegZPos() const
{
  return fBegVtxZ;
}


double TrackSegmentCamAtNu::GetEndZPos() const
{
  return fEndVtxZ;
}


void TrackSegmentCamAtNu::AddAssocSegToBeg(TrackSegmentCamAtNu* seg) //void AddTriToBeg(TrackSegmentCamAtNu* seg);
{
  fBegAssociatedSegList.push_back(seg);
  return;
}


void TrackSegmentCamAtNu::AddAssocSegToEnd(TrackSegmentCamAtNu* seg) //void AddTriToEnd(TrackSegmentCamAtNu* seg);
{
  fEndAssociatedSegList.push_back(seg);
  return;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetAssocSegBeg(unsigned int i)//TrackSegmentCamAtNu* GetTriBeg(unsigned int i);
{
  if(i<fBegAssociatedSegList.size()) {return fBegAssociatedSegList[i];}
  else return 0;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetAssocSegEnd(unsigned int i)//TrackSegmentCamAtNu* GetTriEnd(unsigned int i);
{
  if(i<fEndAssociatedSegList.size()) {return fEndAssociatedSegList[i];}
  else return 0;
}


void TrackSegmentCamAtNu::AddPrefSegToBeg(TrackSegmentCamAtNu* seg)//void AddAssocTriToBeg(TrackSegmentAtNu* segment);
{
  fBegPreferredSegList.push_back(seg);
  return;
}


void TrackSegmentCamAtNu::AddPrefSegToEnd(TrackSegmentCamAtNu* seg)//void AddAssocTriToEnd(TrackSegmentAtNu* segment);
{
  fEndPreferredSegList.push_back(seg);
  return;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetPrefSegBeg(unsigned int i)//TrackSegmentAtNu* GetAssocTriBegAt(Int_t i);
{
  if(i<fBegPreferredSegList.size()) {return fBegPreferredSegList[i];}
  else return 0;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetPrefSegEnd(unsigned int i)//TrackSegmentAtNu* GetAssocTriEndAt(Int_t i);
{
  if(i<fEndPreferredSegList.size())  {return fEndPreferredSegList[i];}
  else return 0;
}


void TrackSegmentCamAtNu::AddMatchSegToBeg(TrackSegmentCamAtNu* seg)//void AddSegToBeg(TrackSegmentAtNu* segment);
{
  fBegMatchedSegList.push_back(seg);
  return;
}


void TrackSegmentCamAtNu::AddMatchSegToEnd(TrackSegmentCamAtNu* seg)//void AddSegToEnd(TrackSegmentAtNu* segment)
{
  fEndMatchedSegList.push_back(seg);
  return;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetMatchSegBeg(unsigned int i)//TrackSegmentAtNu* GetSegBegAt(Int_t i)
{
  if(i<fBegMatchedSegList.size()) {return fBegMatchedSegList[i];}
  else return 0;
}


TrackSegmentCamAtNu* TrackSegmentCamAtNu::GetMatchSegEnd(unsigned int i)//TrackSegmentAtNu* GetSegEndAt(Int_t i)
{
  if(i<fEndMatchedSegList.size()) {return fEndMatchedSegList[i];}
  else return 0;
}


int TrackSegmentCamAtNu::GetPlaneView() const
{
  return fPlaneView;
}


double TrackSegmentCamAtNu::GetScore(vector<TrackSegmentCamAtNu*> *BegSegBank, vector<TrackSegmentCamAtNu*> *EndSegBank)
{
  int begplane2, begplane1, endplane1, endplane2;
  int nplane, plane;
  
  vector<ClusterCamAtNu*> TempContainer;
  bool IsInTemp;

  
  //Store the beginning and end of the segment locally
  begplane1=this->GetBegPlane();
  endplane1=this->GetEndPlane();
  begplane2=this->GetBegPlane();
  endplane2=this->GetEndPlane();
  //
  //  begplane2============begplane1===========endplane1============endplane2
  //            BegSegBank             this              EndSegBank
  //

  // Store pointers to all clusters
  unsigned int nclusters = this->GetEntries();  
  for(unsigned int i=0; i<nclusters; ++i) {
    TempContainer.push_back(this->GetCluster(i));
  }

  unsigned int nsegments;
  unsigned int nclustersintemp;
  int tempplane;
  if(BegSegBank) {
    nsegments = (*BegSegBank).size();
    for(unsigned int j=0; j<nsegments; ++j) {
      //looking for the earliest plane in the Beginning Segments
      tempplane = (*BegSegBank)[j]->GetBegPlane();
      if(tempplane<begplane2) begplane2=tempplane;
      //Loop over clusters in this segment and add any we don't have
      nclusters=(*BegSegBank)[j]->GetEntries();
      for(unsigned int k=0; k<nclusters; ++k) {
        ClusterCamAtNu* clust = (*BegSegBank)[j]->GetCluster(k);
        
        IsInTemp=false;
        nclustersintemp = TempContainer.size();
        for(unsigned int l=0; l<nclustersintemp; ++l) {
          if(TempContainer[l]==clust) {IsInTemp=true; break;}
        }
        if(IsInTemp==false) {TempContainer.push_back(clust);}
      }
    }
  }

  if(EndSegBank) {
    nsegments = (*EndSegBank).size();
    for(unsigned int j=0; j<nsegments; ++j) {
      //looking for the latest plane in the End Segments
      tempplane = (*EndSegBank)[j]->GetEndPlane();
      if(tempplane>endplane2) endplane2=tempplane;
      //Loop over clusters in this segment and add any we don't have
      nclusters = (*EndSegBank)[j]->GetEntries();
      for(unsigned int k=0; k<nclusters; ++k) {
        ClusterCamAtNu* clust = (*EndSegBank)[j]->GetCluster(k);
        
        IsInTemp=false;
        nclustersintemp = TempContainer.size();
        for(unsigned int l=0; l<nclustersintemp; ++l) {
          if(TempContainer[l]==clust) {IsInTemp=true; break;}
        }
        if(IsInTemp==false) {TempContainer.push_back(clust);}
      }
    }
  }

  //Convert planes to a single view co-ordinate system where begplane2 = 0;
  begplane1-=begplane2; begplane1/=2;
  endplane1-=begplane2; endplane1/=2;
  endplane2-=begplane2; endplane2/=2;

  //calculate the number of planes we will be considering
  nplane=1+endplane2;
  double* T = new double[nplane];
  double* Z = new double[nplane];
  double* W = new double[nplane];
  for(int i=0; i<nplane; ++i) {T[i]=0.; Z[i]=0.; W[i]=0.;}

  int km, kp;  
  double m,c;
  double dt2,sn;
  double score, dstraightness, straightness, expected;
  double sw, swz, swt, swzt, swzz;
  unsigned int nhits;

  nclusters = TempContainer.size();

  for(unsigned int k=0; k<nclusters; ++k) {
    ClusterCamAtNu* clust = TempContainer[k];
    //Convert clust->GetPlane() to a single view co-ordinate system where begplane2 = 0;
    plane=(clust->GetPlane()-begplane2)/2;

    if(!(plane<0 || plane>=nplane)) {
      sw=0.; swz=0.; swt=0.;
      nhits = clust->GetEntries();
      for(unsigned int k1=0; k1<nhits; ++k1) {
        HitCamAtNu* hit = clust->GetHit(k1);
        
        swz+=hit->GetCharge()*hit->GetZPos();
        swt+=hit->GetCharge()*hit->GetTPos();
        sw+=hit->GetCharge();
      }
      
      if(sw>0.){
        Z[plane]=swz/sw; T[plane]=swt/sw;
        // Weight segments on planes spanned by seed segment. 
        // Deweight segments on other planes.
        if( plane+1>begplane1 && plane-1<endplane1 ) {W[plane]=5.;} 
        else {W[plane]=0.5;}
      }
    }
  }
    
    
  score=0.; straightness=0.; expected=0.;
  
  
  for(int k=0; k<nplane; ++k) {

    if(W[k]>0.){
      swz=0.; swt=0.; swzz=0.; swzt=0.; sw=0.; sn=0.;
      dstraightness=0.;

      km=k-5; kp=k+5;
            
      if(km<0) {km=0;}
      if(kp>nplane-1) {kp=nplane-1;}
      

      // Fit this section
      for(int k1=km; k1<kp+1; ++k1){
        if(W[k1]>0.) {
          swz+=W[k1]*Z[k1]; swt+=W[k1]*T[k1]; 
          swzz+=W[k1]*Z[k1]*Z[k1]; swzt+=W[k1]*Z[k1]*T[k1]; 
          sw+=W[k1]; sn+=1.;
        }
      }
      
      // Calculate deviation from fit at this plane
      if(sn>2.){
        m=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
        c=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);

        dt2=pow(T[k]-(m*Z[k]+c),2);

        if(dt2<1.e-5) {straightness+=1;}
        else {straightness+=dt2/1.e-5;}
        expected+=1;
      }
    }
  }

  // Protect against divide by zero
  if(expected==0) {expected=1; straightness=1;}

  score = 1.e4 + double(TempContainer.size()) - pow(straightness/expected,0.5);

  
  MSG("AlgTrackCamAtNuList", Msg::kVerbose) << " nclusters " << double(TempContainer.size()) 
                                        << " straightness " << pow(straightness/expected,0.5) 
                                        << " Score " << score 
                                        << " begplane " << this->GetBegPlane() << " endplane " << this->GetEndPlane()
                                        << " begtpos " << this->GetBegTPos() << " endtpos " << this->GetEndTPos() << endl;

//   MSG("AlgTrackCamAtNuList", Msg::kVerbose) << "Clusters involved: " << endl;
  
//   for(unsigned int i=0; i<TempContainer.size(); ++i) {
//     MSG("AlgTrackCamAtNuList", Msg::kVerbose) << "plane " << TempContainer[i]->GetPlane() << " begtpos " << TempContainer[i]->GetBegTPos() 
//                                        << " endtpos " << TempContainer[i]->GetEndTPos() << endl;
//   }

//   MSG("AlgTrackCamAtNuList", Msg::kVerbose) << "----------------------" << endl;


  if(score<0.) {score=0.;}
 
  if(T) {delete [] T;}
  if(Z) {delete [] Z;}
  if(W) {delete [] W;}

  return score;
}

---