BRAT 2.4.5
Class index
Full class index
brahmlib
BRAHMS
ROOT page
//____________________________________________________________________________
//
// BrMrsTrackingModule
//
// Tracking module for the MRS. The module will setup the module that combines
// the local tracks found front and after the D5 magnet. 
//
// Requires: BrDetectorTracks Table in inputnode of the
//           two tracking TPS TPM1 and TPM2.
//
// Note that this new version always get the fields from the database
// Please use it being aware of this.
//
// This may be a real issue when running simulated data - but someone will discover this
// one day - be forwarned. (fv)

//____________________________________________________________________________
//
// $Id: BrMrsTrackingModule.cxx,v 1.10 2002/08/30 17:36:34 videbaek Exp $
// $Author: videbaek $
// $Date: 2002/08/30 17:36:34 $
//


#include <BrIostream.h>
#include <iomanip>

#include <BrMrsTrackingModule.h>
#include <TMath.h>
#include <TDirectory.h>
#include <TString.h>
#include "BrVector3D.h"
#include "BrLine3D.h"
#include "BrPlane3D.h"
#include "BrEvent.h"
#include "BrModuleMatchTrack.h"
#include "BrDetectorTrack.h"
#include "BrGeometryDbManager.h"
#include "BrParameterDbManager.h"
#include "BrRunInfoManager.h"


//______________________________________________________________________ 
ClassImp(BrMrsTrackingModule);

//______________________________________________________________________ 
 BrMrsTrackingModule::BrMrsTrackingModule()
{
  //// default constructor
  SetState(kSetup);
  fCombineModule = 0;
  SetUseDb();
}

//______________________________________________________________________ 
 BrMrsTrackingModule::BrMrsTrackingModule(const Char_t* Name, 
					 const Char_t* Title)
  : BrModule(Name, Title)
{
  //
  // Normal constructor.
  // This Name & title is informational.
  //

  SetState(kSetup); 

  fCombineModule = new BrModuleMatchTrack("MRS","MRS matching",
					  "TPM1","TPM2","D5");
  fCombineModule->SetHistMomentumRange(3.0);
  SetUseDb();
}

//______________________________________________________________________ 
 BrMrsTrackingModule::~BrMrsTrackingModule()
{
  //// default destructor
  if(fCombineModule)
    delete fCombineModule;
}


//__________________________________________________________________________
 void BrMrsTrackingModule::Init(){
  //
  // Clear datastructures
  // request geometry information using the GeometryDbManager.
  //
  SetState(kInit);

  fCombineModule->SetDebugLevel(DebugLevel());
  fCombineModule->SetVerbose(Verbose());
  fCombineModule->Init();

  if(fCombineModule->GetStatus()== BrModule::kAbort){
    Abort("Init","Combine Module Failed");
    return;
  }
 
  fFrontVolume     = fCombineModule->GetFrontVolume();
  fBackVolume      = fCombineModule->GetBackVolume();
  fMagnetVolume    = fCombineModule->GetMagnetVolume();

  BrGeometryDbManager  *geomDb = BrGeometryDbManager::Instance();
  BrParameterDbManager *parDb  = BrParameterDbManager::Instance();

  fTofParams = (BrDetectorParamsTof*)parDb->
    GetDetectorParameters("BrDetectorParamsTof", "TOFW");
  
  fPanelPar = new BrDetectorParamsTof* [fTofParams->GetNoPanels()];
  fPanelVol = new BrDetectorVolume* [fTofParams->GetNoPanels()];
  
  for (Int_t p = 0; p < fTofParams->GetNoPanels(); p++) {
    fPanelPar[p] = 
      (BrDetectorParamsTof*)parDb->
      GetDetectorParameters("BrDetectorParamsTof", Form("TFP%d", p+1));
    fPanelVol[p] = 
      (BrDetectorVolume*)geomDb->
      GetDetectorVolume("BrDetectorVolume", Form("TFP%d", p+1));
  }
}

