BRAT 2.4.5
Class index
Full class index
brahmlib
BRAHMS
ROOT page 
/////////////////////////////////////////////////////////////////////////////////
//  
// BrTPMTrackVertexModule
//
// Vertexing module for the MRS. 
// June 1999
// J.H. Lee
// Requires: at least two local tracks in MTPC1 (or Spectrometer Tracks)
//

//
// $Id: BrTPMTrackVertexModule.cxx,v 1.4 2002/01/03 19:54:02 cholm Exp $
// $Author: cholm $
// $Date: 2002/01/03 19:54:02 $
// $Copyright: 2000 Brahms Collaboration
// 
///////////////////////////////////////////////////////////////////////


#include "BrTPMTrackVertexModule.h"
#ifndef BRAT_BrTableNames
#include "BrTableNames.h"
#endif
#ifndef BRAT_BrTpcTrackCandidate
#include "BrTpcTrackCandidate.h"
#endif
#ifndef ROOT_TDirectory
#include "TDirectory.h"
#endif

ClassImp(BrTPMTrackVertexModule);


 BrTPMTrackVertexModule::BrTPMTrackVertexModule():BrModule()
{
  // Default Constructor.
  SetState(kSetup);  
  fFrontVolume  = 0;
  fVertex_useplane = 1;
}

//__________________________________________________________________________
 BrTPMTrackVertexModule::BrTPMTrackVertexModule(Char_t *Name,Char_t *Title):BrModule(Name,Title)
{  
  SetState(kSetup);  
  fFrontVolume  = 0;
  fVertex_useplane = 1;
  SetDefaultParameters();
}

//___________________________________________________________________
 void BrTPMTrackVertexModule::DefineHistograms(){
  //
  // define diagnostics histograms for track combination module 
  // Book histograms (called from BrModule::Book)

  TDirectory* saveDir = gDirectory; 
  TDirectory* histDir = gDirectory->mkdir("TPM1_TrackVertex");
  histDir->cd();
  
  hVertex_x = new TH1F("hVertex_x", "Vertex position in x", 
		       200, -5., 5.);
  hVertex_y = new TH1F("hVertex_y", "Vertex position in y", 
		       200, -5., 5.);
  hVertex_z = new TH1F("hVertex_z", "Vertex position in z", 
		       200, -5., 5.);
  hVertex_z_wide = new TH1F("hVertex_z_wide", "Vertex position in z", 
			    160, -40., 40.);
  hVertex_ca = new TH1F("hVertex_ca", "Closest Approach", 
			200, 0., 20.);
  hVertex_chisq = new TH1F("hVertex_chisq", "Vertex chisq", 
			   100, 0., 20.);
  hVertex_y_vs_x = new TH2F("hVertex_y_vs_x", "Vertex position y vs x", 
			    80,-5.,5.,80,-5.,5.);
  hVertex_x_vs_z = new TH2F("hVertex_x_vs_z", "Vertex position x vs z", 
			    140,-35.,35.,100,-10.,10.);
  hVertex_x_vs_z_prof = new TProfile("hVertex_x_vs_z_prof", 
				     "Profile Vertex x vs z", 
				     70,-35.,35.,-10.,10.);
  
  fHistOn = kTRUE;


  gDirectory = saveDir;
}

//___________________________________________________________________
 void BrTPMTrackVertexModule::SetDefaultParameters()
{

  // Set/Reset internal vertexing parameters
  // to default values.
  // Initially, put some "loose"(?) cuts
  
  fMinNoTracks = 3;
  fVertex_dist_cut   = 2.; // tracks 10cm away form  
  fVertex_chisq_cut  = 2.; // chisq per d.o.f. cut
 
}

//___________________________________________________________________
 void BrTPMTrackVertexModule::Init()
{
  SetState(kInit);
  // Initialize parameters
  BrGeometryDbManager* gGeomDb = BrGeometryDbManager::Instance();
  fFrontVolume  = (BrDetectorVolume*)
    gGeomDb->GetDetectorVolume("BrDetectorVolume", GetName());

  if (fVertex_useplane) {

    Double_t TPM_theta_deg = fFrontVolume->GetRotMatrix()->GetTheta();
    if (DebugLevel()) cout << " Angle Theta is: " << TPM_theta_deg << endl;
    Double_t TPM_theta_rad = TPM_theta_deg*3.14159/180;

    fVert_plane = new BrPlane3D(sin(TPM_theta_rad),0,cos(TPM_theta_rad),0);
 
  }
}

