CDTracker.cxx

Go to the documentation of this file.

// Program name: CDTracker.cxx
//
// Package: CalDetTracker
//
// Purpose: 
//
// Contact: Chris Smith, Ryan Nichol, Leo Jenner or Jeff Hartnell

#include "TGraphErrors.h"
#include "TF1.h"
#include "TMath.h"

#include "CalDetSI/Helpers.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexSEIdAltL.h"
#include "Plex/PlexStripEndId.h"
#include "Validity/VldContext.h"
#include "Validity/VldTimeStamp.h"

#include "CalDetTracker/CDTracker.h"
#include "CalDetTracker/CDTrackerOptions.h"
#include "CalDetTracker/CDTrackInfo.h"

using std::cout;
using std::endl;
using std::map;
using namespace CalDetConstants;

CVSID("$Id: CDTracker.cxx,v 1.16 2007/01/15 19:52:00 rhatcher Exp $");

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

CDTracker::CDTracker()
{
  MSG("CDTracker",Msg::kVerbose)
    <<"Running default CDTracker constructor..."<<endl;

  fevenresult=0;
  foddresult=0;

  MSG("CDTracker",Msg::kVerbose)
    <<"Finished default CDTracker constructor"<<endl;
}

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

CDTracker::CDTracker(map<Int_t,CandStripHandle> stripmap,
                     CDTrackerOptions *cdto)
{
  fevenresult=-1; //horizontal strips
  foddresult=-1;  //vertical strips
  fnum_cc_hits=-1;
  fStripMap = stripmap;
  
  if(!this->SetTrackOptions(cdto)){
    MSG("CDTracker",Msg::kWarning)<<"No Options set"<<endl;
  }
}

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

Int_t CDTracker::SetMap(map<Int_t,CandStripHandle> stripmap)
{
  fStripMap=stripmap;
  return 1;
}

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

Int_t CDTracker::SetTrackOptions(CDTrackerOptions *cdto)
{
  if(cdto!=0){    
    fcd_to=cdto;
    return 1;
  }
  else return 0;
}

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

map<Int_t,CandStripHandle> CDTracker::GetTrackedStripMap
(PlaneView::EPlaneView view)
{
  if(view==PlaneView::kX) {
    if(fevenresult==-1) fevenresult = this->FindStripTrack(0);    
    return fEvenPlStrips;
  }
  else if(view==PlaneView::kY) {
    if(foddresult==-1) foddresult = this->FindStripTrack(1);    
    return fOddPlStrips;
  }
  else {
    MSG("CDTracker",Msg::kWarning)
      <<"Neither Even nor Odd plane " 
      <<"tracking requested. Not Tracking. Returning entire Strip Map."
      <<endl;
  }
  return fStripMap;
}

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

map<Int_t,CandStripHandle> CDTracker::GetCCStripMap()
{
  if(fnum_cc_hits==-1) fnum_cc_hits = this->LocateCCStrips();
  return fCCStrips;
}

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

CDTrackInfo* CDTracker::GetTrackInfo()
{
  CDTrackInfo *cdti = 0;
  if(fStripMap.size()>0) {
    while(cdti==0) cdti = GetTrackStripInfo();
  }
  else{
    MSG("CDTracker",Msg::kDebug)
      <<"Stripmap is size 0; default CDTrackInfo constructed"<<endl;
    cdti=new CDTrackInfo();
  }
  return cdti;
}

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