//__________________________________________________________________________
 void BrMrsTrackingModule::DefineHistograms()
{
  //
  // define histograms for modules owned by this 'super module'
  //
  TDirectory* savdir = gDirectory;

  TDirectory* hdir = savdir->mkdir("Tracking_MRS");
  if (!hdir) {
    Warning("DefineHistograms", "Could not create histogram subdirectory");
    return;
  }

  hdir->cd();
  fCombineModule->Book();

  hGhost = new TH2F("hGhost", "Tracks after ghostbusting vs Tracks before",
		    6, 0.5, 6.5, 6, 0.5, 6.5);
  hStatus      = new TH1F("hStatus", "Ghost busting status", 10, 0.5, 10.5);
  hTargetTz    = new TH1F("hTz","MRS target image (z)", 200, -50., 50.);
  hTargetTx    = new TH1F("hTx","MRS target image (x)", 200, -10., 10.);
  hTargetTy    = new TH1F("hTy","MRS target image (y)", 200, -10., 10.);
  hP           = new TH1F("hMomentum", "Accepted MRS tracks momentum",320, -4., 4.);
  hTheta       = new TH1F("hTheta","MRS Track Theta", 200, 0.0, 2.0);
  hPhi         = new TH1F("hPhi","MRS Track Phi", 200, -.2, .2);
  hPathLength  = new TH1F("hPathLength","MRS path length",400,300.,500.);
  hPartialPath = new TH1F("hPartialPath","MRS path length from TPM1 to TOFW",
			  400, 100., 500.);
  hSlatPos = 
    new TH1F("hSlatPos","Slat no for projected tracks", 130,-0.5, 129.5);
  
  hTrkBbDist = 
    new TH1F("hTrkDistBb", "Distance between track projection to BB "
	     "plane normal to MRS", 400, -10, 10);
  
  hTrkBbXY = 
    new TH2F("hTrkBbXY", "Track projection to BB vertex (small tubes)"
	     "plane normal to MRS (in local coord.)", 
	     200, -10, 10, 200, -10, 10);

  hTrkBbZ = 
    new TH1F("hTrkBbZ", "BB Vertex (small tubes) - Track intersection with beam line", 
	     500, -25, 25);

 // Restore directory

  savdir->cd();
}

//___________________________________________________________________
 void BrMrsTrackingModule::SetMagnetField()
{
  // ----------------------------------------------------
  // get field information from database
  // update info each time a field is needed (event loop)
  // ----------------------------------------------------

  BrRunInfoManager* runMan = BrRunInfoManager::Instance();
  
  if (!runMan) {
    Abort("SetMagnetField", "Could not get an instance of BrRunInfoManager");
    return;
  }
  
  const BrRunInfo* run = runMan->GetCurrentRun();
  if(!run){
    Abort("SetMagnetField", "Could not get current run");
    return;
  }

  if (run->GetRunNo() == -1) {
    Abort("SetMagnetField", "Run number is -1, check it out");
    return;
  }
  
  // get magnet setting
  fMagnetVolume->SetCurrent(run->GetD5Set(), run->GetD5Pol());

  if (Verbose() > 10) 
    cout << " *** BrMrsTrackingModule::SetMagnetField() -- > D5:  " 
	 << fMagnetVolume->GetField() << endl;
}

//___________________________________________________________________________
 void BrMrsTrackingModule::Begin()
{
  // Set magnet field and pick up the platform angle
  // this is done at begin time since one might scan sequences 
  // from diff. runs with diff. settings (maybe not a good thing to do)
  if(fUseDb)
    SetMagnetField();
}