//___________________________________________________________________________________
 BrTPMTrackVertexModule::~BrTPMTrackVertexModule()
{
  //Destructor
}

//____________________________________________________________________________________
 void BrTPMTrackVertexModule::Event(BrEventNode* InputNode,BrEventNode* OutputNode)
{

  SetState(kEvent);
  Float_t   Vertex[3];
  for (Int_t i = 0;i<=2;i++) Vertex[i] = 999;

  Float_t Vertex_chisq = 999;

  fVertex_found = 0;

  if(DebugLevel()>2){
    cout<< "Vertexing routine***" << endl;
  }if(!fFrontVolume){
    Error("Event","Module not properly setup");
    return;
  }
  if(DebugLevel()>2){
    cout<< "Input table name Vertexing routine:" << endl;
    InputNode->ListObjects(); 
  }

  BrDataTable* FrontDetectorTracks = InputNode->GetDataTable(Form("%s %s", BRTABLENAMES kTpcTrack, GetName()));
  if (!FrontDetectorTracks) {
     FrontDetectorTracks = OutputNode->GetDataTable(Form("%s %s", BRTABLENAMES kTpcTrack, GetName()));
    if (!FrontDetectorTracks) {
      Warning("Event","TPM1 Tracktable not found");
      return;
    }

  } 

  BrVertex* vtxdata = new BrVertex("TPM1 Track VtxData", "TPM1 Track Vertex Data Object");
 
  if (FrontDetectorTracks) {
    FindVertex(FrontDetectorTracks,vtxdata);
  }
  else {
    if(DebugLevel()) cout << "Detector Tracks not set - vertexing not done!" << endl;
    vtxdata = 0;
  }

  OutputNode->AddObject(vtxdata); 

}

