//$Id: BrTPMClusterVertexModule.cxx,v 1.3 2002/01/03 19:53:54 cholm Exp $

/////////////////////////////////////////////////////////////////////////////
// BrTPMClusterVertexModule
//
// Vertexing module for the MRS, using only clusters from  MTPC1 
// May 2000
// bjornhs
// 
//
// Requires: a minimum number of clusters in MTPC1
//
// This vertexing-method assumes x=0, then combines all pairs of clusters
// in TPM1 and checks where the line between them intersects the yz-plane.
// The intersections are histogrammed, and the y,z-coord is estmated by 
// zooming in on and gauss-fitting this histogram
//
// Usage:
//
// 1) Do BrTPMClusterVertexModule* vertex = new BrTPMClusterVertexModule("TPM1","VertexFinder");
// 2) vertex->Event(inputnode,outputnode) puts a BrVertex-object
//    (BrDataObject) with name TPM1 ClusterVtxData in outputnode.
//
///////////////////////////////////////////////////////////////////////////////

#include <BrTPMClusterVertexModule.h>
#ifndef BRAT_BrTableNames
#include "BrTableNames.h"
#endif

ClassImp(BrTPMClusterVertexModule)

 BrTPMClusterVertexModule::BrTPMClusterVertexModule():BrModule() 