CDTrackInfo* CDTracker::GetTrackStripInfo()
{
  //this method determines information such as:
  //total pe
  //pe per plane in track
  //number of hits in track  
  //number of odd and even hits

  MSG("CDTracker",Msg::kDebug)
      <<"Running GetTrackStripInfo method..."<<endl;

  if(fevenresult==-1) fevenresult = this->FindStripTrack(0);
  if(foddresult==-1) foddresult = this->FindStripTrack(1);
  if(fnum_cc_hits==-1) fnum_cc_hits = this->LocateCCStrips();

  Float_t stripPos[60] = {};

  Float_t weights[60] = {};
  Float_t weightvar[60] = {};
  Float_t sum[60] = {};
  Int_t num[60] = {};

  Int_t striphitlog[60][24] = {};
  Int_t nhits_even = 0;
  Int_t nhits_odd = 0;

  Float_t tot_pe = 0;
  Float_t trk_pe = 0;
  
  map<Int_t,CandStripHandle>::iterator stripmapIter = fStripMap.begin();
  map<Int_t,CandStripHandle>::iterator stripmapend = fStripMap.end();
  
  //sum up the total pe
  while(stripmapIter!=stripmapend){
    
    Int_t pln = stripmapIter->second.GetPlane();
    Int_t strp = stripmapIter->second.GetStrip();
    
    if(pln>=0 && pln<=59 && strp<24 && strp>=0)  
      tot_pe+=stripmapIter->second.GetCharge(CalDigitType::kPE);
    
    stripmapIter++;
  }
  
  MSG("CDTracker",Msg::kDebug)<<"Total charge in pe = "<<tot_pe<<endl;
  
  //get iterators to even hits
  map<Int_t,CandStripHandle>::iterator starteven = fEvenPlStrips.begin();
  map<Int_t,CandStripHandle>::iterator endeven = fEvenPlStrips.end();

  //While filling the arrays, also check for gaps between "tracked" hit strips
  //on the same plane. If a gap is found, take the biggest hit to be the real
  //track hit.
  float hits[60][24] = {};
  int keys[60][24] = {};
  int lastHitPlane = 0;

  while(starteven!=endeven){

    Int_t pln=starteven->second.GetPlane();
    Int_t strp=starteven->second.GetStrip();

    if(pln>=0 && pln<=59 && strp<24 && strp>=0) {

      //Note: not using PEs anymore to calculate angle. Instead use SigCorr
      Float_t temp_pe=starteven->second.GetCharge(CalDigitType::kSigCorr);
      
      hits[pln][strp] = temp_pe;
      keys[pln][strp] = starteven->first;

      if(pln>lastHitPlane) lastHitPlane = pln;

      num[pln]++;
      stripPos[pln] += strp*temp_pe;
      weights[pln] += strp*4.1*temp_pe;
      weightvar[pln]+=strp*strp*4.1*4.1*temp_pe;
      sum[pln] += temp_pe;
      trk_pe += starteven->second.GetCharge(CalDigitType::kPE);
      nhits_even += 1;
      striphitlog[pln][strp] = 1;
    }
    starteven++;
  }

  MSG("CDTracker",Msg::kDebug)
      <<"Number even planes hit = "<<nhits_even<<endl;

  //get iterators to odd hits
  map<Int_t,CandStripHandle>::iterator startodd = fOddPlStrips.begin();
  map<Int_t,CandStripHandle>::iterator endodd = fOddPlStrips.end();

  while(startodd!=endodd){

    Int_t pln = startodd->second.GetPlane();
    Int_t strp = startodd->second.GetStrip();
    
    if(pln>=0 && pln<=59 && strp<24 && strp>=0) {
      
      Float_t temp_pe=startodd->second.GetCharge(CalDigitType::kSigCorr);
      
      hits[pln][strp] = temp_pe;
      keys[pln][strp] = startodd->first;
      
      if(pln>lastHitPlane) lastHitPlane = pln;

      num[pln]++;
      stripPos[pln] += strp*temp_pe;
      weights[pln] += strp*4.1*temp_pe;
      weightvar[pln]+=strp*strp*4.1*4.1*temp_pe;
      sum[pln] += temp_pe;
      trk_pe += startodd->second.GetCharge(CalDigitType::kPE);
      nhits_odd += 1;
      striphitlog[pln][strp] = 1;
    }
    startodd++;
  }
  
  MSG("CDTracker",Msg::kDebug)
    <<"Number odd planes hit = "<<nhits_odd<<endl;

  //Start looking for problems with the track:
  bool anyChanges = false;

  if(!fcd_to->Cosmic()){ //if beam tracker is being used
    //Check for plane gaps:
    //Look for _first_ gap in track >3 planes long
    //can remove all hits downstream of first track segment 
    //since we only care about beam related events

    int gapCount = -1;
    int removeDownstreamOfThisPlane = 0;
    for(int i=0;i<=lastHitPlane;i++){
      MSG("CDTracker",Msg::kVerbose)
        <<"gapCount="<<gapCount<<", plane="<<i<<", sum="<<sum[i]
        <<endl;
      if(sum[i]>0) {                   //on first plane hit
        gapCount = 0;                  //set gapCount to zero
        removeDownstreamOfThisPlane = i;
      }
      else if(gapCount!=-1) {
        gapCount+=1; //only increment gapCount 
                     //after first hit plane
      }
      if(gapCount>3) { //any gap larger than 3 planes can't be good
        anyChanges=true; //throw out all hits below this
        break;
      }
    }
    
    if(anyChanges) {
      MSG("CDTracker",Msg::kDebug)
        <<"Found large gap after plane "<<removeDownstreamOfThisPlane
        <<", lastHitPlane="<<lastHitPlane<<endl;
      
      for(int i=removeDownstreamOfThisPlane+1;i<=lastHitPlane;i++){
        for(int j=0;j<24;j++){
          if(hits[i][j]!=0){
            if(i%2==0) fEvenPlStrips.erase(keys[i][j]);
            else fOddPlStrips.erase(keys[i][j]);
          }
        }
      }  
      CDTrackInfo *cdti = 0;
      return cdti;  
    }
  }

  //Check for strip gaps: 
  //only possible to have one kind of gap topology with the trackers
  //hit-gap-hit, so look for this pattern:
  for(int i=0;i<60;i++){
    int hitstp = -1;
    float stppe = 0;
    for(int j=0;j<24;j++){
      if(hits[i][j]>0) {
        if(hitstp==-1) { //found first strip in plane
          hitstp=j;
          stppe=hits[i][j];
        }
        else if (j!=hitstp+1){ //second strip hit in plane 
          //if it's not next to hitstp:

          //first check to see if surrounding planes are both hit
          //if so, pick hit most consistent
          if(i-2>=0&&i+2<PLANECONST&&num[i-2]>0&&num[i+2]>0){

            float xtra = 0;
            if(int(stripPos[i-2]/sum[i-2])<int(stripPos[i+2]/sum[i+2])) 
              xtra=0.5;
            else if(int(stripPos[i-2]/sum[i-2])>int(stripPos[i+2]/sum[i+2])) 
              xtra=-0.5;
            else xtra=0.;
            
            int predPos = int((stripPos[i-2]/sum[i-2]
                               +stripPos[i+2]/sum[i+2])/2.+xtra);

            if(predPos-hitstp==0) {
              if(i%2==0) fEvenPlStrips.erase(keys[i][j]);
              else fOddPlStrips.erase(keys[i][j]);
            }
            else if(predPos-j==0) {
              if(i%2==0) fEvenPlStrips.erase(keys[i][hitstp]);
              else fOddPlStrips.erase(keys[i][hitstp]);
            }
            else {
            //if neither seem to be where expected, resort to throwing out
            //the smaller hit from the track
              if(stppe>hits[i][j]){
                if(i%2==0) fEvenPlStrips.erase(keys[i][j]);
                else fOddPlStrips.erase(keys[i][j]);
              }
              else {
                if(i%2==0) fEvenPlStrips.erase(keys[i][hitstp]);
                else fOddPlStrips.erase(keys[i][hitstp]);
              }     
            }
          }

          //else get rid of the smaller hit from the track
          else{
            if(stppe>hits[i][j]){
              if(i%2==0) fEvenPlStrips.erase(keys[i][j]);
              else fOddPlStrips.erase(keys[i][j]);
            }
            else {
              if(i%2==0) fEvenPlStrips.erase(keys[i][hitstp]);
              else fOddPlStrips.erase(keys[i][hitstp]);
            }
          }
          anyChanges=true;
        }
      }
    }
  }

  if(anyChanges) {
    CDTrackInfo *cdti = 0;
    return cdti;  
  }


  //Do quick Sanity Check.
  //The Even and Odd Planes must overlap.
  
  Int_t firstEven=-1;
  Int_t lastEven=-1;
  Int_t firstOdd=-1;
  Int_t lastOdd=-1;
  //get the first and last planes
  for(Int_t pln=0;pln<60;pln++) {
    if(num[pln]>0) {
      if(firstEven==-1) firstEven=pln;
      lastEven=pln;
    }
    if(num[pln]>0) {
      if(firstOdd==-1) firstOdd=pln;
      lastOdd=pln;
    }
  }
  
  MSG("CDTracker",Msg::kDebug)
    <<"First even plane="<<firstEven
    <<", odd="<<firstOdd<<endl
    <<"Last even plane="<<lastEven
    <<", odd="<<lastOdd
    <<endl;
  
  //work out if the tracks overlap
  if(firstEven==-1 || firstOdd==-1) {
    //if there is not a hit in one (or both) views than set result to 0
    if(firstEven==-1) fevenresult=0;
    if(firstOdd==-1) foddresult=0;
  }
  else if(firstEven<firstOdd) {
    if(lastEven<firstOdd) {
      //They don't overlap;
      fevenresult=5;
      foddresult=5;
    }
    //else they do overlap
  }
  else {
    if(lastOdd<firstEven) {
      //They don't overlap;
      fevenresult=5;
      foddresult=5;
    }
    //else they do overlap
  }
  
  Float_t vertex1[2] = {0};
  Float_t vertex2[2] = {0};
  vertex1[0] = -1;//the lowest plane hit
  vertex1[1] = 0;//the sum of strip*pe in the lowest plane hit
  vertex2[0] = -1;//the highest plane hit
  vertex2[1] = 0;//the sum of strip*pe in the highest plane hit
  
  Float_t evenZ[30] = {0};
  Float_t evenZerr[30] = {0};
  Float_t evenT[30] = {0};
  Float_t evenTerr[30] = {0};
  Int_t evencnt = 0;//number of even planes hit
  Float_t oddZ[30] = {0};
  Float_t oddZerr[30] = {0};
  Float_t oddT[30] = {0};
  Float_t oddTerr[30] = {0};
  Int_t oddcnt = 0;
  
  for(Int_t i=0;i<60;i++){
    if(sum[i]!=0) {
      
      Int_t numstrps = 0; //number of *different* strips hit per plane
      for(Int_t lp = 0;lp<24;lp++) {
        if(striphitlog[i][lp]==1) numstrps+=1;
      }
      
      stripPos[i]/=sum[i];
      weights[i]/=sum[i];
      weightvar[i]/=sum[i];
      if(numstrps==1) weightvar[i] = 4.1/sqrt(12.);
      else {
        float temp = weightvar[i]-weights[i]*weights[i];
        if(temp>0) weightvar[i] = sqrt(temp);
        else weightvar[i] = 4.1/sqrt(12.);
      }
      if(weightvar[i]<4.1/sqrt(12.)) weightvar[i] = 4.1/sqrt(12.);

      //Fill vertex arrays
      if(vertex1[0]==-1) {
        vertex1[0] = i;
        vertex1[1] = stripPos[i];
      }
      vertex2[0] = i;
      vertex2[1] = stripPos[i];

      //Fill arrays for fitting
      if(i%2==0) {
        evenZ[evencnt] = Float_t(i)*6.0;
        evenZerr[evencnt] = 0.;
        evenT[evencnt] = weights[i];
        evenTerr[evencnt] = weightvar[i];
        evencnt++;
      }
      else {
        oddZ[oddcnt] = Float_t(i)*6.0;
        oddZerr[oddcnt] = 0.;
        oddT[oddcnt] = weights[i];
        oddTerr[oddcnt] = weightvar[i];
        oddcnt++;
      }
    }
  }
  
  TF1* form=new TF1("form","pol1",0,60);
  
  Float_t evenangle[2] = {};
  Float_t dydz[2] = {};
  Float_t oddangle[2] = {};
  Float_t dxdz[2] = {};

  if(evencnt>=3){
    
    //only loop if debug is set
    if (MsgService::Instance()->IsActive("CDTracker",Msg::kDebug)){
      for (Int_t i=0;i<evencnt;i++){
        MSG("CDTracker",Msg::kDebug)
          <<"evencnt = "<<evencnt
          <<", plane = "<<evenZ[i]
          <<", strip = "<<evenT[i]<<" +/- " << evenTerr[i]<<endl;
      } 
    }
    
    TGraphErrors* evengr=new TGraphErrors(evencnt,evenZ,evenT,
                                          evenZerr,evenTerr);  
    evengr->Fit("form","QN");
    
    MSG("CDTracker",Msg::kDebug)
      <<"Even fit parameters: icept = "<<form->GetParameter(0)
      <<", grad = "<<form->GetParameter(1)<<endl;
    
    evenangle[0]=form->GetParameter(0)/4.1; //in units of "strip"
    evenangle[1]=180.0*TMath::ATan(form->GetParameter(1))/TMath::Pi();
    dydz[0] = form->GetParameter(0)/4.1; //in units of "strip"
    dydz[1] = 6.0*form->GetParameter(1)/4.1;
    delete evengr;
  }
  else evenangle[0]=evenangle[1]=0;
  
  if(oddcnt>=3){
    //only loop if debug is set
    if (MsgService::Instance()->IsActive("CDTracker",Msg::kDebug)){
      for (Int_t i=0;i<oddcnt;i++){
        MSG("CDTracker",Msg::kDebug)
          <<"oddcnt = "<<oddcnt
          <<", plane = "<<oddZ[i]
          <<", strip = "<<oddT[i]<<" +/- "<<oddTerr[i]<<endl;
      } 
    }

    TGraphErrors* oddgr=new TGraphErrors(oddcnt,oddZ,oddT,oddZerr,oddTerr);
    oddgr->Fit("form","QN");

    MSG("CDTracker",Msg::kDebug)
      <<"Odd fit parameters: icept = "<<form->GetParameter(0)
      <<", grad = "<<form->GetParameter(1)<<endl;
    
    oddangle[0]=form->GetParameter(0)/4.1;
    oddangle[1]=180.0*TMath::ATan(form->GetParameter(1))/TMath::Pi();
    dxdz[0] = form->GetParameter(0)/4.1; //in units of "strip"
    dxdz[1] = 6.0*form->GetParameter(1)/4.1;
    delete oddgr;
  }
  else oddangle[0]=oddangle[1]=0;

  delete form;

  //check now that hits are not out of the detector for cosmic tracker
  if(fcd_to->Cosmic()){
    int Ecnt = 0;
    int Ocnt = 0;
    int *badEvenKeys = new int[fEvenPlStrips.size()];  
    int *badOddKeys = new int[fOddPlStrips.size()];
    
    starteven=fEvenPlStrips.begin(); 
    while(starteven!=endeven){
      Int_t pln=starteven->second.GetPlane();
      
      if(pln>=0 && pln<=59) {
        float transpos = dxdz[1]*float(pln) + dxdz[0];
        float transplus = dxdz[1]*float(pln+(1./12.)) + dxdz[0];
        float transminus = dxdz[1]*float(pln-(1./12.)) + dxdz[0];
        if(transpos<-0.8){
          if(transplus>-0.8||transminus>-0.8) break;
          else {
            badEvenKeys[Ecnt] = starteven->first;
            Ecnt++;
          }
        }
        else if(transpos>23.8){
          if(transplus<23.8||transminus<23.8) break;
          else {
            badEvenKeys[Ecnt] = starteven->first;
            Ecnt++;
          }     
        }
      }
      starteven++;
    }
    for(int iE=0;iE<Ecnt;iE++) fEvenPlStrips.erase(badEvenKeys[iE]);
    
    startodd=fOddPlStrips.begin(); 
    while(startodd!=endodd){
      Int_t pln=startodd->second.GetPlane();
      if(pln>=0 && pln<=59) {
        float transpos = dydz[1]*float(pln) + dydz[0];
        float transplus = dydz[1]*float(pln+(1./12.)) + dydz[0];
        float transminus = dydz[1]*float(pln-(1./12.)) + dydz[0];
        
        if(transpos<-0.8){
          if(transplus>-0.8||transminus>-0.8) break;
          else {
            badOddKeys[Ocnt] = startodd->first;
            Ocnt++;
          }
        }
        else if(transpos>23.8){
          if(transplus<23.8||transminus<23.8) break;
          else {
            badOddKeys[Ocnt] = startodd->first;
            Ocnt++;
          }     
        }
      }
      startodd++;
    }
    for(int iO=0;iO<Ocnt;iO++) fOddPlStrips.erase(badOddKeys[iO]);
    
    if(Ocnt>0||Ecnt>0) {
      CDTrackInfo *cdti = 0;
      return cdti;
    }
    delete badEvenKeys;
    delete badOddKeys;  
  }
  
  double time=double(fStripMap.begin()->second.GetVldContext()
                     ->GetTimeStamp().GetSec());
  

  //Following need to be set in TrackerModule::Ana
  double abstime = 0;  //abstime from cdlh. leave as zero for now.
  Int_t trigtime = 0;    //external trigger time. leave as zero for now.
  Float_t trk_range = 0; //track range in scint from track angle(s)
  
  CDTrackInfo *cdti=new CDTrackInfo(fevenresult,foddresult,
                                    nhits_even,nhits_odd,
                                    time,tot_pe,trk_pe,
                                    vertex1,vertex2,
                                    evenangle,oddangle,trk_range,
                                    abstime,trigtime,
                                    this->CCStripsOnTrack(evenangle,
                                                        oddangle));
  MSG("CDTracker",Msg::kDebug)
      <<"Finished GetTrackStripInfo method"<<endl;
  return cdti;
}

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

Int_t CDTracker::LocateCCStrips()
{
  Int_t num_cc_hits = 0;

  map<Int_t,CandStripHandle>::iterator stripmapIter=fStripMap.begin();
  map<Int_t,CandStripHandle>::iterator stripmapend=fStripMap.end();
  
  while(stripmapIter!=stripmapend){
    Int_t plane=stripmapIter->second.GetPlane();
    if(plane>=61 && plane<=64) {
      fCCStrips[stripmapIter->first] = stripmapIter->second;      
      num_cc_hits+=1;
    }
    stripmapIter++;
  }
  return num_cc_hits;
}

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

bool CDTracker::CCStripsOnTrack(Float_t* evenangle,Float_t* oddangle)
{
  //get rid of warnings for now JJH
  Float_t* hack1;
  Float_t* hack2;
  hack1=evenangle;
  hack2=oddangle;

  if(fnum_cc_hits==-1) fnum_cc_hits = this->LocateCCStrips();

  bool cchitsontrack = false;
  
  if(fnum_cc_hits==0){
    return cchitsontrack;
  }
  else cchitsontrack=true;
  
  return cchitsontrack;
}

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

---