//___________________________________________________________________________________
 void BrTPMTrackVertexModule::FindVertex(BrDataTable *FrontDetectorTracks, BrVertex* vtxdata)
  // Algorithm based on the CERN vertexing routines I used for the
  // E866 - JH Lee
  // Input: local tracks in TPC1, but can be spectrometer tracks
  //
  // Get Vertex with every tracks 
  //  select tracks with chisq and/or distance from the vertex
{

// Get vertex position, closest approach distances

  float pdotn,ssq,hicda,dist;
  float t11,t12,t13,t22,t23,t33,det,nvhi,sq,wsum,previous_chisq;
  int i,k,n,l,nloop;
  const Int_t max_track_for_vtx=100;
  
  //don't deal with events which have more than 100 tracks  
  //const int max_track=100;
  
  float vec[5][max_track_for_vtx],cda[4][max_track_for_vtx],ww[max_track_for_vtx],vert[3];
  float dd[3],am[3][3],cosv[3];

  Float_t   Vertex[3];
  for (Int_t i = 0;i<=2;i++) Vertex[i] = 999;

  Float_t Vertex_chisq = 999;
  Int_t Vertex_number_of_tracks = 0;

  Float_t   Vertex_ca[100];                //  Closest distance of approac
  //
  //  ***LOOP ON TRACKS
  //Let's first start with local tracks
  //It can also be from combined tracks
  
  //first check if there are tracks in MTPC1
  
  BrVector3D xg, ag, fIntersect;
  BrTpcTrackCandidate* FrontTrack_p;
  
  int NumFrontTr;
  int number_of_out_of_range;
  int ntrack_used;
  
  if(!FrontDetectorTracks){
    if(DebugLevel()>2){
      cout<< "Can't vertex tracks because Detector Tracks not present ";
      cout<< fFrontVolume->GetName() << endl;
    }
    NumFrontTr = 0;
  }
  else
    NumFrontTr = FrontDetectorTracks->GetEntries();
  
  //Not enough tracks to do vertexing, there should be at least two tracks for
  //vertexing
  
   if(DebugLevel()>2){
     cout<< "Number of Front Tracks From FindVertex:" << NumFrontTr << endl;
   }

   if(NumFrontTr < fMinNoTracks) return;
   
   for(i=0;i<NumFrontTr;i++) {
     FrontTrack_p = (BrTpcTrackCandidate*) FrontDetectorTracks->At(i);
   
     //
     // The local track coordinates are given in the local coordinates
     // for TPM1.  Convert to Global 
     //   
     
     if(DebugLevel()>2){
       cout<< "volume name:" << fFrontVolume->GetName() << endl;
       cout<< "Vertex0: position" << FrontTrack_p->GetPos() << endl;
       cout<< "Vertex0: alpha" << FrontTrack_p->GetSlope() << endl;
     }

     fFrontVolume->LocalToGlobal(FrontTrack_p->GetPos(),xg,0);
     //     fFrontVolume->LocalToGlobal(FrontTrack_p->GetAlpha(),ag,1);
     fFrontVolume->LocalToGlobal(FrontTrack_p->GetSlope(),ag,1);
     
     if(DebugLevel()>2){
       cout<< "Vertex0: xg[0]" << xg[0] << endl;
       cout<< "Vertex0: xg[1]" << xg[1] << endl;
       cout<< "Vertex0: xg[2]" << xg[2] << endl;
       cout<< "Vertex0: ag[0]" << ag[0] << endl;
       cout<< "Vertex0: ag[1]" << ag[1] << endl;
       cout<< "Vertex0: ag[2]" << ag[2] << endl;
     }
     
     if (fVertex_useplane) { 
       // Project to plane paralel to front of MTPC1
       //     BrLine3D *fTrack = new BrLine3D(xg,ag);
       //     BrVector3D* fIntersect;
       //     fIntersect = fVert_plane->GetIntersectionWithLine(fTrack);
       // Project to plane - taken from GetIntersectionWithLine in BrLine3D
       Double_t u; 
       
       Double_t fA = fVert_plane->GetA();
       Double_t fB = fVert_plane->GetB();
       Double_t fC = fVert_plane->GetC();
       Double_t fD = fVert_plane->GetD();
       
       u = - ( fA * xg[0]   + fB * xg[1]   + fC * xg[2]   + fD ) / 
	 ( fA * ag[0] + fB * ag[1] + fC * ag[2]      );
       
       fIntersect[0] = xg[0] + u*ag[0];
       fIntersect[1] = xg[1] + u*ag[1];
       fIntersect[2] = xg[2] + u*ag[2];
       
       vec[0][i]=fIntersect[0];
       vec[1][i]=fIntersect[1];
       vec[2][i]=fIntersect[2];
       vec[3][i]=ag[0]/ag[2];
       vec[4][i]=ag[1]/ag[2];

     }
     else {       // Project to vertex (z=0)
       
       float t = -xg[2] / ag[2];
       float xt = xg[0] + t * ag[0];
       float yt = xg[1] + t * ag[1];
       
       vec[0][i]=xt;
       vec[1][i]=yt;
       vec[2][i]=0.;
       vec[3][i]=ag[0]/ag[2];
       vec[4][i]=ag[1]/ag[2];
       
     }
     
     if(DebugLevel()>2){
       cout << "Vertex:projection:x " << vec[0][i] << endl;
       cout << "Vertex:projection:y " << vec[1][i] << endl;
       cout << "Vertex:projection:z " << vec[2][i] << endl;
     }
     
     Vertex_ca[i] = 0.;
     
     // Initialize all weights to 1
     ww[i]=1.;           
     
   }
   
   Vertex_number_of_tracks = 0;
   Vertex_chisq = 100.;
   nloop = 1;
   
   if(DebugLevel()>0){
     cout<< "Total number of tracks " << NumFrontTr << endl;
   }
   //loop until track < 4 or outot range track = 0
   do {
     
     for(k=0;k<=2;k++){
       am[k][0]=0.;
       am[k][1]=0.;
       am[k][2]=0.;
       dd[k]=0;
     }
     
     previous_chisq = Vertex_chisq ; 
     
     ntrack_used=0;  
     
     for(n=0; n < NumFrontTr; n++) {
       if(Vertex_ca[n] >= 0.){
	 if((vec[3][n]*vec[3][n]+vec[4][n]*vec[4][n]+1.)<0.0){
	   if(DebugLevel()>0) cout << "1: VERTEX MATRIX SINGULAR VERTEX AND CDA SET TO 0.";
	   ssq = 100000000.;
	   return;
	 }
	 ntrack_used=ntrack_used+1;
	 
	 // Set track weight to inverse of closest approach
	 if (nloop > 1) ww[n]=1./Vertex_ca[n]; 
	 
	 cosv[2]=sqrt(1./(vec[3][n]*vec[3][n]+vec[4][n]*vec[4][n]+1.));
	 cosv[0]=cosv[2]*vec[3][n];
	 cosv[1]=cosv[2]*vec[4][n];
	 
	 am[0][0]=am[0][0]+(cosv[0]*cosv[0]-1.)*ww[n];
	 am[1][1]=am[1][1]+(cosv[1]*cosv[1]-1.)*ww[n];
	 am[2][2]=am[2][2]+(cosv[2]*cosv[2]-1.)*ww[n];
	 am[0][1]=am[0][1]+cosv[0]*cosv[1]*ww[n];
	 am[0][2]=am[0][2]+cosv[0]*cosv[2]*ww[n];
	 am[1][2]=am[1][2]+cosv[1]*cosv[2]*ww[n];
	 
	 if(DebugLevel()>2){
	   cout<< "Vertex: vec[0][n]" << vec[0][n] << endl;
	   cout<< "Vertex: vec[1][n]" << vec[1][n] << endl;
	   cout<< "Vertex: vec[2][n]" << vec[2][n] << endl;
	   cout<< "Vertex: vec[3][n]" << vec[3][n] << endl;
	   cout<< "Vertex: vec[4][n]" << vec[4][n] << endl;
          }
	 
	 pdotn=vec[0][n]*cosv[0]+vec[1][n]*cosv[1]+vec[2][n]*cosv[2];
	 
	 dd[0]=dd[0]+(pdotn*cosv[0]-vec[0][n])*ww[n];
	 dd[1]=dd[1]+(pdotn*cosv[1]-vec[1][n])*ww[n];
	 dd[2]=dd[2]+(pdotn*cosv[2]-vec[2][n])*ww[n];
       }
     }
     
     t11=am[1][1]*am[2][2]-am[1][2]*am[1][2];
     t12=am[1][2]*am[0][2]-am[0][1]*am[2][2];
     t13=am[0][1]*am[1][2]-am[1][1]*am[0][2];
     t22=am[0][0]*am[2][2]-am[0][2]*am[0][2];
     t23=am[0][1]*am[0][2]-am[0][0]*am[1][2];
     t33=am[0][0]*am[1][1]-am[0][1]*am[0][1];
     
      det=am[0][0]*t11+am[0][1]*t12+am[0][2]*t13;
      
      if(det == 0.0){
	if(DebugLevel()>0) cout << "2: VERTEX MATRIX SINGULAR VERTEX AND CDA SET TO 0.";
	ssq = 100000000.;
	return;
      }
      
      det=1./det;
      
      vert[0]=det*(dd[0]*t11+dd[1]*t12+dd[2]*t13);
      vert[1]=det*(dd[0]*t12+dd[1]*t22+dd[2]*t23);
      vert[2]=det*(dd[0]*t13+dd[1]*t23+dd[2]*t33);
      
      //      fVertex[0]=vert[0];
      //      fVertex[1]=vert[1];
      //      fVertex[2]=vert[2];
      
      if(DebugLevel()>0){
	cout<< "Vertex: det " << det << endl; 
	cout<< "Vertex: dd0 " << dd[0] << endl; 
	cout<< "Vertex: dd1 " << dd[1] << endl; 
	cout<< "Vertex: dd2 " << dd[2] << endl; 
	cout<< "Vertex: position0 " << vert[0] << " " 
	    << vert[1] << " " << vert[2] <<  endl;
      }
      //
      //   ***   COMPUTATION OF C.D.A AND CHISQ
      //
      ssq = 0.;
      nvhi = 1;
      hicda = 0.;
      number_of_out_of_range = 0;
      Vertex_number_of_tracks = 0;

      
      for(n=0;n<NumFrontTr;n++) { 
	if((vec[3][n]*vec[3][n]+vec[4][n]*vec[4][n]+1.)<0.0){
	  if(DebugLevel()>0) cout << "3: VERTEX MATRIX SINGULAR VERTEX AND CDA SET TO 0.";
	  ssq = 100000000.;
	  return;
	}
	
	Vertex_ca[n] = -100.;

	if(ww[n] != 0.0) {
	  Vertex_number_of_tracks = Vertex_number_of_tracks + 1;   
	  cosv[2]=sqrt(1./(vec[3][n]*vec[3][n]+vec[4][n]*vec[4][n]+1.));
	  cosv[0]=cosv[2]*vec[3][n];
	  cosv[1]=cosv[2]*vec[4][n];
          
	  dist=0.;
	  
	  for(l=0;l<=2;l++) {
	    dd[l]=vert[l]-vec[l][n];
	    dist=dist+dd[l]*cosv[l];
	  }
	  
	  for(l=0;l<=2;l++) cda[l][n]=vert[l]-vec[l][n]-dist*cosv[l];
	  
	  sq = cda[0][n]*cda[0][n] + cda[1][n]*cda[1][n] + cda[2][n]*cda[2][n];
	  
	  cda[3][n] = sqrt(sq); 
	  
	  if (sq < hicda) ssq = ssq + sq*ww[n];
	  else{
	    hicda = sq;
	    nvhi = n;
	    ssq = ssq+sq*ww[n];
	  }
	  wsum = wsum + ww[n];
	  Vertex_ca[n] = cda[3][n];
	  if(fHistOn){
	    hVertex_ca->Fill(Vertex_ca[n]);
	  } 
	  //If the distance is bigger than the cut value, make it negative 
	  //make a weight 0 also
	  if (cda[3][n] > fVertex_dist_cut) {
	    cda[3][n] = cda[3][n]*(-1.);
	    Vertex_ca[n] = cda[3][n];
	    if(DebugLevel()>0)
	      cout <<"ca (out of range) "<< Vertex_ca[n] 
		   << " for " <<n << endl; 
	    ww[n]=0.;
	    number_of_out_of_range = number_of_out_of_range  + 1;
          }
	}
      }
      
      for(l=0;l<=2;l++) Vertex[l]=vert[l];
      Vertex_chisq = ssq / wsum;
      
      // Will be some looping to select tracks more likely associated
      // with the "vertex" ex: dist should be smaller than 1cm??
      // (depend on angle...) needs some naming changes.... 
      // 
      nloop = nloop + 1;
      if(DebugLevel()>0){
	cout<< "Vertex: nloop " << nloop << endl;
	cout<< "Vertex: position " << vert[0] << " " 
	    << vert[1] << " " << vert[2] <<  endl;
	cout<< "Vertex: chi-sq " << Vertex_chisq  <<  endl; 
	cout<< "number of out of range " << number_of_out_of_range << endl;
	cout<< "number of track used " << ntrack_used << endl;
      }
   } while (Vertex_number_of_tracks > 4  && 
	    number_of_out_of_range > 0 && 
	    nloop <= 10 && 
	    Vertex_chisq < previous_chisq);
   
   fVertex_found = 1;
   
   if(fHistOn){
     hVertex_x->Fill(vert[0]);
     hVertex_y->Fill(vert[1]);
     hVertex_z->Fill(vert[2]);
     hVertex_z_wide->Fill(vert[2]);
     hVertex_chisq->Fill(Vertex_chisq);
     hVertex_y_vs_x->Fill(vert[0],vert[1]);
     hVertex_x_vs_z->Fill(vert[2],vert[0]);
     hVertex_x_vs_z_prof->Fill(vert[2],vert[0],1);
   } 
   //stop looping when track <= 4 or all track are in the range  

  // Fill vertex object
   vtxdata->SetX(Vertex[0]);
   vtxdata->SetY(Vertex[1]);
   vtxdata->SetZ(Vertex[2]);
   vtxdata->SetVertexChisq(Vertex_chisq);
   vtxdata->SetVertexFound(fVertex_found);
   vtxdata->SetVertexMethodUsed(2);
   

} 