// Default Constructor.
{
  SetState(kSetup);  
  fFrontVolume = 0;
  fAllHits = 0;
  fHitSlice = 0;
}

 BrTPMClusterVertexModule::BrTPMClusterVertexModule(Text_t *Name,Char_t *Title)
  :BrModule(Name,Title)
{
  // Use this creator
  SetState(kSetup);  
  fFrontVolume = 0;
  fAllHits = 0;
  fHitSlice = 0;
  SetDefaultParameters();

}

 void BrTPMClusterVertexModule::SetDefaultParameters() 
{
  // Set/Reset internal vertexing parameters to default values.
  fNoSlices = 16;
  fMinNoHits = 50;
  fMaxNoHits = 1000;
  fNHistBins = 50;
  fDraw_Histogram = 0;
  fUseHistoZoomY = 0;
  fUseHistoZoomZ = 1;

  kRemoveBackgrounds = kFALSE;

  fZpos_RMS = 5;
  fYpos_RMS = 2;

  fYZoomInterval = 10;
  fZZoomInterval = 20;

  fHistbinsDivisionFactor = 1;

}


 BrTPMClusterVertexModule::~BrTPMClusterVertexModule() //Destructor
{
  //if(fFrontVolume) delete fFrontVolume;
  if (fHitSlice) {
    for (Int_t j = 0 ; j < fNoSlices ; j++) {
      fHitSlice[j].Clear();
    }
    delete [] fHitSlice;
  }

  if (fAllHits){
    fAllHits->Delete();
    delete fAllHits;
  }
}

 void BrTPMClusterVertexModule::DefineHistograms()
{
  //
  // define diagnostics histograms
  // Book histograms (called from BrModule::Book)

  hVertex_z = new TH1F("hVertex_z", "Vertex position in z", 800, -40, 40);
  hVertex_sigma = new TH1F("hVertex_sigma", "Sigma of det. histogram", 100, 0, 10);

}

 void BrTPMClusterVertexModule::Init()
{
  // Get Parameters from databases (geometry and parameters)
  //
  // Run this before event!

  SetState(kInit);

  BrGeometryDbManager* gGeomDb = BrGeometryDbManager::Instance();
  fFrontVolume = (BrDetectorVolume*)gGeomDb->GetDetectorVolume("BrDetectorVolume", GetName());

  BrParameterDbManager* gParamDb = BrParameterDbManager::Instance();
  fParams_p = (BrDetectorParamsTPC*)gParamDb->GetDetectorParameters("BrDetectorParamsTPC", GetName());

  
  // Vertex not found yet!
  //  fVertex_found = 0;

  // Hit arrays
  fAllHits = new TObjArray();
  fHitSlice = new TObjArray[fNoSlices];

  // Histograms
  //  ypos = new TH1F("ypos", "y-coord of vertex", fNHistBins, -20, 20);
  //  zpos = new TH1F("zpos", "z-coord of vertex", fNHistBins, -300, 300);

}


 void BrTPMClusterVertexModule::Event(BrEventNode* InputNode, BrEventNode* OutputNode)
{
  //
  // Take clusters from TPM1, find primary vertex
  //

  SetState(kEvent);
  fVertex_found = 0;

  if(DebugLevel()>2) cout<< "Histogram/cluster-based vertexing routine***" << endl;

  if(!fFrontVolume){
    Error("Event","Module not properly setup");
    return;
  }

  BrDataTable* tpm1clusters = InputNode->GetDataTable(Form("%s %s", BRTABLENAMES kTpcHit, GetName()));

  if (!tpm1clusters) {
    tpm1clusters = OutputNode->GetDataTable(Form("%s %s", BRTABLENAMES kTpcHit, GetName()));
    if (!tpm1clusters) {
      Warning("Event","BrDataTable %s %s not found in this event!",BRTABLENAMES kTpcHit, GetName());
      return;
    }

  } 

  if(tpm1clusters) {

    BrVertex* vtxdata = new BrVertex("TPM1 Cluster VtxData", "TPM1 Cluster Vertex Data Object");

    FindVertex(tpm1clusters,vtxdata);

    OutputNode->AddObject(vtxdata);
    if (DebugLevel()) {
      OutputNode->ListObjects();
      cout << " Vertex z: " << vtxdata->GetZ() << " VarianceZ: " << vtxdata->GetVarianceZ() << " VertexZSigma: " << vtxdata->GetVertexZSigma() << endl;
      cout << " Vertex y: " << vtxdata->GetY() << " VarianceY: " << vtxdata->GetVarianceY() << " VertexYSigma: " << vtxdata->GetVertexYSigma() << endl;

    }
  }

}

 void BrTPMClusterVertexModule::FindVertex(BrDataTable* tpm1clusters, BrVertex* vtxdata)
{
  // This does the work - called by Event()
  
  if (DebugLevel()>2) cout << "TPM1 hits > " << tpm1clusters->GetEntries() << endl;
  ClustersToHits(tpm1clusters); 
  // Go from BrTpcHits to BrVector3Ds

  Float_t Vertex_loc_y;               //  Vertex y(cm)
  Float_t Vertex_loc_z;               //  Vertex z(cm)
  Float_t Vertex_fit_sigma_y;         // Sigma of gaussian used to find y
  Float_t Vertex_fit_sigma_z;         // Sigma of gaussian used to find z
  Float_t VarianceZ;                  // Error in z-pos = sqrt(VarianceZ)
  Float_t VarianceY;                  // Error in y-pos = sqrt(VarianceY)

  // Transform to global coords
  Int_t number_of_hits = 0;
  BrVector3D* hit_b = 0;
  Float_t h_l[3];
  Float_t h_g[3];
  for (Int_t i=0; i<fAllHits->GetLast()+1; i++){
    hit_b = (BrVector3D*)fAllHits->At(i);
    h_l[0]=hit_b->GetX();
    h_l[1]=hit_b->GetY();
    h_l[2]=hit_b->GetZ();
    fFrontVolume->LocalToGlobal(h_l, h_g, 0);
    Float_t yoverx = h_g[1]/h_g[0];
    for (Int_t j = 0 ; j < fNoSlices ; j++) {
      Float_t lower_lim = -0.08+j*(0.16/(fNoSlices));
      Float_t upper_lim = -0.08+(j+1)*(0.16/(fNoSlices));
      if (yoverx < upper_lim & yoverx > lower_lim) {
	fHitSlice[j].Add(hit_b);
	number_of_hits++;
      }
    }
  }

  ypos = new TH1F("ypos", "y-coord of vertex", fNHistBins, -20, 20);
  zpos = new TH1F("zpos", "z-coord of vertex", fNHistBins, -300, 300);

  // Test
  /*
  TH1F* ypos = new TH1F("ypos", "y-coord of vertex", fNHistBins, -20, 20);
  TH1F* zpos = new TH1F("zpos", "z-coord of vertex", fNHistBins, -300, 300);

  TH1F* zpos_orig = 0;
  TH1F* ypos_orig = 0;
  */
  // Test

  //
  // Test if there are enough clusters to work with
  //
  if (number_of_hits > fMinNoHits && number_of_hits < fMaxNoHits) {

    // Combine hits, project to yz-plane, increment histograms
    for (Int_t j = 0 ; j < fNoSlices ; j++) {
      CombineHits(fHitSlice[j], zpos, ypos);
    }

    if (fDraw_Histogram) {
      zpos_orig = (TH1F*)zpos->Clone();
      ypos_orig = (TH1F*)ypos->Clone();
    }

    if (kRemoveBackgrounds) {
      // Remove linear backgrounds
      TH1F *zpos_subtr = (TH1F*)RemoveLinearBackground(zpos);
      TH1F *ypos_subtr = (TH1F*)RemoveLinearBackground(ypos);
      zpos->Delete();
      ypos->Delete();
      zpos = zpos_subtr;
      ypos = ypos_subtr;     
    }

    // Find the max-value-points of zpos and ypos 
    Float_t zpos_mean = GetMaxPoint(zpos);
    Float_t ypos_mean = GetMaxPoint(ypos);

    if (fUseHistoZoomY || fUseHistoZoomZ) {
      // Zoom in on the max-points to get a better determination
      TH1F *zpos3 = new TH1F("zpos3","zpos3",fNHistBins/2,zpos_mean-fZZoomInterval/2,zpos_mean+fZZoomInterval/2);
      TH1F *ypos3 = new TH1F("ypos3","ypos3",fNHistBins/2,ypos_mean-fYZoomInterval/2,ypos_mean+fYZoomInterval/2);

      // Combine hits, increment zoomed histograms
      for (Int_t j = 0 ; j < fNoSlices ; j++) {
	CombineHits(fHitSlice[j], zpos3, ypos3); 
      }    

      if (fUseHistoZoomZ) {
	zpos->Delete();
	zpos=zpos3;
      }
      if (fUseHistoZoomY) {
      ypos->Delete();
      ypos=ypos3;
      }

      if (!fUseHistoZoomY) delete ypos3;
      if (!fUseHistoZoomZ) delete zpos3;

      if (kRemoveBackgrounds) {
	// Remove linear backgrounds
	TH1F *zpos_subtr_z = (TH1F*)RemoveLinearBackground(zpos);
	TH1F *ypos_subtr_z = (TH1F*)RemoveLinearBackground(ypos);
	zpos->Delete();
	ypos->Delete();
	zpos = zpos_subtr_z;
	ypos = ypos_subtr_z;
      }

      // Find the max-value-points of zpos and ypos 
      if (fUseHistoZoomZ) zpos_mean = GetMaxPoint(zpos);
      if (fUseHistoZoomY) ypos_mean = GetMaxPoint(ypos);

    }

    if (DebugLevel()>1) {
      cout << "Y maxpoint: " << ypos_mean << endl;
      cout << "Z maxpoint: " << zpos_mean << endl;
    }

    // Do a gaussian fit of the z-histogram
    TF1 *f1 = new TF1("f1","gaus",zpos_mean-(fZpos_RMS),zpos_mean+(fZpos_RMS));
    f1->SetParameters(0,zpos_mean,1);
    f1->SetParLimits(1,zpos_mean-fZpos_RMS,zpos_mean+fZpos_RMS);
    f1->SetParLimits(2,0,5);
    zpos->Fit("f1","q0r");
    if (DebugLevel()>2) 
      {
	cout << "Mean from hist > " << zpos_mean << endl;
	cout << "Mean from fit  > " << f1->GetParameter(1) << endl;
	cout << "Sigma from fit  > " << f1->GetParameter(2) << endl;
      }
    if (f1->GetParameter(1)!=0 && f1->GetParameter(2)!=0) {
      Vertex_loc_z = f1->GetParameter(1);
      Vertex_fit_sigma_z = f1->GetParameter(2);
      VarianceZ = FindVariance(zpos, Vertex_loc_z, Vertex_fit_sigma_z);
    }
    else {
      Vertex_loc_z = 9999;
      Vertex_fit_sigma_z = 9999;
      VarianceZ = 9999;
    }

    // Do a gaussian fit of the y-histogram
    TF1 *f2 = new TF1("f2","gaus"); 
    //,ypos_mean-(ypos_RMS/2),ypos_mean+(ypos_RMS/2));
    ypos->Fit("f2","q0");
    if(DebugLevel()>2){
      cout << "Mean from hist > " << ypos_mean << endl;
      cout << "Mean from fit  > " << f2->GetParameter(1) << endl;
      cout << "Sigma from fit  > " << f2->GetParameter(2) << endl;
    }
    Vertex_loc_y = f2->GetParameter(1);
    Vertex_fit_sigma_y = f2->GetParameter(2);
    VarianceY = FindVariance(ypos, Vertex_loc_y, Vertex_fit_sigma_y);

    if (fDraw_Histogram) {
      TCanvas *det_hist = new TCanvas("det_hist","Vert.det. histogram",0,0,500,500);
      det_hist->Divide(2,2);
      det_hist->cd(1);
      ypos_orig->Draw();
      det_hist->cd(2);
      zpos_orig->Draw();
      det_hist->cd(3);
      ypos->Draw();
      det_hist->cd(4);
      zpos->Draw();
    }

    // Is vertex found?
    if (f1->GetParameter(2) != 0 && f1->GetParameter(2)!= 0) fVertex_found = 1;


    // Delete fitting-functions 
    delete f1; 
    delete f2;  

    if(fHistOn){
      hVertex_z->Fill(Vertex_loc_z);
      hVertex_sigma->Fill(Vertex_fit_sigma_z);
    }

  }
  else {
    if (DebugLevel()>2) cout << "Too few or too many hits in TPM1" << endl;
    Vertex_loc_z = 999; // Set vertex to something weird if it isn't located   
    Vertex_loc_y = 999;
    fVertex_found = 0;
  }

  // Fill vertex object
  vtxdata->SetX(0);
  vtxdata->SetY(Vertex_loc_y);
  vtxdata->SetZ(Vertex_loc_z); 
  vtxdata->SetVertexYSigma(Vertex_fit_sigma_y);
  vtxdata->SetVertexZSigma(Vertex_fit_sigma_z);
  vtxdata->SetVertexVarianceY(VarianceY);
  vtxdata->SetVertexVarianceZ(VarianceZ);
  vtxdata->SetVertexFound(fVertex_found);
  vtxdata->SetVertexMethodUsed(1);

  // Cleanup
  if (fDraw_Histogram==0){
    delete zpos;
    delete ypos;
  }
  for (Int_t j = 0 ; j < fNoSlices ; j++) {
    fHitSlice[j].Clear();
  }

  fAllHits->Delete();
  
}


 void BrTPMClusterVertexModule::ListVertexParameters()
{
  // List cut values used in vertexing.
  cout << "Parameters in BrTPMClusterVertexModule:" << endl;
  cout << "Number of TPC-slices: " << fNoSlices << endl;
  cout << "Minimum number of hits in MTPC1: " << fMinNoHits << endl;
  cout << "Maximum number of hits to use: " << fMaxNoHits << endl;
  cout << "Number of bins in the determinig histogram: " << fNHistBins << endl;
  
}

 Float_t BrTPMClusterVertexModule::GetMaxPoint(TH1 *hist)
{
  // Returns the point where hist has its maximum 
  Float_t max_value = 0; 
  Float_t bin_value;
  Int_t max_bin;
  {
    for (Int_t i = 1; i <= hist->GetNbinsX(); i++) {
      bin_value = hist->GetBinContent(i);
      if (bin_value >= max_value) {
	max_value = bin_value;
	max_bin = i;
      }
    }
  }    
  Float_t hist_max = hist->GetBinCenter(max_bin);
  return hist_max;
}

 Float_t BrTPMClusterVertexModule::FindVariance(TH1 *hist, Float_t mean, Float_t sigma)
{
    Float_t min_lim =  mean-3*sigma;
    Float_t max_lim =  mean+3*sigma;
    Float_t sum_entries;
    {
      for (Int_t i = 1;i<=hist->GetNbinsX();i++) {
	if (hist->GetBinCenter(i)>min_lim && hist->GetBinCenter(i)<max_lim) 
	  sum_entries = sum_entries+hist->GetBinContent(i);
      }
    }
    Float_t Variance = (sigma*sigma)/sum_entries;
    return Variance;
}

 TH1 *BrTPMClusterVertexModule::RemoveLinearBackground(TH1 *hist)
{
  // Removes a linear background from hist
  TAxis *tx1 = hist->GetXaxis();
  TF1 *linear = new TF1("linear","[0]*x+[1]",tx1->GetXmin(),tx1->GetXmax());  
  hist->Fit("linear","q0");
  TH1F *hist2=(TH1F*)hist->Clone();
  hist2->Reset();
  for (Int_t bin = 1;bin<=hist->GetNbinsX();bin++) {
    Float_t x = hist->GetBinCenter(bin);
    Double_t fval = linear->Eval(x);
    Double_t d = hist->GetBinContent(bin)-fval;
    if (d>0) hist2->Fill(x,d);
  } 
  delete linear;
  return hist2;
}

 void BrTPMClusterVertexModule::ClustersToHits(BrDataTable* tpm1clusters)
{
  // Converts the middle position of BrDetectorHits (clusters) in tpm1clusters to
  // BrVector3D's in allhits

  BrTpcHit* dethit = 0;
  //  const Float_t *hitloc;

  for (Int_t i=0;i<tpm1clusters->GetEntries();i++) {
    dethit = (BrTpcHit*)tpm1clusters->At(i);
    //    hitloc = dethit->GetPos();
    //    BrVector3D* hit = new BrVector3D(hitloc[0], hitloc[1], hitloc[2]); 
    BrVector3D* hit = new BrVector3D(dethit->GetPos()); 
    fAllHits->Add(hit);
  }

} 

 void BrTPMClusterVertexModule::CombineHits(TObjArray& hits, TH1F* histogram_z, TH1F* histogram_y)
{
  // Combines two and two clusters, projects the line between them back to x=0
  // (global), and increments the histogram histogram_z
  TObjArray* FrontHits = new TObjArray;
  TObjArray* BackHits = new TObjArray;
  Float_t b_l[3], f_l[3],b_g[3], f_g[3];
  Float_t a, ay, y, z;
  BrVector3D* hit;
  BrVector3D* hit_b;  
  BrVector3D* hit_f;
  {
    for (Int_t i=0; i<hits.GetLast()+1; i++){
      hit = (BrVector3D*)hits.At(i);
      if (hit->GetZ() < 0) FrontHits->Add(hit); 
      // These lines decide where to divide
      else if (hit->GetZ() >= 0) BackHits->Add(hit); // the TPC
    }
  }

  for (Int_t i=0; i<BackHits->GetLast()+1; i++){
    hit_b = (BrVector3D*)BackHits->At(i);
    b_l[0] = hit_b->GetX();
    b_l[1] = hit_b->GetY();
    b_l[2] = hit_b->GetZ();
    fFrontVolume->LocalToGlobal(b_l, b_g, 0); 
    for (Int_t j=0; j<FrontHits->GetLast()+1; j++){
      hit_f = (BrVector3D*)FrontHits->At(j);
      f_l[0] = hit_f->GetX();
      f_l[1] = hit_f->GetY();
      f_l[2] = hit_f->GetZ();
      fFrontVolume->LocalToGlobal(f_l, f_g, 0);
      a = (b_g[0]-f_g[0])/(b_g[2]-f_g[2]);
      ay = (b_g[1]-f_g[1])/(b_g[0]-f_g[0]);
      y = f_g[1]-f_g[0]*ay;
      z = f_g[2]-(f_g[0]/a);
      histogram_z->Fill(z);
      histogram_y->Fill(y);      
    }
  }
  delete FrontHits;
  delete BackHits;
}

