HoughTransNCPi0.cxx

Go to the documentation of this file.

//
// HoughTransCB
//
// Mark's hough transform.
//
// Created:  M. Dorman -- May, 2005
//
// $Author: dja25 $ 
//
// $Revision: 1.3 $
// 
// $Name:  $
//
// $Id: HoughTransNCPi0.cxx,v 1.3 2008/04/14 15:03:49 dja25 Exp $
//

#include "HoughTransNCPi0.h"
#include <iostream>
#include <cmath>
#include "TCanvas.h"

HoughTransNCPi0::HoughTransNCPi0()//:houghSpace("houghSpace","houghSpace",80,-2.0,2.0,80,-4.0,4.0),vtxz(0)
{
  houghSpace = new TH2F("houghSpace","houghSpace",80,-2.0,2.0,80,-4.0,4.0);
  houghSpace->SetDirectory(0);
  vtxz = 0.0;
}

HoughTransNCPi0::~HoughTransNCPi0(){
  // Mark, you need to delete histograms created with new!!!
  delete houghSpace; houghSpace=0;

}

HoughTransNCPi0::HoughTransNCPi0(Int_t numSamps,Float_t lowGrad,Float_t highGrad)//:
{
  
  houghSpace = new TH2F("houghSpace","houghSpace",numSamps,lowGrad,highGrad,numSamps,-4.0,4.0);
  houghSpace->SetDirectory(0);
  vtxz = 0.0;
}

HoughTransNCPi0::HoughTransNCPi0(Int_t numSamps,Float_t lowGrad,Float_t highGrad,Float_t lowIcept,Float_t highIcept)//:
{
  houghSpace = new TH2F("houghSpace","houghSpace",numSamps,lowGrad,highGrad,numSamps,lowIcept,highIcept);
  houghSpace->SetDirectory(0);
  vtxz = 0.0;
}

void HoughTransNCPi0::SetVtxz(Float_t vertz)
{
  vtxz=vertz;
}

Float_t HoughTransNCPi0::GetVtxz()
{
  return vtxz;
}

void HoughTransNCPi0::FillHough(NtpSRStrip* stp)
{
  FillHough(stp->z,stp->tpos);
}

void HoughTransNCPi0::FillHough(Float_t zPos, Float_t tPos)
{       
  if(zPos<(vtxz)) return;
  for(Int_t i=1;i<=houghSpace->GetNbinsX();++i){
    Float_t cos=(houghSpace->GetXaxis())->GetBinCenter(i);
    Float_t icept=tPos-sqrt(pow((1/cos),2)-1)*(zPos-vtxz);
    if(icept>4.0 || icept<-4.0) continue; //max MINOS transverse dimensions
    houghSpace->Fill(cos,icept);
    // cout << "Filling Hough Space with grad " << gradient << " and icept " << icept << endl;
  }
}
void HoughTransNCPi0::FillHoughweight(Float_t zPos, Float_t tPos, Float_t energy)
{
  if(zPos<vtxz) return;
  for(Int_t i=1;i<houghSpace->GetNbinsX();++i){
    Float_t gradient=(houghSpace->GetXaxis())->GetBinCenter(i);
    Float_t icept=tPos-gradient*(zPos-vtxz);
    Float_t weight= energy;
    if(icept>4.0 || icept<-4.0) continue; //max MINOS transverse dimensions
    houghSpace->Fill(gradient,icept,weight);
    // cout << "Filling Hough Space with grad " << gradient << " and icept " << icept << endl;
  }
}
void HoughTransNCPi0::DrawHough()
{
  //TFile f("hough.root",new");
  houghSpace->Draw("LEGO2");
  //  houghSpace->Write();
}

Int_t HoughTransNCPi0::EntriesHough()
{
  Int_t counts = static_cast<Int_t>(houghSpace->Integral());
  return counts;
}

void HoughTransNCPi0::ResetHough()
{
  houghSpace->Reset();
}

void HoughTransNCPi0::DeleteHough()
{
  delete houghSpace;
}

Int_t HoughTransNCPi0::GetPeakHeight()
{
  Int_t binmax = static_cast<Int_t>(houghSpace->GetBinContent(houghSpace->GetMaximumBin()));
  return binmax;
}

Int_t HoughTransNCPi0::GetPeakGradBin()
{
  Int_t x=0;
  Int_t y=0;
  Int_t z=0;
  houghSpace->GetMaximumBin(x,y,z);
  return x;
}

Int_t HoughTransNCPi0::GetPeakIceptBin()
{
  Int_t x=0;
  Int_t y=0;
  Int_t z=0;
  houghSpace->GetMaximumBin(x,y,z);
  return y;
}

Float_t HoughTransNCPi0::GetPeakGradVal()
{
  Int_t x=0;
  Int_t y=0;
  Int_t z=0;
  houghSpace->GetMaximumBin(x,y,z);
  return (houghSpace->GetXaxis())->GetBinCenter(x);
}