//___________________________________________________________________________
 void BrMrsTrackingModule::Event(BrEventNode* InputNode, 
				BrEventNode* OutputNode)
{
  //
  //  Normal Event member function for MRS tracking
  // 
  //  Requires: Two instances of BrDataTable named "DetectorTrack TPM1" and 
  //            "DetectorTrack TPM2" in either InputNode or OutputNode. 
  //            These tables have to be filled with instances of BrDetectorTrack
  //            from their respective TPCs.
  //  Result:   BrDataTable "MrsTracks" added to outputnode.
  //            The elements of this table are BrMrsTrack.
  //
  //  The module has the following steps:
  //  a) Combine Front (TPM1) and Back (TPM2) tracks
  //  b) Track combined tracks back to nominal vertex
  //  c) Check goodness of tracks and perform cleanup
  //     of ghost tracks, i.e. cases where the same track
  //     are used several times. The tracks are not removed
  //     but are marked as 'ghosts'.
  //  d) Add tracks to OutputNode (good tracks and potential
  //     ghosts if requested). 
  //


  SetState(kEvent);
  SetStatus(kOk);

  fCombineModule->Clear();

  // --------------------------------------------
  // --        Combining local tracks          --
  // --------------------------------------------

  BrEventNode* node = InputNode;
  // YACK!!
  if (!InputNode->GetDataTable("DetectorTrack TPM1") ||
      !InputNode->GetDataTable("DetectorTrack TPM2"))
    node = OutputNode;

  fCombineModule->CombineFrontBack(node);
  
  if(DebugLevel() > 2)
    cout << "<BrMrsTrackingModule::Event>: Get array of matched tracks" << endl;
  
  // ----------------------------------------
  // --        Get matched tracks          --
  // ----------------------------------------

  fTracksFrontBack  = fCombineModule->GetMatchedTracks();
    
  if(DebugLevel() > 1)
    cout << "<BrMrsTrackingModule::Event>: Number of matched tracks = "
	 << fTracksFrontBack->GetEntries() << endl;

  if (!fTracksFrontBack->GetEntries())
    return;

  // ----------------------------------------
  // --        Prepare output table        --
  // ----------------------------------------

  BrDataTable* mrstable = new BrDataTable("MrsTracks");
  OutputNode->AddObject(mrstable);

  // ----------------------------------------
  // --    Track down to primary vertex    --
  // ----------------------------------------

  TrackToVertex(mrstable);
  
  // -----------------------------------------------------------------------------
  // --  Pick up BB vertex (the new one) if required                            --
  // --  (Note the increasing usage of the check both in input and output node  --
  // --  for data objects. This shows the wrong usage of the event model          --
  // --  and bratmain. No good.)                                                --
  // -----------------------------------------------------------------------------

  if (HistOn()) {
    BrBbVertex* bbVtx = (BrBbVertex*)InputNode->GetObject("BB Vertex");
    if (!bbVtx)
      bbVtx = (BrBbVertex*)OutputNode->GetObject("BB Vertex");
    
    // check only with small tubes vtx
    // well maybe not for pp !!
    if (bbVtx)
      if (bbVtx->GetMethod() == 2)
	TrackToBbVertex(mrstable, bbVtx);
  }

  // ----------------------------------------
  // --     Do some selected print out     --
  // ----------------------------------------
  
  for (Int_t i = 0; i < mrstable->GetEntries(); i++) {
    if (Verbose() > 25)
      ((BrMrsTrack*)mrstable->At(i))->Print("A");
    else if (Verbose() > 20)
      ((BrMrsTrack*)mrstable->At(i))->Print("PD");
    else if (Verbose() > 15)
      ((BrMrsTrack*)mrstable->At(i))->Print("P");
    else if (Verbose() > 10)
      ((BrMrsTrack*)mrstable->At(i))->Print("");
  }


}
//___________________________________________________________________
 void BrMrsTrackingModule::Finish()
{
  // Finisn entry method
  // Call combine module finish
  fCombineModule->Finish();
}
//___________________________________________________________________
 void BrMrsTrackingModule::TrackToVertex(BrDataTable* tab)
{
  // This routine will loop over all tracks found using their position
  // in T1 as a starting point for backward tracking.
  // The target Tx ,Ty and Tz are in the IP coordinate system ie at
  // x = 0
  // Calculate scattering angles from vertex and position
  // in tpm1. This should be considered a first approximation for
  // vertex , and particle theta, phi. It should be redone when
  // checked against vertex positions determined by the vertex modules.
  // It also calculate the partial tracklength from the vertex to TPM1 
  // position.
  // The vertex assumed is that of each individual track not any global
  // vertexing; this should be done later.
  // Some diagnostics plots are also filled.
  // The front back tracks are added to table at this point.

  if(DebugLevel() > 0)
    cout << " TrackToVertex: " << endl
	 << "  #Tracks = "     << fTracksFrontBack->GetEntries() << endl;

  //
  // Setup the beamline plane use to project each track
  // to find the tvertex for this track.
  //

  BrPlane3D beamPlane(0,0,0, 0,0,1, 0,1,0);

  const BrCoordinateSystem* tpm1system = fFrontVolume->GetCoordinateSystem();

  Int_t nmatch = fTracksFrontBack->GetEntries();
  Int_t nGood  = 0;

  for(Int_t it = 0; it < nmatch; it++) {
    // pick up matched track
    BrMatchedTrack *matchTrk = (BrMatchedTrack*) fTracksFrontBack->At(it);

    if (HistOn())
      hStatus->Fill(matchTrk->GetStatus());

    // check ghostbusting status
    // BUT do not reject tracks !!!
    //
    if (matchTrk->GetStatus() == 1)
      nGood++;

    BrDetectorTrack *frTrk = matchTrk->GetFrontTrack();
 
    BrVector3D tpm1position = tpm1system->TransformToMaster(frTrk->GetPos());
    BrLine3D frontGlobalTrack(tpm1position,
			      tpm1system->RotateToMaster(frTrk->GetAlpha()));
    
    BrVector3D vertex = beamPlane.GetIntersectionWithLine(frontGlobalTrack);
    
    BrVector3D trackdir = tpm1position - vertex;
    Double_t theta = trackdir.Theta();
    Double_t phi   = trackdir.Phi();

    if(HistOn()) {
      hTargetTx->Fill(vertex(0));
      hTargetTy->Fill(vertex(1));
      hTargetTz->Fill(vertex(2));
      hPhi->Fill(phi);
      hTheta->Fill(theta);
    }
    
    //   
    // Create a new MRS track
    
    BrMrsTrack* mrstrack = new BrMrsTrack();
    
    mrstrack->SetTrackId(it);
    mrstrack->SetMatchedTrack(matchTrk);
    mrstrack->SetTrackVertex(vertex);
    mrstrack->SetTheta(theta);
    mrstrack->SetPhi(phi);

    // check now the other end (TOFW stuff)
    // Similar to previous trackToTofw
    Float_t length = PartialPath(mrstrack);

    mrstrack->SetPartialPath(length);
    mrstrack->SetPathLength(length + trackdir.Norm());

    tab->Add(mrstrack);

    if (HistOn()) {
      hPathLength->Fill(mrstrack->GetPathLength());
      hPartialPath->Fill(mrstrack->GetPartialPath());
      hP->Fill(mrstrack->GetMomentum());
    }
  }

  if (HistOn())
    hGhost->Fill(nmatch, nGood);
}

