|
// $Id: BrLocalTrack.cxx,v 1.1.1.1 2001/06/21 14:55:03 hagel Exp $ // $Author: hagel $ // $Date: 2001/06/21 14:55:03 $ // // BrLocalTrack is a BRAHMS data class for storing information for // Localtrack i.e. those created temporarely during local tracking. // The class contains information on associated BrDetectorHit (s) and // Trackgroup (i.e. ghost candidates) // #include <BrLocalTrack.h> #include <BrDetectorHit.h> #include <BrMath.h> #include <iostream.h> #include <iomanip.h> //___________________________________________________________________ ClassImp(BrLocalTrack); //___________________________________________________________________ BrLocalTrack::BrLocalTrack() { fPos = new Float_t [3]; fPosError = new Float_t [3]; fVec = new Float_t [3]; fVecError = new Float_t [3]; for(Int_t i=0;i<3;i++) { fPos[i] = 0.0; fPosError[i] = 0.0; fVec[i] = 0.0; fVecError[i] = 0.0; } fCovar = 0; fVtrack = 0; fTr1 = 0; fTr2 = 0; fNhit = 0; fFlg = 0; fStatus = kOk; fChisq = 0.0; } //___________________________________________________________________ BrLocalTrack::BrLocalTrack(BrLocalTrack* loctra_p) { // // Copy constructor // Int_t i; SetPos(loctra_p->GetPos()); SetVec(loctra_p->GetVec()); SetNhit(loctra_p->GetNhit()); SetFlg(loctra_p->GetFlg()); SetStat(loctra_p->GetStat()); SetChisq(loctra_p->GetChisq()); SetVtrack(loctra_p->GetVtrack()); SetTr1(loctra_p->GetTr1()); SetTr2(loctra_p->GetTr2()); Int_t numtrahit = loctra_p->GetTrackEntries(); for(i=0;i<numtrahit;i++) { AddTrackHit(loctra_p->GetTrackHitAt(i)); } } //___________________________________________________________________ BrLocalTrack::~BrLocalTrack() { //Just clear these hit arrays because they are the responsibility of //someone else to delete them. fDetectorHits.Clear(); fTrahit.Clear(); delete [] fPos; delete [] fPosError; delete [] fVec; delete [] fVecError; if(fCovar) { for(Int_t i = 0;i < 4;i++) delete [] fCovar[i]; delete fCovar; fCovar = 0; } } //________________________________________________________________________ ostream& operator<< (ostream & os, BrLocalTrack *loctra_p) { os<<"Local track; Number of Hits ="<<loctra_p->GetNhit(); os<<"Flag, Stat = "<<loctra_p->GetFlg()<<","<<loctra_p->GetStat()<<endl; os<<"Track origin = ("; os<<setw(8)<<loctra_p->GetPos()[0]<<","; os<<setw(8)<<loctra_p->GetPos()[1]<<","; os<<setw(8)<<loctra_p->GetPos()[2]<<")"; os<<endl; os<<"Track vector = ("; os<<setw(8)<<loctra_p->GetVec()[0]<<","; os<<setw(8)<<loctra_p->GetVec()[1]<<","; os<<setw(8)<<loctra_p->GetVec()[2]<<")"; os<<endl; os <<"Track Hits " << endl; Int_t i = 0; for(i = 0; i < loctra_p->GetNhit(); i++) { const BrDetectorHit *const hit_p = (BrDetectorHit*) loctra_p->GetDetectorHitAt( i ); os << setw(8) << hit_p->GetPos()[0] << setw(8) << hit_p->GetPos()[1] << setw(8) << hit_p->GetPos()[2] << endl; } return os; } //________________________________________________________________________ void BrLocalTrack::Fit() { // Fit the localtrack to the hits // Fill postions and vector in track as well as chisq and covariance matrix // 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 i=0;i<4;i++) { beta[i][0] = (Double_t)0.0; for(int j=0;j<4;j++) { alpha[i][j] = (Double_t)0.0; } } //Is there a more elegant way to do this? for(int i=0;i<4;i++) { alphaa[i] = &alpha[i][0]; betaa [i] = &beta[i][0]; } //Setup matrix equation if(fNhit<2) { cerr<<"Too few points for fitting; this ought not to happen"<<endl; return; } Int_t i = 0; for(i = 0; i < fNhit; i++) { const BrDetectorHit *const hit_p = (BrDetectorHit*) fDetectorHits.At( i ); a[0] = hit_p->GetPos()[2]; a[1] = 1; a[2] = hit_p->GetPos()[2]; a[3] = 1; Float_t dx = hit_p->GetDpos()[0]; if( dx < 0.0001 ) { // This unfortunately do happen at the moment dx = 0.02; // 18May00 Warning( "Fit" , "Hit Width dx == 0" ); } Float_t dy = hit_p->GetDpos()[1]; if( dy < 0.0001 ) { // This unfortunately do happen at the moment dy = 0.02; // 18May00 Warning( "Fit" , "Hit Width dy == 0" ); } for(int i=0;i<2;i++) { for(int j=0;j<2;j++) alpha[i][j] += a[i]*a[j]/dx/dx; beta[i][0] += a[i] * hit_p->GetPos()[0]/dx/dx; } for(int i=2;i<4;i++) { for(int j=2;j<4;j++) alpha[i][j] += a[i]*a[j]/dy/dy; beta[i][0] += a[i] * hit_p->GetPos()[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); // // Update found position... // if ( ifail != 0 ) Error( "BrLocalTrack::Fit", "Inversion of matrix failed!!!!!" ); fPos[0] = beta[1][0]; fPos[1] = beta[3][0]; fPos[2] = 0.0; fVec[0] = beta[0][0]; fVec[1] = beta[2][0]; fVec[2] = 1.0; //Calculate chisqr fChisq = (Float_t)0.0; for(int ic=0; ic< fNhit; ic++) { const BrDetectorHit *const hit_p = (BrDetectorHit*)fDetectorHits.At(ic); a[0] = hit_p->GetPos()[2]; a[1] = 1; a[2] = hit_p->GetPos()[2]; a[3] = 1; Float_t dx = hit_p->GetDpos()[0]; Float_t dy = hit_p->GetDpos()[1]; Float_t dchi1 = (Float_t)(beta[1][0] + beta[0][0]*a[0] - hit_p->GetPos()[0]); Float_t dchi2 = (Float_t)(beta[3][0] + beta[2][0]*a[2] - hit_p->GetPos()[1]); fChisq += dchi1*dchi1/(dx*dx) + dchi2*dchi2/(dy*dy); } if(fNhit == 2) fChisq = 0.0; else fChisq = fChisq/(fNhit-2); if(!fCovar) { fCovar = new Float_t* [4]; for(i=0;i<4;i++) fCovar[i] = new Float_t [4]; } for(i = 0; i < 4; i++) for(Int_t j = 0; j < 4; j++) fCovar[i][j] = alpha[i][j]; // The calculation below is uncessary // The covariance matrix already contains this information and // more , namely the correlations. At this point the data members // are not removed but the calculations skipped and replace by // a simple assignment to the covar matrix elements. // // Calculate the errors on the fit according to Bevington p 108-9 // Bevington, Robinson : // "Data Reduction and Error Analysis for the Physical Sciences" // if x = az + b and ey is known for each point then // // ea = 1/delta * sum z_i^2 / ex_i^2 // eb = 1/delta * sum 1 / ex_i^2 // where // ( sum 1 / ex_i^2 sum z_i / ex_i^2 ) // delta = det ( ) // ( sum z_i / ex_i^2 sum z_i^2 / ex_i^2 ) // // The same method will be employed here. First for x (ipos = 0) // and then for y (ipos = 1) /* for(Int_t ipos = 0; ipos < 2; ipos++) { Double_t deltaMatrix[2][2] = { {0, 0}, {0, 0} }; for(int ihit=0; ihit< fNhit; ihit++) { const BrDetectorHit *const hit_p = (BrDetectorHit*)fDetectorHits.At(ihit); Double_t z = hit_p->GetPos()[2]; Double_t error_2 = TMath::Power(hit_p->GetDpos()[ipos], 2); deltaMatrix[0][0] += 1 / error_2; deltaMatrix[0][1] += z / error_2; deltaMatrix[1][0] += z / error_2; deltaMatrix[1][1] += TMath::Power(z,2) / error_2; } Double_t delta = deltaMatrix[0][0] * deltaMatrix[1][1] - deltaMatrix[0][1] * deltaMatrix[1][0]; Double_t ea = 1 / delta * deltaMatrix[1][1]; Double_t eb = 1/delta * deltaMatrix[0][0]; ea = TMath::Sqrt(1 / delta * deltaMatrix[1][1]); eb = TMath::Sqrt(1/delta * deltaMatrix[0][0]); fVecError[ipos] = ea; fPosError[ipos] = eb; } */ fPosError[0] = TMath::Sqrt(fCovar[0][0]); fVecError[0] = TMath::Sqrt(fCovar[1][1]); fPosError[1] = TMath::Sqrt(fCovar[2][2]); fVecError[1] = TMath::Sqrt(fCovar[3][3]); fVecError[2] = 0.0; fPosError[2] = 0.0; } //________________________________________________________________________ Float_t BrLocalTrack::RelativeOverlap(BrLocalTrack* t2, Float_t zl, Float_t r) { // // Calculate the relative overlap between "this" track and the given track // // Parameters: // In: dt_p pointer to the other track // zl length of detector (in local z direction) // r radius of cylinder around tracks // // Algorithm: // A cylinder of radius r is put around each track. // The relative overlap is defined as the ratio between // the volume of the overlap region of two cylinders // and // the volume of the cylinder around "this" track is calculated. // The cylinders are distorted in order to have a circular intersection // between the cylinders and planes at arbitrary z. // // Requires: // The tracks are parametrized as: // x = ax * z + bx // y = ay * z + by // const Float_t twor = 2*r; const Float_t halfl = zl/2; // Difference between the two tracks const Float_t dax = GetVec()[0] - t2->GetVec()[0]; const Float_t day = GetVec()[1] - t2->GetVec()[1]; const Float_t dbx = GetPos()[0] - t2->GetPos()[0]; const Float_t dby = GetPos()[1] - t2->GetPos()[1]; Float_t overlap = 0; const Float_t a = dax*dax + day*day; if (a == 0) { // Tracks are parallell const Float_t d = BrMath::Sqrt(dbx*dbx+dby*dby)/twor; if (d < 1) overlap = BrMath::ACos(d)-d*BrMath::Sqrt(1-d*d); } else { const Float_t b = (dax*dbx+day*dby)/a; const Float_t c = (dbx*dbx+dby*dby)/a; Float_t d = b*b-c+twor*twor/a; if (d > 0) { // Cylinders intersect d = BrMath::Sqrt(d); Float_t zn = -b - d; Float_t zx = -b + d; if (zn < -halfl) zn = -halfl; if (zx > halfl) zx = halfl; if (zn < zx) { // Cylinders intersect inside the detector // // Integrate // acos(d(z)/(2r))-d(z)/(2r)*sqrt(1-(d(z)/(2r))**2) // where // d(z)/(2r)=sqrt(a*((z+2*b)*z+c))/(2*r) // from zn to zx // // Quick and dirty implementation: Simple N-point rule const Int_t n = 5; const Float_t dz = (zx-zn)/n; for (Int_t i = 0; i < n; i++) { Float_t z = zn + dz*(i+0.5); d = BrMath::Sqrt(a*((z+2*b)*z+c))/twor; overlap += BrMath::ACos(d)-d*BrMath::Sqrt(1-d*d); } overlap *= dz/zl; } } } return overlap/(BrMath::Pi/2); } //________________________________________________________________________ BrVector3D BrLocalTrack::GetDeviationFromFit(const BrDetectorHit *const detHit) const { // // Returns a vector with the deviation in x and y of the hit from // the fitted line. // The deviation is calculated for the same z // and is given as fit - hit. if( fPos == 0 || fVec == 0 ) { Warning("GetDeviatioFromFit", "No fit data"); return BrVector3D(); } BrVector3D origin( fPos ); BrVector3D slope( fVec ); BrVector3D hit( detHit->GetPos() ); return BrVector3D((hit[2] - origin[2])*slope + origin - hit); } //________________________________________________________________________ Bool_t BrLocalTrack::HasHitInRow(Int_t rowNo) const { // // Look through the hits to see if there is hit in the given // row. returns kTRUE if it finds a hit for( Int_t i = 0; i < fNhit; i++ ) { BrDetectorHit* hit_p = (BrDetectorHit*) fDetectorHits.At( i ); if(hit_p->GetRow() == rowNo) return kTRUE; } return kFALSE; } //________________________________________________________________________ Bool_t BrLocalTrack::HasAllHitsBetween(Int_t index,Float_t low, Float_t high) const { // // Look through the hits to see if the index component of all // position coordinates is between low and high // returns kTRUE if it does not find anyone outside // Useful if you only want tracks in a certain region of the TPC, no // matter if it is pad, time or row. if( index < 0 || index > 2 ) { Warning("HasAllHitsBetween", "Not a valid index!"); return kFALSE; } for( Int_t i = 0; i < fNhit; i++ ) { BrDetectorHit* hit_p = (BrDetectorHit*) fDetectorHits.At( i ); if(hit_p->GetPos()[index] < low || hit_p->GetPos()[index] > high) return kFALSE; } return kTRUE; } //________________________________________________________________________ void BrLocalTrack::RemoveDetectorHit(BrDetectorHit *detHit) { if(!fDetectorHits.Remove(detHit)) return; fDetectorHits.Compress(); fNhit--; } //________________________________________________________________________ // $Log: BrLocalTrack.cxx,v $ // Revision 1.1.1.1 2001/06/21 14:55:03 hagel // Initial revision of brat2 // // Revision 1.16 2001/04/14 22:18:58 videbaek // Fix couple of problems: // chisqs might be undefined since fitting could be called with fNhits=2. // The fCovar was not initialized. // Bypassed not necessary calculation of errors - used covar instead. // Added hits to output list of in os << method // // Revision 1.15 2001/02/09 18:13:43 pchristi // Added remove detector hit. // // Revision 1.14 2001/02/01 17:25:46 pchristi // Added calculations of errors to the local track module and changed some // of the member structure. // // Revision 1.13 2001/01/30 10:31:52 pchristi // Added some utility methods to BrLocalTrack. // // Revision 1.11 2000/09/04 14:56:26 videbaek // Add some redundant brackets to make visual C++ ver 6.0 happy. // // Revision 1.10 2000/05/24 07:23:32 pchristi // Method Fit : If hit width = 0, width = 0.02 and Warning. // // Revision 1.9 2000/05/05 09:50:30 bjornhs // Added RelativeOverlap in BrLocalTrack // // Revision 1.8 2000/05/01 14:19:10 videbaek // Fix serious error. fPos and fVec not set. // Thanks to T.Tveter. // // Revision 1.7 2000/04/28 21:18:58 videbaek // Updates to BrLocalTrack. Added fit method; Status is changed. Uses // MarkAsBad, IsBad instead of fixed 999. // Cleanup of BrModuleMatchtrack. Added histohgrams. // Changed LocalTrack and LocalTracking to use IsBad etc. // // Revision 1.6 1999/04/28 15:31:40 hagel // Change destructor to clear fDetectorHits and fTraHit // // Revision 1.5 1999/03/07 22:38:03 videbaek // Added DetectorHit and BrLocalTrackingModule as a more general way to // deal with local tracks and hits. Initial insertion. Not fully debugged. // Cosmetic changes to other modules. // // Revision 1.4 1999/01/15 16:37:27 videbaek // Working version of MRS tracking of mtp1,mtp2 and new general track // classes. Changes also added to BrModuleMatchTrack for this reason. // Updated makeNt for non-cygnus win95 systems // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|