Float_t HoughTransNCPi0::GetPeakIceptVal()
{
  Int_t x=0;
  Int_t y=0;
  Int_t z=0;
  houghSpace->GetMaximumBin(x,y,z);
  return (houghSpace->GetYaxis())->GetBinCenter(y);
}

Float_t HoughTransNCPi0::GetCos() 
{
  Int_t binx, biny;
  Float_t x = 0.;
  Float_t y = 0;
  Float_t peak = 0.;
  
  peak = houghSpace->GetBinContent(houghSpace->GetMaximumBin());
  
  for (biny=houghSpace->GetYaxis()->GetFirst();biny<=houghSpace->GetYaxis()->GetLast();++biny) {
    for (binx=houghSpace->GetXaxis()->GetFirst();binx<=houghSpace->GetXaxis()->GetLast();++binx) {
      if (houghSpace->GetBinContent(binx,biny) == peak) {
        x = houghSpace->GetXaxis()->GetBinCenter(binx);
        y = houghSpace->GetYaxis()->GetBinCenter(biny);
        //      continue; 
      }
    }
  }
   

  return x;
}
Float_t HoughTransNCPi0::GetGradientweight()
{
  Int_t binx, biny;
  Float_t x =0.;
  Float_t y=0.;
  Float_t peak = 0.;
  
  peak = houghSpace->GetBinContent(houghSpace->GetMaximumBin());
  
  for (biny=houghSpace->GetYaxis()->GetFirst();biny<=houghSpace->GetYaxis()->GetLast();++biny) {
    for (binx=houghSpace->GetXaxis()->GetFirst();binx<=houghSpace->GetXaxis()->GetLast();++binx) {
      if (houghSpace->GetBinContent(binx,biny) == peak) {
        x = houghSpace->GetXaxis()->GetBinCenter(binx);
        y = houghSpace->GetYaxis()->GetBinCenter(biny);
        //      continue; 
      }
    }
  }
  
  
  return x;
}

Float_t HoughTransNCPi0::GetTransverseVtxCoord() 
{
  Int_t binx, biny;
  Float_t x =0.;
  Float_t y =0.;
  Float_t peak = 0.;
  
  peak = houghSpace->GetBinContent(houghSpace->GetMaximumBin());
  
  for (biny=houghSpace->GetYaxis()->GetFirst();biny<=houghSpace->GetYaxis()->GetLast();++biny) {
     for (binx=houghSpace->GetXaxis()->GetFirst();binx<=houghSpace->GetXaxis()->GetLast();++binx) {
       if (houghSpace->GetBinContent(binx,biny) == peak) {
         x = houghSpace->GetXaxis()->GetBinCenter(binx);
         y = houghSpace->GetYaxis()->GetBinCenter(biny);
         //     continue; 
       }
     }
  }
  
  
  return y;
}

Float_t HoughTransNCPi0::GetIntercept() 
{
  
  Float_t intercept = GetTransverseVtxCoord() - sqrt(pow((1/GetCos()),2)-1)*GetVtxz();
  
  return intercept;
}

Float_t HoughTransNCPi0::GetRms(Float_t frac)
{
  Int_t bincount = 0;
  Int_t mRmsInput[6400];
  Int_t cRmsInput[6400];
  Float_t peakheight = GetPeakHeight();
  Float_t cos = GetCos();
  Float_t interc = GetIntercept();
  Float_t mpeak = 0;
  Float_t cpeak = 0;
  Float_t contents = 0;
  
  for(Int_t mBin=1;mBin<houghSpace->GetNbinsX();++mBin){
    for(Int_t cBin=1;cBin<houghSpace->GetNbinsY();++cBin){
      if(houghSpace->GetBinContent(mBin,cBin)>=frac*peakheight){
        mpeak = (houghSpace->GetXaxis())->GetBinCenter(mBin);
        contents = houghSpace->GetBinContent(mBin,cBin);
        cpeak = (houghSpace->GetYaxis())->GetBinCenter(cBin);
        //cout<< "mpeak " << mpeak << endl;
        //cout<< "contents = " << contents << endl;
        //cout<< "cpeak " << cpeak << endl;
        //cout << "m_av " << grad << endl;
        //cout << "c_av " << interc << endl;
        
        mRmsInput[bincount]=static_cast<Int_t>(mpeak);
        cRmsInput[bincount]=static_cast<Int_t>(cpeak); 
        bincount++;
      }
    }
  }
  if(bincount==0) return 0;
  Float_t squaresum = 0;
  for (Int_t count=1;count<bincount;++count) {
    squaresum+=(pow((mRmsInput[count]-cos),2)+pow((cRmsInput[count]-interc),2));
    
  }
  Float_t Rms = sqrt(squaresum/bincount);
  return Rms;
}

---