//___________________________________________________________________
 void BrMrsTrackingModule::TrackToBbVertex(BrDataTable* tab,
					  BrBbVertex* bbVtx)
{
  // protected
  // create plane centered on Z0 and normal to TPM1 front plane
  // in global coord.!!
  // (fv) Why! The relevant resolution of the BbVertex is along the
  // Z axis so differences should be measured along this axis. On the
  // normal plane the resolution goes as sigma(z)/sin(mrsangle)
  //
  BrVector3D normal(fFrontVolume->GetFrontPlane().GetA(),
		    fFrontVolume->GetFrontPlane().GetB(),
		    fFrontVolume->GetFrontPlane().GetC());

  BrVector3D beamDir(0, 0, 1);
  BrVector3D xDir(1, 0, 0);
  BrVector3D bb(0, 0, bbVtx->GetZ0());
  BrVector3D ip(0, 0, 0);

  // BB vertex in TPM1 front plane coord.
  BrVector3D bbLocal;
  fFrontVolume->GlobalToLocal(bb, bbLocal, 0);

  BrPlane3D bbPlane(bb, normal); // parallel to TPM1 front plane
  BrPlane3D longPlane(ip, xDir); // sits at X (global) = 0
  
  for (Int_t t = 0; t < tab->GetEntries(); t++) {
    BrMrsTrack* mrstrack = (BrMrsTrack*)tab->At(t);
    BrDetectorTrack* trk = mrstrack->GetFrontTrack();

    // project track line on BB plane normal to MRS axis    
    BrLine3D trkLine(fFrontVolume->LocalToGlobal(trk->GetTrackLine()));
    BrVector3D inters, gInters;
    gInters = bbPlane.GetIntersectionWithLine(trkLine);

    // calculate distance from intersection to this point:
    Float_t dist = (bb - gInters).Norm();
    hTrkBbDist->Fill(dist);

    // convert into plane coord. system 
    fFrontVolume->
      GlobalToLocal(gInters, inters, 0);

    hTrkBbXY->Fill(inters(0) - bbLocal(0), inters(1));

    // projection on plane normal to global X axis
    BrVector3D proj = longPlane.GetIntersectionWithLine(trkLine);
    hTrkBbZ->Fill(bb(2) - proj(2));
  }
}