//___________________________________________________________________________________
 void BrTPMTrackVertexModule::ListVertexParameters()
{
  // List cut values used in vertexing.
  //
  printf("**Vertexing Cut Parameter Listing: n");
  cout << "Distance cut  " << fVertex_dist_cut << " cmn";
  cout << "Chi-sq Cut    " << fVertex_chisq_cut << " (per d.o.f.)n";
  cout << endl;
}

//___________________________________________________________________________________
 void BrTPMTrackVertexModule::ListEventStatistics() {
  BrModule::ListEventStatistics();
}


// $Log: BrTPMTrackVertexModule.cxx,v $
// Revision 1.4  2002/01/03 19:54:02  cholm
// Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names
//
// Revision 1.3  2001/10/07 13:35:07  ouerdane
// Put histos in a TDirectory TPM1_TrackVertex
//
// Revision 1.2  2001/09/07 18:35:15  trine
// Changed the files BrTPMTrackVertexModule.h and .cxx to make
// the module work with new TPC code / table names and with bratmain.
// Changed class version to 0, plus a little code cleanup.
// Also made the minimum required number of tracks settable.
//
// Revision 1.1.1.1  2001/06/21 14:55:14  hagel
// Initial revision of brat2
//
// Revision 1.6  2000/11/22 21:21:25  videbaek
// Changes due to better handling of Init().
//
// Revision 1.5  2000/11/17 07:54:43  bjornhs
// Added some if (DebugLevel())'s to get rid of recurring error-messages
//
// Revision 1.4  2000/11/07 16:40:24  bjornhs
// Added BrVertexModule to CVS - superclass for vertex-finding.
// Major change in TPC-vertex-modules:
// BrTPMTrackVertexModule adds BrVertex-object "TPM1 Track VtxData" to outputnode.
// BrTPMClusterVertexModule has been thorougly cleaned up to work better and faster. Adds BrVertex-object "TPM1 Cluster VtxData" to outputnode.
// BrVertex is slightly edited to reflect these changes.
//
// Revision 1.3  2000/06/27 10:03:00  bjornhs
// Added Init() to fix more initialization-problems
//
// Revision 1.2  2000/06/21 03:07:43  jhlee
// initialize  fVertex_number_of_tracks=0
//
// Revision 1.1  2000/06/10 14:43:42  bjornhs
// 1) Renamed vertexing-classes.
// Old BrVertex is now BrTPMTrackVertexModule
// Old BrHCVertex is now BrTPMClusterVertexModule
// Old BrVertexData is now BrVertex
//
// 2) Made the syntax more universal. Common usage:
//    - Generate vertexmodule-object
//    - In event-loop, use SetDetectorTracks() or SetDetectorHits()
//      to pass data to the vertexmodule
//    - vertexmodule->Event(inputnode, outputnode) adds a BrVertex-object
//      called "VtxData" to outputnode
//    - use BrVertex->GetVertex_found to determine if vertex was successfully
//      located
//
// Revision 1.5  2000/05/10 15:56:54  nbi
// Removed some multiline comments that gave annoying warning messages.
//
// Revision 1.4  2000/04/26 19:57:26  videbaek
// add Id
//
// Revision 1.3  2000/02/29 12:59:47  bjornhs
// Added histogram-based vertexfinder (BrHCVertex) and new data-class
// (BrVertexData). Modified old BrVertex to use BrVertexData
//
// Revision 1.2  1999/10/21 21:02:09  jhlee
// adding ID and Log
//

This page automatically generated by script docBrat by Christian Holm

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

Validate HTML
Validate CSS