//$Log: BrTPMClusterVertexModule.cxx,v $
//Revision 1.3  2002/01/03 19:53:54  cholm
//Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names
//
//Revision 1.2  2001/09/06 19:30:47  trine
//Revised BrTPMClusterModule to work with new TPC code and table names,
//and to run with bratmain.
//Class version changed to 0 as proper for a module.
//Small general code cleanup that should not alter functionality.
//
//Revision 1.1.1.1  2001/06/21 14:55:14  hagel
//Initial revision of brat2
//
//Revision 1.11  2001/01/30 10:36:58  pchristi
//Added some const pointers to reflect changes to const methods in BrDetectorHit.
//
//Revision 1.10  2000/11/22 21:21:25  videbaek
//Changes due to better handling of Init().
//
//Revision 1.9  2000/11/07 16:40:23  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.8  2000/09/21 20:52:10  cholm
//Several changes so as it may compile on Windoze.
//
//Revision 1.7  2000/09/15 12:32:51  pchristi
//Removed delete of non owned object in destructor.
//
//Revision 1.6  2000/09/07 08:51:46  bjornhs
//Added Y-uncertainties. Removed some unnecesary CPU-usage.
//
//Revision 1.5  2000/09/04 16:51:37  videbaek
//Add some redundent brackets to make WIN32 happy.
//
//Revision 1.4  2000/08/30 09:17:21  bjornhs
//Nasty mem-leaks removed. Thanks, Trine.
//
//Revision 1.3  2000/07/15 14:09:03  oslo
//Added prelim. VarianceZ() to BrTPM1ClusterVertexModule
//
//Revision 1.2  2000/06/27 10:07:46  bjornhs
//Added y-determination, plus minor bugfixes
//
//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.4  2000/05/12 07:06:10  bjornhs
//Event is now Event(BrEventNode* Inputnode, BrEventNode* outputnode);
//Adds BrVertex-object named "VtxData" to outputnode
//Fixed minor bugs
//

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

Validate HTML
Validate CSS