//___________________________________________________________________
 Float_t BrMrsTrackingModule::PartialPath(BrMrsTrack* track)
{
  // Projects to tofw and
  // calculate the length from the track position in TPM1 
  // to the time of flight detector.
  // If the track does not point within any panel A slatNo of 0 isassigned
  // and the position is undefined.

  if(DebugLevel() > 0)
    cout << " In PartialLength() " << endl;

  // initial stuff
  BrDetectorTrack* dbtrack = track->GetBackTrack();
  BrDetectorTrack* dftrack = track->GetFrontTrack();
    
  Float_t length = 0;
  Int_t     slat = 0;
  Int_t    panel = -1;

  BrVector3D local;    
  BrVector3D trkProj;
  // -------------------
  
  for (Int_t panelNo = 0; panelNo < fTofParams->GetNoPanels(); panelNo++) {
    BrPlane3D frontPanel(fPanelVol[panelNo]->GetFrontPlane()); 
    BrLine3D
      trkLine(fBackVolume->LocalToGlobal(dbtrack->GetTrackLine()));
    
    trkProj = frontPanel.GetIntersectionWithLine(trkLine);
    fPanelVol[panelNo]->GlobalToLocal(trkProj, local, 0);

    if(TMath::Abs(local(0)) < fPanelVol[panelNo]->GetSize()[0]/2. &&
       TMath::Abs(local(1)) < fPanelVol[panelNo]->GetSize()[1]/2.) {
	panel = panelNo;
	break;
    }
  }

  if (panel < 0) {
    BrVector3D wrong(-999, -999, -999);
    track->SetPointedSlat(0);
    track->SetPointedPanel(panel);
    track->SetProjOnTof(wrong);
    return 0;
  }
  
  slat = fPanelPar[panel]->GetSlatNo(local(0));
  
  if (panel > 0)
    for(Int_t j = 0; j < panel;j ++)
      slat += fPanelPar[j]->GetNoSlats();
      
  // get matched track length (from TPM1 to TPM2 track pos)
  Double_t mtLength = track->GetMatchedTrack()->GetMatchedTrackLength(); 
      
  // back track position in TPM2 (global c.s.)
  BrVector3D bkPos;
  fBackVolume->LocalToGlobal(dbtrack->GetPos(), bkPos, 0);
      
  // TPM2 to TOFW
  BrVector3D backTrackToTof(bkPos - trkProj);
  length =  mtLength + backTrackToTof.Norm();

  track->SetPointedPanel(panel);
  track->SetPointedSlat(slat);
  track->SetProjOnTof(trkProj); // global coord. sys.

  if (HistOn())
    hSlatPos->Fill(slat);

  return length;
}


//_________________________________________________________________________
//
// $Log: BrMrsTrackingModule.cxx,v $
// Revision 1.10  2002/08/30 17:36:34  videbaek
// Corrected the orderof verbose selection 10,15,20,25..
//
// Revision 1.9  2002/06/13 00:42:22  videbaek
// Change pathlength in histogram
//
// Revision 1.8  2002/02/26 08:58:20  cholm
// Added const'ness to many methods - needed by ISO C++
//
// Revision 1.7  2002/01/27 21:15:31  videbaek
// Remove references to hTriggers. This was only used by me for
// so checking (could and should have been done outside this module).
//
// Revision 1.6  2002/01/25 17:03:33  videbaek
// Take care of errors from combine module
//
// Revision 1.5  2001/12/07 21:09:22  videbaek
//   Add SetUseDb()  method to control if magnet filed is picked from Db or Not.
// default is to do so, but there is a needed particular for MC studies not to use
// the rundatabase. This is a lesson for other modules too.
//
// Added finish method()
//
// Revision 1.4  2001/12/06 19:53:58  jens
// Updated BrMrsTrackingModule::Event to look for the TPC track tables in the output node if not found in the input node. 
// Also cleaned up this function a little bit.
//
// Revision 1.3  2001/11/19 15:07:40  ouerdane
// removed useless member fMrsAngle
//
// Revision 1.2  2001/11/08 21:39:47  videbaek
// Do not skip adding tracks even if ghost- should not be taken out-
// Just count good. Track removal should be in later stage of analysis.
//
// Revision 1.1  2001/11/05 06:42:21  ouerdane
// added new classes replacing older spectrometer track modules
//



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