//____________________________________________________________________
// 
// Class for Tpc track candidates. 
// In addition to the inherited information it also holds information
// on the average energy pr hit in the track 
// 

//____________________________________________________________________
//
// $Id: BrTpcTrackCandidate.cxx,v 1.7 2002/03/22 21:55:59 videbaek Exp $
// $Author: videbaek $
// $Date: 2002/03/22 21:55:59 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//
#ifndef BRAT_BrTpcTrackCandidate
#include "BrTpcTrackCandidate.h"
#endif
#ifndef BRAT_BrMath
#include "BrMath.h"
#endif
#ifndef WIN32
#include <iostream>
#include <iomanip>
#else
#include <iostream.h>
#include <iomanip.h>
#endif

//____________________________________________________________________
ClassImp(BrTpcTrackCandidate);

//____________________________________________________________________
 BrTpcTrackCandidate::BrTpcTrackCandidate()
  : BrTrackCandidate()
{
  // Default constructor - empty
  fLastRow       = 0;
  //fResiduals   = new TClonesArray("BrTrackResidual");
  fResiduals     = 0;
}

//____________________________________________________________________
 BrTpcTrackCandidate::~BrTpcTrackCandidate()
{
  if(fResiduals){
    fResiduals->Delete();
    delete fResiduals;
  }
}

//____________________________________________________________________
 BrTpcTrackCandidate::BrTpcTrackCandidate(const BrTpcTrackCandidate &track)
  : BrTrackCandidate(track)
{
  // Copy constructor
  fLastRow       = track.GetLastRow();
  fResiduals     = 0;

   if(track.GetResidualArray()) {

    const Int_t nHits = track.GetNhit();
    
    fResiduals = new TClonesArray("BrTrackResidual", nHits);
    
    for(Int_t i = 0; i < nHits; i++) {
      
      new((*fResiduals)[i]) BrTrackResidual(*(track.GetResidualAt(i)));
    }
  }
}
//________________________________________________________________________
 void BrTpcTrackCandidate::FillResiduals()
{
  // Fill clone array with residuals
  // This is also where the clones array is created 

  if(!GetHitList()) {

    Warning("FillResiduals", 
	    "No hits are associated with this track");
    return;
  }

  const Int_t nHits = BrTrackCandidate::GetNhit();
  
  if(!fResiduals)
    fResiduals = new TClonesArray("BrTrackResidual", nHits);
  else
    fResiduals->Delete();
  
  for(Int_t i = 0; i < nHits; i++) {
    
    BrTpcHit* hit = GetTpcHitAt(i);
    
    Float_t x, y;
    GetXYPositionAtZ(x, y, hit->GetPos()[2]);
    
    BrTrackResidual *res = new((*fResiduals)[i]) BrTrackResidual();
    res->SetRow(hit->GetRow());
    res->SetAdcSum(hit->GetAdcSum());
    res->SetDx(x-hit->GetPos()[0]);
    res->SetDy(y-hit->GetPos()[1]);
  }
}

//________________________________________________________________________
 void BrTpcTrackCandidate::Fit()
{
  // Fit the localtrack to the hits
  // Fill postions and vector in track as well as chisq and covariance matrix
  //
  const Int_t nHits = BrTrackCandidate::GetNhit();
  
  if(nHits < 2) {
    
    Warning("Fit", "Too few points for fitting. SHOULD NOT HAPPEN");
    return;
  } 

  Double_t a[4];
  Double_t beta[4][1];
  Double_t alpha[4][4];
  Double_t *alphaa[4],*betaa[4];
  Float_t work[4];
  
  for(Int_t i = 0; i < 4; i++) {

    beta[i][0] = (Double_t)0.0;
    alphaa[i] = &alpha[i][0];
    betaa [i] = &beta[i][0];

    for(Int_t j = 0; j < 4; j++) {

      alpha[i][j] = (Double_t)0.0;
    }
  }

  //Setup matrix equation
  for(Int_t i = 0; i < nHits; i++) {

    const BrVector3D pos      = GetTpcHitAt(i)->GetPos();
    const BrVector3D posError = GetTpcHitAt(i)->GetPosError();
    
    a[0] = pos[2];
    a[1] = 1;
    
    a[2] = pos[2];
    a[3] = 1;
    
    Float_t dx   = posError[0];
    if( dx == 0 ) {  
      dx   = 0.0001; 
      Warning( "Fit" , "Hit Width dx == 0" );
    }
    
    Float_t dy   = posError[1];
    if( dy == 0 ) {  
      dy   = 0.0001; 
      Warning( "Fit" , "Hit Width dy == 0" );
    }
    
    for(Int_t j = 0;j < 2; j++) {

      for(Int_t k = 0; k < 2; k++) 
	alpha[j][k] += a[j]*a[k]/dx/dx;

      beta[j][0] += a[j] * pos[0]/dx/dx;
    }
    
    for(Int_t j = 2; j < 4; j++) {

      for(Int_t k = 2; k < 4; k++) 
	alpha[j][k] += a[j]*a[k]/dy/dy;

      beta[j][0] += a[j] * pos[1]/dy/dy;
    }
  }

  //Solve 4x4 matrix for xc,yc and ex,ey vector
  Int_t ifail;
  BrMath::DEQINV(4,alphaa,4,work,ifail,1,betaa);

  if ( ifail != 0 ) 
    Error("Fit", "Inversion of matrix failed!!!!!" );
  
  //
  // Update found position...
  BrTrack::SetPos(beta[1][0], beta[3][0], 0);
  BrTrack::SetSlope(beta[0][0], beta[2][0], 1);

  //Calculate chisqr
  Float_t chiSquare = 0;

  for(Int_t i = 0; i < nHits; i++) {
    
    const BrVector3D pos      = GetTpcHitAt(i)->GetPos();
    const BrVector3D posError = GetTpcHitAt(i)->GetPosError();
    
    Float_t x, y;
    BrTrack::GetXYPositionAtZ(x, y, pos[2]);
    
    x -= pos[0];
    y -= pos[1];
    
    Float_t dx   = posError[0];
    if( dx == 0 ) {  
      dx   = 0.0001; 
    }
      
    Float_t dy   = posError[1];
    if( dy == 0 ) {  
      dy   = 0.0001; 
    }
      
    chiSquare += x*x/dx/dx + y*y/dy/dy;
    
  }

  // reduced chisquare (2 * (nHits -2 ) DOF)
  chiSquare /= 2*(nHits-2);
  
  BrTrack::SetQuality(chiSquare);
  
  // Set covariance matrix
  for(Int_t i = 0; i < 4; i++)
    for(Int_t j = 0; j < 4; j++)
      BrTrack::SetCovariance(i, j, alpha[i][j]);
  
  // finito
}
//____________________________________________________________________
Bool_t
 BrTpcTrackCandidate::HasHitInRow(Int_t rowNo) const
{
  // Return kTRUE if there is a hit in row rowNo
  // (kFALSE otherwise)
  
  const Int_t n      = BrTrackCandidate::GetNhit();

  for(Int_t i = 0; i < n; i++) {
    
    BrTpcHit *hit = GetTpcHitAt(i);

    if(hit->GetRow() == rowNo)
      return kTRUE;
  }

  return kFALSE;
}
//____________________________________________________________________
void
 BrTpcTrackCandidate::IncrementHitsUsed()
{
  // Increment used for all hits in track
  // (kFALSE otherwise)
  
  const Int_t n      = BrTrackCandidate::GetNhit();

  for(Int_t i = 0; i < n; i++) {
    
    BrTpcHit *hit = GetTpcHitAt(i);

    hit->IncrementUsed();
  }
}

//____________________________________________________________________
 void BrTpcTrackCandidate::Print(Option_t* option = "") const
{
  // Print out information 
  BrTrackCandidate::Print(option);
  cout << "Last hit in row : " << setw(10) << fLastRow << endl;
}

  

//____________________________________________________________________
//
// $Log: BrTpcTrackCandidate.cxx,v $
// Revision 1.7  2002/03/22 21:55:59  videbaek
// Removed the setting up of the ClonesArray in the constructor. This fails
// in the copy constructor in actual code. The class quires a match between elements
// in the clonesarray and the GetHits from the track elements.
// This shows potential for inconsistencies, as well as be awre cosmetic changes
// They may break real code too.
//
// Revision 1.6  2002/03/20 00:41:06  cholm
// Some cosmetics
//
// Revision 1.5  2001/11/13 12:27:24  pchristi
// Made residuals persistent and fixed bugs in copyconstructor.
//
// Revision 1.4  2001/11/13 10:57:32  pchristi
// Added destructor to close memory leak (shame on me) and a copy constructor.
//
// Revision 1.3  2001/11/02 13:32:51  pchristi
// Added new class BrTrackResidual. This class deals with storing track
// residuals for good tracks.
// Removed adc calculations and informations from BrTpcTrackCandidate since
// this was done in a useless way. The adcsum information is now stored with
// the residual for each hit and it makes it possible to later make a
// truncated mean if one wishes stored in the track residual
//
// Revision 1.2  2001/08/14 19:26:23  pchristi
// Removed detructors and copy constructors where the default was ok.
//
// Revision 1.1.1.1  2001/06/21 14:55:03  hagel
// Initial revision of brat2
//
// Revision 1.1  2001/06/17 17:41:37  pchristi
// The new tracking classes.
//
//

Copyright ; 2002 BRAHMS Collaboration <brahmlib@rcf.rhic.bnl.gov>
Last Update on by

Validate HTML
Validate CSS