[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTest.cxx

// @(#) $Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ALICE HLT Group

#include "AliL3StandardIncludes.h"
#include "AliL3HoughTest.h"
#include "AliL3ModelTrack.h"
#include "AliL3Transform.h"
#include "AliL3Histogram.h"
#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"

#include <TRandom.h>
#include <TMath.h>
#include <TH2.h>
#include <TH3.h>
#include <TAxis.h>

#if __GNUC__ == 3
using namespace std;
#endif

//_____________________________________________________________
// AliL3HoughTest


ClassImp(AliL3HoughTest)

AliL3HoughTest::AliL3HoughTest()
{
  fData=0;
}


AliL3HoughTest::~AliL3HoughTest()
{
  if(fData)
    delete [] fData;
}

Bool_t AliL3HoughTest::GenerateTrackData(Double_t pt,Double_t psi,Double_t tgl,Int_t sign,Int_t patch,Int_t minhits)
{
  fCurrentPatch=patch;
  if(fData)
    delete fData;
  fData = new SimData[AliL3Transform::GetNRows(patch)];
  memset(fData,0,AliL3Transform::GetNRows(patch)*sizeof(SimData));
  
  AliL3ModelTrack *track = new AliL3ModelTrack();
  track->Init(0,patch);
  track->SetPt(pt);
  track->SetPsi(psi);
  track->SetTgl(tgl);
  track->SetCharge(sign);
  track->SetFirstPoint(0,0,0);
  track->CalculateHelix();

  Int_t temp[200];
  Int_t temp2[AliL3Transform::GetNTimeBins()];
  Int_t entries=100;
  Int_t clustercharge=100;
  Int_t hitcounter=0;
  for(Int_t i=AliL3Transform::GetFirstRow(patch); i<=AliL3Transform::GetLastRow(patch); i++)
    {
      
      Float_t xyz[3];
      if(!track->GetCrossingPoint(i,xyz))
	continue;

      Int_t rowindex = i - AliL3Transform::GetFirstRow(patch);
      Int_t sector,row;
      AliL3Transform::Slice2Sector(0,i,sector,row);
      AliL3Transform::Local2Raw(xyz,sector,row);

      if(xyz[1] < 0 || xyz[1] >= AliL3Transform::GetNPads(i) || xyz[2] < 0 || xyz[2] >= AliL3Transform::GetNTimeBins())
	continue;
      hitcounter++;
      track->SetPadHit(i,xyz[1]);
      track->SetTimeHit(i,xyz[2]);

      Double_t beta = track->GetCrossingAngle(i);
      track->SetCrossingAngleLUT(i,beta);
      track->CalculateClusterWidths(i);
      
      memset(temp,0,200*sizeof(Int_t));
      memset(temp2,0,AliL3Transform::GetNTimeBins()*sizeof(Int_t));
      Double_t xysigma = sqrt(track->GetParSigmaY2(i));
      Double_t zsigma = sqrt(track->GetParSigmaZ2(i));
      
      Int_t minpad=200,j;
      Int_t mintime = 1000;
      for(j=0; j<entries; j++)
	{
	  Int_t pad = TMath::Nint(gRandom->Gaus(xyz[1],xysigma));
	  Int_t time = TMath::Nint(gRandom->Gaus(xyz[2],zsigma));
	  if(pad < 0 || pad >= AliL3Transform::GetNPads(i) || time < 0 || time >= AliL3Transform::GetNTimeBins())
	    continue;
	  temp[pad]++;
	  temp2[time]++;
	  if(pad < minpad)
	    minpad=pad;
	  if(time < mintime)
	    mintime=time;
	}
      Int_t npads=0;
      for(j=0; j<200; j++)
	{
	  if(temp[j]==0) continue;
	  
	  Int_t index = j - minpad;
	  
	  if(index < 0 || index >= 10)
	    {
	      cerr<<"AliL3HoughTest::GenerateTrackData : Wrong index "<<index<<endl;
	      exit(5);
	    }
	  npads++;
	  Int_t seq_charge = clustercharge*temp[j]/entries;
	  Int_t ntimes=0;
	  for(Int_t k=0; k<AliL3Transform::GetNTimeBins(); k++)
	    {
	      if(temp2[k]==0) continue;
	      Int_t tindex = k - mintime;
	      if(tindex < 0 || tindex >= 10)
		{
		  cerr<<"AliL3HoughTest::GenerateTrackData : Wrong timeindex "<<tindex<<" "<<k<<" "<<mintime<<endl;
		  exit(5);
		}
	      Int_t charge = seq_charge*temp2[k]/entries;
	      if(charge < 3 ) 
		continue;
	      if(charge > 1023)
		charge=1023;
	      //cout<<"row "<<i<<" pad "<<j<<" time "<<k<<" charge "<<charge<<endl;
	      ntimes++;
	      fData[rowindex].pads[index][tindex]=charge;
	    }
	}
      fData[rowindex].minpad=minpad;
      fData[rowindex].mintime=mintime;
      fData[rowindex].npads=npads;
    }
  delete track;
  if(hitcounter < minhits)
    return kFALSE;
  return kTRUE;
}

void AliL3HoughTest::Transform2Circle(AliL3Histogram *hist)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::Transform : No data"<<endl;
      return;
    }
  Float_t R,phi,phi0,kappa,xyz[3];
  Int_t pad,time,charge,sector,row;
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      AliL3Transform::Slice2Sector(0,i,sector,row);
      for(Int_t j=0; j<fData[rowindex].npads; j++)
	{
	  pad = j + fData[rowindex].minpad;
	  for(Int_t k=0; k<10; k++)
	    {
	      time = k + fData[rowindex].mintime;
	      charge = fData[rowindex].pads[j][k];
	      if(charge == 0) continue;
	      AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	      
	      R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
	      
	      phi = AliL3Transform::GetPhi(xyz);
	      
	      for(Int_t k=hist->GetFirstYbin(); k<=hist->GetLastYbin(); k++)
		{
		  phi0 = hist->GetBinCenterY(k);
		  kappa = 2*sin(phi-phi0)/R;
		  hist->Fill(kappa,phi0,charge);
		}
	    }
	}
    }
}

void AliL3HoughTest::Transform2CircleC(AliL3Histogram *hist)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
      return;
    }
  Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
  Float_t r1,r2,phi1,phi2,phi_0,kappa,hit[3],hit2[3];
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	{
	  pad1 = d1 + fData[rowindex1].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time1 = j + fData[rowindex1].mintime;
	      charge1 = fData[rowindex1].pads[d1][j];
	      if(charge1==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	      r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
	      phi1 = atan2(hit[1],hit[0]);
	      
	      for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
		{
		  Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
		  for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
		    {
		      pad2 = d2 + fData[rowindex2].minpad;
		      for(Int_t k=0; k<10; k++)
			{
			  time2 = k + fData[rowindex2].mintime;
			  charge2 = fData[rowindex2].pads[d2][k];
			  if(charge2==0) continue;
			  AliL3Transform::Slice2Sector(0,j,sector,row);
			  AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
			  r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
			  phi2 = atan2(hit2[1],hit2[0]);
			  phi_0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
			  
			  kappa = 2*sin(phi1-phi_0) / r1;
			  hist->Fill(kappa,phi_0,charge1+charge2);
			}
		    }
		}
	      return;
	    }
	}
    }
}

void AliL3HoughTest::Transform2CircleF(AliL3Histogram *hist)
{
  //Fix one point in the middle of the tpc
  
  if(!fData)
    {
      cerr<<"AliL3HoughTest::TransformF : No data"<<endl;
      return;
    }
  Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
  Float_t r1,r2,phi1,phi2,phi_0,kappa,hit[3],hit2[3];
  
  Int_t rowindex1 = 80 - AliL3Transform::GetFirstRow(fCurrentPatch);
  for(Int_t d1=1; d1<fData[rowindex1].npads; d1++)
    {
      pad1 = d1 + fData[rowindex1].minpad;
      
      AliL3Transform::Slice2Sector(0,80,sector,row);
      AliL3Transform::Raw2Local(hit,sector,row,pad1,0);
      r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
      phi1 = atan2(hit[1],hit[0]);
      
      for(Int_t j=0; j<10; j++)
	{
	  time1 = j + fData[rowindex1].mintime;
	  charge1 = fData[rowindex1].pads[d1][j];
	  if(charge1==0) continue;
	  
	  for(Int_t j=0; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
	    {
	      if(j==80) continue;
	      Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
	      for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
		{
		  pad2 = d2 + fData[rowindex2].minpad;
		  for(Int_t k=0; k<10; k++)
		    {
		      time2 = k + fData[rowindex2].mintime;
		      charge2 = fData[rowindex2].pads[d2][k];
		      if(charge2==0) continue;
		      AliL3Transform::Slice2Sector(0,j,sector,row);
		      AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
		      r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
		      phi2 = atan2(hit2[1],hit2[0]);
		      phi_0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
		      
		      kappa = 2*sin(phi1-phi_0) / r1;
		      hist->Fill(kappa,phi_0,charge1+charge2);
		      //cout<<"Filling "<<kappa<<" psi "<<phi_0<<" charge "<<charge1<<endl;
		    }
		}
	    }
	}
      return;
    }
}

void AliL3HoughTest::Transform2Line(AliL3Histogram *hist,Int_t *rowrange)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
      return;
    }
  
  Int_t pad,time,charge,sector,row;
  Float_t hit[3],theta,rho;
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      if(i<rowrange[0])
	continue;
      if(i>rowrange[1])
	break;
      Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d=0; d<fData[rowindex].npads; d++)
	{
	  pad = d + fData[rowindex].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time = j + fData[rowindex].mintime;
	      charge = fData[rowindex].pads[d][j];
	      if(charge==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad,time);
	      
	      hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	      
	      for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
		{
		  theta = hist->GetBinCenterX(xbin);
		  rho = hit[0]*cos(theta) + hit[1]*sin(theta);
		  hist->Fill(theta,rho,charge);
		}
	    }
	}
    }
}

void AliL3HoughTest::Transform2LineC(AliL3Histogram *hist,Int_t *rowrange)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
      return;
    }
  
  Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
  Float_t theta,rho,hit[3],hit2[3];
  for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
    {
      Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	{
	  pad1 = d1 + fData[rowindex1].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time1 = j + fData[rowindex1].mintime;
	      charge1 = fData[rowindex1].pads[d1][j];
	      if(charge1==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	      
	      hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	      
	      for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
		{
		  Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
		  for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
		    {
		      pad2 = d2 + fData[rowindex2].minpad;
		      for(Int_t k=0; k<10; k++)
			{
			  time2 = k + fData[rowindex2].mintime;
			  charge2 = fData[rowindex2].pads[d2][k];
			  if(charge2==0) continue;
			  AliL3Transform::Slice2Sector(0,i2,sector,row);
			  AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
			  
			  hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
			  
			  theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
			  rho = hit[0]*cos(theta)+hit[1]*sin(theta);
			  hist->Fill(theta,rho,charge1+charge2);
			}
		    }
		}
	    }
	}
    }
}

void AliL3HoughTest::FillImage(TH2 *hist,Int_t displayrow)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::FillImage : No data to fill"<<endl;
      return;
    }
  
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      if(displayrow >=0)
	if(i != displayrow) continue;

      //cout<<"row "<<i<<" npads "<<fData[rowindex].npads<<endl;
      for(Int_t j=0; j<fData[rowindex].npads; j++)
	{
	  Int_t pad = j + fData[rowindex].minpad;
	  for(Int_t k=0; k<10; k++)
	    {
	      Int_t time = k + fData[rowindex].mintime;
	      Int_t charge = fData[rowindex].pads[j][k];
	      if(charge==0) continue;
	      //cout<<i<<" "<<pad<<" "<<time<<" "<<charge<<endl;
	      Float_t xyz[3];
	      Int_t sector,row;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	      if(displayrow>=0)
		hist->Fill(pad,time,charge);
	      else
		hist->Fill(xyz[0],xyz[1],charge);
	    }
	}
      if(displayrow>=0)
	break;
    }
}

void AliL3HoughTest::Transform2Line3D(TH3 *hist,Int_t *rowrange,Float_t *phirange)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
      return;
    }
  
  Int_t pad,time,charge,sector,row;
  Float_t hit[3],theta,rho,R,delta;
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      
      if(i<rowrange[0])
	continue;
      if(i>rowrange[1])
	break;
      
      Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d=0; d<fData[rowindex].npads; d++)
	{
	  pad = d + fData[rowindex].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time = j + fData[rowindex].mintime;
	      charge = fData[rowindex].pads[d][j];
	      if(charge==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad,time);
	      
	      Float_t phi = AliL3Transform::GetPhi(hit);
	      if(phi < phirange[0] || phi > phirange[1])
		continue;
	      
	      hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	      Float_t x = hit[0] + AliL3Transform::Row2X(rowrange[0]);
	      R = sqrt(x*x + hit[1]*hit[1]);
	      delta = atan(hit[2]/R);
	      
	      for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
		{
		  theta = hist->GetXaxis()->GetBinCenter(xbin);
		  rho = hit[0]*cos(theta) + hit[1]*sin(theta);
		  hist->Fill(theta,rho,delta,charge);
		}
	    }
	}
    }
}

void AliL3HoughTest::Transform2LineC3D(TH3 *hist,Int_t *rowrange)
{
  if(!fData)
    {
      cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
      return;
    }
  
  Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
  Float_t theta,rho,hit[3],hit2[3],R1,R2,delta,delta1,delta2;
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      if(i<rowrange[0])
	continue;
      if(i>rowrange[1])
	break;
      Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	{
	  pad1 = d1 + fData[rowindex1].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time1 = j + fData[rowindex1].mintime;
	      charge1 = fData[rowindex1].pads[d1][j];
	      if(charge1==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	      R1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
	      delta1 = atan(hit[2]/R1);
	      hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	      
	      for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
		{
		  Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
		  for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
		    {
		      pad2 = d2 + fData[rowindex2].minpad;
		      for(Int_t k=0; k<10; k++)
			{
			  time2 = k + fData[rowindex2].mintime;
			  charge2 = fData[rowindex2].pads[d2][k];
			  if(charge2==0) continue;
			  AliL3Transform::Slice2Sector(0,i2,sector,row);
			  AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
			  R2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
			  delta2 = atan(hit2[2]/R2);
			  delta = (charge1*delta1 + charge2*delta2)/(charge1+charge2);
			  hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
			  
			  theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
			  rho = hit[0]*cos(theta)+hit[1]*sin(theta);
			  hist->Fill(theta,rho,delta,charge1+charge2);
			}
		    }
		}
	    }
	}
    }  
}

void AliL3HoughTest::TransformLines2Circle(TH3 *hist,AliL3TrackArray *tracks)
{
  
  for(Int_t i=0; i<tracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!tr) continue;
      Int_t middlerow = (Int_t)(tr->GetFirstRow()+(tr->GetLastRow()-tr->GetFirstRow())/2);
      Float_t hit[3];
      tr->GetLineCrossingPoint(middlerow,hit);
      hit[0] += AliL3Transform::Row2X(tr->GetFirstRow());
      Float_t R = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
      hit[2] = R*tr->GetTgl();
      Float_t phi = atan2(hit[1],hit[0]);
      Float_t theta = tr->GetPsiLine() - AliL3Transform::Pi()/2;
      Float_t psi = 2*phi - theta;
      Float_t kappa = 2/R*sin(phi-psi);
      Float_t delta = atan(hit[2]/R);
      hist->Fill(kappa,psi,delta,tr->GetWeight());
      
    }
}

void AliL3HoughTest::Transform2Center(AliL3Histogram *hist)
{
  //Choose parameter space to be center of curvature.
  
  if(!fData)
    {
      cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
      return;
    }
  Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
  Float_t r1,r2,phi1,phi2,hit[3],hit2[3];
  //Float_t phi_0,kappa;
  for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
    {
      Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
      for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	{
	  pad1 = d1 + fData[rowindex1].minpad;
	  for(Int_t j=0; j<10; j++)
	    {
	      time1 = j + fData[rowindex1].mintime;
	      charge1 = fData[rowindex1].pads[d1][j];
	      if(charge1==0) continue;
	      AliL3Transform::Slice2Sector(0,i,sector,row);
	      AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	      r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1])/2;
	      phi1 = atan2(hit[1],hit[0]);
	      
	      for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
		{
		  Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
		  for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
		    {
		      pad2 = d2 + fData[rowindex2].minpad;
		      for(Int_t k=0; k<10; k++)
			{
			  time2 = k + fData[rowindex2].mintime;
			  charge2 = fData[rowindex2].pads[d2][k];
			  if(charge2==0) continue;
			  AliL3Transform::Slice2Sector(0,j,sector,row);
			  AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
			  r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1])/2;
			  phi2 = atan2(hit2[1],hit2[0]);
			  Float_t yc = (r2-(r1/cos(phi1))*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
			  Float_t xc = r1/cos(phi1) - (r2*tan(phi1)-r1*sin(phi1)*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
			  hist->Fill(xc,yc,charge1+charge2);
			}
		    }
		}
	    }
	}
    }
}


void AliL3HoughTest::FindAbsMaxima(TH3 *hist,Int_t zsearch,Float_t &max_x,Float_t &max_y,Float_t &max_z,Int_t &maxvalue)
{
  
  TH1 *h1 = hist->Project3D("z");
      
  TAxis *z = hist->GetZaxis();
  Float_t zpeak[50];
  Int_t n=0,i,zbin[50];
  for(i=z->GetFirst()+1; i<z->GetLast()-1; i++)
    {
      int bin1 = h1->GetBin(i-1);
      int bin2 = h1->GetBin(i);
      int bin3 = h1->GetBin(i+1);
      if(h1->GetBinContent(bin2) > h1->GetBinContent(bin1) && h1->GetBinContent(bin2) > h1->GetBinContent(bin3))
	{
	  zbin[n]=bin2;
	  zpeak[n++]=h1->GetBinCenter(bin2);
	}
    }
  
  Int_t zrange[2] = {z->GetFirst(),z->GetLast()};
  z->SetRange(zbin[0]-zsearch,zbin[0]+zsearch);

  TH1 *h2 = hist->Project3D("yx");
  z->SetRange(zrange[0],zrange[1]);

  TAxis *x = h2->GetXaxis();
  TAxis *y = h2->GetYaxis();
  maxvalue=0;
  for(i=0; i<n; i++)
    {
      for(Int_t xbin=x->GetFirst(); xbin<=x->GetLast(); xbin++)
	{
	  Float_t xvalue = x->GetBinCenter(xbin);
	  for(Int_t ybin=y->GetFirst(); ybin<=y->GetLast(); ybin++)
	    {
	      Float_t yvalue = y->GetBinCenter(ybin);
	      Int_t value = (Int_t)h2->GetBinContent(xbin,ybin);
	      if(value > maxvalue)
		{
		  maxvalue=value;
		  max_x = xvalue;
		  max_y = yvalue;
		  max_z = zpeak[i];
		}
	    }
	}
    }
  
  delete h1;
  delete h2;
}
Commit	Line	Data
3e87ef69	1	// @(#) $Id$
7a1f1432	2
7a1f1432	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
3e87ef69	4	//*-- Copyright &copy ALICE HLT Group
7a1f1432	5
	6	#include "AliL3StandardIncludes.h"
	7	#include "AliL3HoughTest.h"
	8	#include "AliL3ModelTrack.h"
	9	#include "AliL3Transform.h"
	10	#include "AliL3Histogram.h"
3e87ef69	11	#include "AliL3TrackArray.h"
	12	#include "AliL3HoughTrack.h"
	13
	14	#include <TRandom.h>
	15	#include <TMath.h>
	16	#include <TH2.h>
	17	#include <TH3.h>
	18	#include <TAxis.h>
7a1f1432	19
0bd0c1ef	20	#if __GNUC__ == 3
7a1f1432	21	using namespace std;
	22	#endif
	23
	24	//_____________________________________________________________
	25	// AliL3HoughTest
	26
	27
	28	ClassImp(AliL3HoughTest)
	29
	30	AliL3HoughTest::AliL3HoughTest()
	31	{
	32	fData=0;
	33	}
	34
	35
	36	AliL3HoughTest::~AliL3HoughTest()
	37	{
	38	if(fData)
	39	delete [] fData;
	40	}
	41
c3853fab	42	Bool_t AliL3HoughTest::GenerateTrackData(Double_t pt,Double_t psi,Double_t tgl,Int_t sign,Int_t patch,Int_t minhits)
7a1f1432	43	{
	44	fCurrentPatch=patch;
	45	if(fData)
	46	delete fData;
	47	fData = new SimData[AliL3Transform::GetNRows(patch)];
	48	memset(fData,0,AliL3Transform::GetNRows(patch)*sizeof(SimData));
	49
	50	AliL3ModelTrack *track = new AliL3ModelTrack();
	51	track->Init(0,patch);
	52	track->SetPt(pt);
	53	track->SetPsi(psi);
c3853fab	54	track->SetTgl(tgl);
d7ab65da	55	track->SetCharge(sign);
d7ab65da	56	track->SetFirstPoint(0,0,0);
7a1f1432	57	track->CalculateHelix();
7a1f1432	58
7a1f1432	59	Int_t temp[200];
c3853fab	60	Int_t temp2[AliL3Transform::GetNTimeBins()];
7a1f1432	61	Int_t entries=100;
7a1f1432	62	Int_t clustercharge=100;
d7ab65da	63	Int_t hitcounter=0;
7a1f1432	64	for(Int_t i=AliL3Transform::GetFirstRow(patch); i<=AliL3Transform::GetLastRow(patch); i++)
7a1f1432	65	{
7a1f1432	66
3e87ef69	67	Float_t xyz[3];
7a1f1432	68	if(!track->GetCrossingPoint(i,xyz))
7a1f1432	69	continue;
3e87ef69	70
d7ab65da	71	Int_t rowindex = i - AliL3Transform::GetFirstRow(patch);
7a1f1432	72	Int_t sector,row;
	73	AliL3Transform::Slice2Sector(0,i,sector,row);
	74	AliL3Transform::Local2Raw(xyz,sector,row);
	75
c3853fab	76	if(xyz[1] < 0 \|\| xyz[1] >= AliL3Transform::GetNPads(i) \|\| xyz[2] < 0 \|\| xyz[2] >= AliL3Transform::GetNTimeBins())
7a1f1432	77	continue;
d7ab65da	78	hitcounter++;
d7ab65da	79	track->SetPadHit(i,xyz[1]);
c3853fab	80	track->SetTimeHit(i,xyz[2]);
3e87ef69	81
	82	Double_t beta = track->GetCrossingAngle(i);
	83	track->SetCrossingAngleLUT(i,beta);
	84	track->CalculateClusterWidths(i);
c3853fab	85
7a1f1432	86	memset(temp,0,200*sizeof(Int_t));
c3853fab	87	memset(temp2,0,AliL3Transform::GetNTimeBins()*sizeof(Int_t));
7a1f1432	88	Double_t xysigma = sqrt(track->GetParSigmaY2(i));
c3853fab	89	Double_t zsigma = sqrt(track->GetParSigmaZ2(i));
3e87ef69	90
7a1f1432	91	Int_t minpad=200,j;
c3853fab	92	Int_t mintime = 1000;
7a1f1432	93	for(j=0; j<entries; j++)
	94	{
	95	Int_t pad = TMath::Nint(gRandom->Gaus(xyz[1],xysigma));
c3853fab	96	Int_t time = TMath::Nint(gRandom->Gaus(xyz[2],zsigma));
c3853fab	97	if(pad < 0 \|\| pad >= AliL3Transform::GetNPads(i) \|\| time < 0 \|\| time >= AliL3Transform::GetNTimeBins())
7a1f1432	98	continue;
7a1f1432	99	temp[pad]++;
c3853fab	100	temp2[time]++;
7a1f1432	101	if(pad < minpad)
7a1f1432	102	minpad=pad;
c3853fab	103	if(time < mintime)
c3853fab	104	mintime=time;
7a1f1432	105	}
c3853fab	106	Int_t npads=0;
7a1f1432	107	for(j=0; j<200; j++)
	108	{
	109	if(temp[j]==0) continue;
c3853fab	110
7a1f1432	111	Int_t index = j - minpad;
c3853fab	112
7a1f1432	113	if(index < 0 \|\| index >= 10)
	114	{
	115	cerr<<"AliL3HoughTest::GenerateTrackData : Wrong index "<<index<<endl;
	116	exit(5);
	117	}
c3853fab	118	npads++;
	119	Int_t seq_charge = clustercharge*temp[j]/entries;
	120	Int_t ntimes=0;
	121	for(Int_t k=0; k<AliL3Transform::GetNTimeBins(); k++)
	122	{
	123	if(temp2[k]==0) continue;
	124	Int_t tindex = k - mintime;
	125	if(tindex < 0 \|\| tindex >= 10)
	126	{
	127	cerr<<"AliL3HoughTest::GenerateTrackData : Wrong timeindex "<<tindex<<" "<<k<<" "<<mintime<<endl;
	128	exit(5);
	129	}
	130	Int_t charge = seq_charge*temp2[k]/entries;
	131	if(charge < 3 )
	132	continue;
	133	if(charge > 1023)
	134	charge=1023;
	135	//cout<<"row "<<i<<" pad "<<j<<" time "<<k<<" charge "<<charge<<endl;
	136	ntimes++;
	137	fData[rowindex].pads[index][tindex]=charge;
	138	}
7a1f1432	139	}
d7ab65da	140	fData[rowindex].minpad=minpad;
c3853fab	141	fData[rowindex].mintime=mintime;
c3853fab	142	fData[rowindex].npads=npads;
7a1f1432	143	}
7a1f1432	144	delete track;
d7ab65da	145	if(hitcounter < minhits)
	146	return kFALSE;
	147	return kTRUE;
7a1f1432	148	}
7a1f1432	149
d7ab65da	150	void AliL3HoughTest::Transform2Circle(AliL3Histogram *hist)
7a1f1432	151	{
	152	if(!fData)
	153	{
	154	cerr<<"AliL3HoughTest::Transform : No data"<<endl;
	155	return;
	156	}
	157	Float_t R,phi,phi0,kappa,xyz[3];
c3853fab	158	Int_t pad,time,charge,sector,row;
7a1f1432	159	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
7a1f1432	160	{
d7ab65da	161	Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
7a1f1432	162	AliL3Transform::Slice2Sector(0,i,sector,row);
d7ab65da	163	for(Int_t j=0; j<fData[rowindex].npads; j++)
7a1f1432	164	{
d7ab65da	165	pad = j + fData[rowindex].minpad;
c3853fab	166	for(Int_t k=0; k<10; k++)
7a1f1432	167	{
c3853fab	168	time = k + fData[rowindex].mintime;
	169	charge = fData[rowindex].pads[j][k];
	170	if(charge == 0) continue;
	171	AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	172
	173	R = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	174
	175	phi = AliL3Transform::GetPhi(xyz);
	176
	177	for(Int_t k=hist->GetFirstYbin(); k<=hist->GetLastYbin(); k++)
	178	{
	179	phi0 = hist->GetBinCenterY(k);
	180	kappa = 2*sin(phi-phi0)/R;
	181	hist->Fill(kappa,phi0,charge);
	182	}
7a1f1432	183	}
	184	}
	185	}
	186	}
	187
d7ab65da	188	void AliL3HoughTest::Transform2CircleC(AliL3Histogram *hist)
7a1f1432	189	{
	190	if(!fData)
	191	{
	192	cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
	193	return;
	194	}
c3853fab	195	Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
7a1f1432	196	Float_t r1,r2,phi1,phi2,phi_0,kappa,hit[3],hit2[3];
	197	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
	198	{
d7ab65da	199	Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
d7ab65da	200	for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
7a1f1432	201	{
d7ab65da	202	pad1 = d1 + fData[rowindex1].minpad;
c3853fab	203	for(Int_t j=0; j<10; j++)
7a1f1432	204	{
c3853fab	205	time1 = j + fData[rowindex1].mintime;
	206	charge1 = fData[rowindex1].pads[d1][j];
	207	if(charge1==0) continue;
	208	AliL3Transform::Slice2Sector(0,i,sector,row);
	209	AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	210	r1 = sqrt(hit[0]hit[0]+hit[1]hit[1]);
	211	phi1 = atan2(hit[1],hit[0]);
	212
	213	for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
7a1f1432	214	{
c3853fab	215	Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
	216	for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
	217	{
	218	pad2 = d2 + fData[rowindex2].minpad;
	219	for(Int_t k=0; k<10; k++)
	220	{
	221	time2 = k + fData[rowindex2].mintime;
	222	charge2 = fData[rowindex2].pads[d2][k];
	223	if(charge2==0) continue;
	224	AliL3Transform::Slice2Sector(0,j,sector,row);
	225	AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
	226	r2 = sqrt(hit2[0]hit2[0]+hit2[1]hit2[1]);
	227	phi2 = atan2(hit2[1],hit2[0]);
	228	phi_0 = atan( (r2sin(phi1) - r1sin(phi2)) / (r2cos(phi1) - r1cos(phi2)) );
	229
	230	kappa = 2*sin(phi1-phi_0) / r1;
	231	hist->Fill(kappa,phi_0,charge1+charge2);
c3853fab	232	}
c3853fab	233	}
7a1f1432	234	}
3e87ef69	235	return;
7a1f1432	236	}
	237	}
	238	}
	239	}
	240
3e87ef69	241	void AliL3HoughTest::Transform2CircleF(AliL3Histogram *hist)
	242	{
	243	//Fix one point in the middle of the tpc
	244
	245	if(!fData)
	246	{
	247	cerr<<"AliL3HoughTest::TransformF : No data"<<endl;
	248	return;
	249	}
	250	Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
	251	Float_t r1,r2,phi1,phi2,phi_0,kappa,hit[3],hit2[3];
	252
	253	Int_t rowindex1 = 80 - AliL3Transform::GetFirstRow(fCurrentPatch);
	254	for(Int_t d1=1; d1<fData[rowindex1].npads; d1++)
	255	{
	256	pad1 = d1 + fData[rowindex1].minpad;
	257
	258	AliL3Transform::Slice2Sector(0,80,sector,row);
	259	AliL3Transform::Raw2Local(hit,sector,row,pad1,0);
	260	r1 = sqrt(hit[0]hit[0]+hit[1]hit[1]);
	261	phi1 = atan2(hit[1],hit[0]);
	262
	263	for(Int_t j=0; j<10; j++)
	264	{
	265	time1 = j + fData[rowindex1].mintime;
	266	charge1 = fData[rowindex1].pads[d1][j];
	267	if(charge1==0) continue;
	268
	269	for(Int_t j=0; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
	270	{
	271	if(j==80) continue;
	272	Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
	273	for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
	274	{
	275	pad2 = d2 + fData[rowindex2].minpad;
	276	for(Int_t k=0; k<10; k++)
	277	{
	278	time2 = k + fData[rowindex2].mintime;
	279	charge2 = fData[rowindex2].pads[d2][k];
	280	if(charge2==0) continue;
	281	AliL3Transform::Slice2Sector(0,j,sector,row);
	282	AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
	283	r2 = sqrt(hit2[0]hit2[0]+hit2[1]hit2[1]);
	284	phi2 = atan2(hit2[1],hit2[0]);
	285	phi_0 = atan( (r2sin(phi1) - r1sin(phi2)) / (r2cos(phi1) - r1cos(phi2)) );
	286
	287	kappa = 2*sin(phi1-phi_0) / r1;
	288	hist->Fill(kappa,phi_0,charge1+charge2);
	289	//cout<<"Filling "<<kappa<<" psi "<<phi_0<<" charge "<<charge1<<endl;
	290	}
	291	}
	292	}
	293	}
	294	return;
	295	}
	296	}
	297
d7ab65da	298	void AliL3HoughTest::Transform2Line(AliL3Histogram hist,Int_t rowrange)
	299	{
	300	if(!fData)
	301	{
	302	cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
	303	return;
	304	}
	305
c3853fab	306	Int_t pad,time,charge,sector,row;
d7ab65da	307	Float_t hit[3],theta,rho;
3e87ef69	308	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
d7ab65da	309	{
3e87ef69	310	if(i<rowrange[0])
	311	continue;
	312	if(i>rowrange[1])
	313	break;
d7ab65da	314	Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
	315	for(Int_t d=0; d<fData[rowindex].npads; d++)
	316	{
	317	pad = d + fData[rowindex].minpad;
c3853fab	318	for(Int_t j=0; j<10; j++)
d7ab65da	319	{
c3853fab	320	time = j + fData[rowindex].mintime;
	321	charge = fData[rowindex].pads[d][j];
	322	if(charge==0) continue;
	323	AliL3Transform::Slice2Sector(0,i,sector,row);
	324	AliL3Transform::Raw2Local(hit,sector,row,pad,time);
	325
	326	hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	327
	328	for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
	329	{
	330	theta = hist->GetBinCenterX(xbin);
	331	rho = hit[0]cos(theta) + hit[1]sin(theta);
	332	hist->Fill(theta,rho,charge);
	333	}
d7ab65da	334	}
	335	}
	336	}
	337	}
	338
	339	void AliL3HoughTest::Transform2LineC(AliL3Histogram hist,Int_t rowrange)
	340	{
	341	if(!fData)
	342	{
	343	cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
	344	return;
	345	}
	346
c3853fab	347	Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
d7ab65da	348	Float_t theta,rho,hit[3],hit2[3];
	349	for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
	350	{
	351	Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
	352	for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	353	{
	354	pad1 = d1 + fData[rowindex1].minpad;
c3853fab	355	for(Int_t j=0; j<10; j++)
d7ab65da	356	{
c3853fab	357	time1 = j + fData[rowindex1].mintime;
	358	charge1 = fData[rowindex1].pads[d1][j];
	359	if(charge1==0) continue;
	360	AliL3Transform::Slice2Sector(0,i,sector,row);
	361	AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	362
	363	hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	364
3e87ef69	365	for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
d7ab65da	366	{
3e87ef69	367	Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
c3853fab	368	for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
	369	{
	370	pad2 = d2 + fData[rowindex2].minpad;
	371	for(Int_t k=0; k<10; k++)
	372	{
	373	time2 = k + fData[rowindex2].mintime;
	374	charge2 = fData[rowindex2].pads[d2][k];
	375	if(charge2==0) continue;
3e87ef69	376	AliL3Transform::Slice2Sector(0,i2,sector,row);
c3853fab	377	AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
	378
	379	hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
	380
	381	theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
	382	rho = hit[0]cos(theta)+hit[1]sin(theta);
	383	hist->Fill(theta,rho,charge1+charge2);
	384	}
	385	}
d7ab65da	386	}
	387	}
	388	}
	389	}
	390	}
	391
3e87ef69	392	void AliL3HoughTest::FillImage(TH2 *hist,Int_t displayrow)
7a1f1432	393	{
	394	if(!fData)
	395	{
	396	cerr<<"AliL3HoughTest::FillImage : No data to fill"<<endl;
	397	return;
	398	}
c3853fab	399
7a1f1432	400	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
7a1f1432	401	{
d7ab65da	402	Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
3e87ef69	403	if(displayrow >=0)
3e87ef69	404	if(i != displayrow) continue;
c3853fab	405
c3853fab	406	//cout<<"row "<<i<<" npads "<<fData[rowindex].npads<<endl;
d7ab65da	407	for(Int_t j=0; j<fData[rowindex].npads; j++)
7a1f1432	408	{
d7ab65da	409	Int_t pad = j + fData[rowindex].minpad;
c3853fab	410	for(Int_t k=0; k<10; k++)
	411	{
	412	Int_t time = k + fData[rowindex].mintime;
	413	Int_t charge = fData[rowindex].pads[j][k];
	414	if(charge==0) continue;
	415	//cout<<i<<" "<<pad<<" "<<time<<" "<<charge<<endl;
	416	Float_t xyz[3];
	417	Int_t sector,row;
	418	AliL3Transform::Slice2Sector(0,i,sector,row);
	419	AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
3e87ef69	420	if(displayrow>=0)
c3853fab	421	hist->Fill(pad,time,charge);
	422	else
	423	hist->Fill(xyz[0],xyz[1],charge);
	424	}
7a1f1432	425	}
3e87ef69	426	if(displayrow>=0)
c3853fab	427	break;
7a1f1432	428	}
7a1f1432	429	}
c3853fab	430
3e87ef69	431	void AliL3HoughTest::Transform2Line3D(TH3 hist,Int_t rowrange,Float_t *phirange)
c3853fab	432	{
	433	if(!fData)
	434	{
	435	cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
	436	return;
	437	}
	438
	439	Int_t pad,time,charge,sector,row;
	440	Float_t hit[3],theta,rho,R,delta;
3e87ef69	441	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
c3853fab	442	{
3e87ef69	443
	444	if(i<rowrange[0])
	445	continue;
	446	if(i>rowrange[1])
	447	break;
	448
c3853fab	449	Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
	450	for(Int_t d=0; d<fData[rowindex].npads; d++)
	451	{
	452	pad = d + fData[rowindex].minpad;
	453	for(Int_t j=0; j<10; j++)
	454	{
	455	time = j + fData[rowindex].mintime;
	456	charge = fData[rowindex].pads[d][j];
	457	if(charge==0) continue;
	458	AliL3Transform::Slice2Sector(0,i,sector,row);
	459	AliL3Transform::Raw2Local(hit,sector,row,pad,time);
	460
3e87ef69	461	Float_t phi = AliL3Transform::GetPhi(hit);
	462	if(phi < phirange[0] \|\| phi > phirange[1])
	463	continue;
	464
c3853fab	465	hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
3e87ef69	466	Float_t x = hit[0] + AliL3Transform::Row2X(rowrange[0]);
	467	R = sqrt(xx + hit[1]hit[1]);
	468	delta = atan(hit[2]/R);
c3853fab	469
	470	for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
	471	{
	472	theta = hist->GetXaxis()->GetBinCenter(xbin);
	473	rho = hit[0]cos(theta) + hit[1]sin(theta);
c3853fab	474	hist->Fill(theta,rho,delta,charge);
	475	}
	476	}
	477	}
	478	}
	479	}
	480
	481	void AliL3HoughTest::Transform2LineC3D(TH3 hist,Int_t rowrange)
	482	{
	483	if(!fData)
	484	{
	485	cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
	486	return;
	487	}
	488
	489	Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
	490	Float_t theta,rho,hit[3],hit2[3],R1,R2,delta,delta1,delta2;
3e87ef69	491	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
c3853fab	492	{
3e87ef69	493	if(i<rowrange[0])
	494	continue;
	495	if(i>rowrange[1])
	496	break;
c3853fab	497	Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
	498	for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	499	{
	500	pad1 = d1 + fData[rowindex1].minpad;
	501	for(Int_t j=0; j<10; j++)
	502	{
	503	time1 = j + fData[rowindex1].mintime;
	504	charge1 = fData[rowindex1].pads[d1][j];
	505	if(charge1==0) continue;
	506	AliL3Transform::Slice2Sector(0,i,sector,row);
	507	AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	508	R1 = sqrt(hit[0]hit[0]+hit[1]hit[1]);
	509	delta1 = atan(hit[2]/R1);
	510	hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
	511
3e87ef69	512	for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
c3853fab	513	{
3e87ef69	514	Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
c3853fab	515	for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
	516	{
	517	pad2 = d2 + fData[rowindex2].minpad;
	518	for(Int_t k=0; k<10; k++)
	519	{
	520	time2 = k + fData[rowindex2].mintime;
	521	charge2 = fData[rowindex2].pads[d2][k];
	522	if(charge2==0) continue;
3e87ef69	523	AliL3Transform::Slice2Sector(0,i2,sector,row);
c3853fab	524	AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
	525	R2 = sqrt(hit2[0]hit2[0]+hit2[1]hit2[1]);
	526	delta2 = atan(hit2[2]/R2);
	527	delta = (charge1delta1 + charge2delta2)/(charge1+charge2);
	528	hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
	529
	530	theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
	531	rho = hit[0]cos(theta)+hit[1]sin(theta);
	532	hist->Fill(theta,rho,delta,charge1+charge2);
	533	}
	534	}
	535	}
	536	}
	537	}
	538	}
	539	}
3e87ef69	540
	541	void AliL3HoughTest::TransformLines2Circle(TH3 hist,AliL3TrackArray tracks)
	542	{
	543
	544	for(Int_t i=0; i<tracks->GetNTracks(); i++)
	545	{
	546	AliL3HoughTrack tr = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	547	if(!tr) continue;
	548	Int_t middlerow = (Int_t)(tr->GetFirstRow()+(tr->GetLastRow()-tr->GetFirstRow())/2);
	549	Float_t hit[3];
	550	tr->GetLineCrossingPoint(middlerow,hit);
	551	hit[0] += AliL3Transform::Row2X(tr->GetFirstRow());
	552	Float_t R = sqrt(hit[0]hit[0] + hit[1]hit[1]);
	553	hit[2] = R*tr->GetTgl();
	554	Float_t phi = atan2(hit[1],hit[0]);
	555	Float_t theta = tr->GetPsiLine() - AliL3Transform::Pi()/2;
	556	Float_t psi = 2*phi - theta;
	557	Float_t kappa = 2/R*sin(phi-psi);
	558	Float_t delta = atan(hit[2]/R);
	559	hist->Fill(kappa,psi,delta,tr->GetWeight());
	560
	561	}
	562	}
	563
	564	void AliL3HoughTest::Transform2Center(AliL3Histogram *hist)
	565	{
	566	//Choose parameter space to be center of curvature.
	567
	568	if(!fData)
	569	{
	570	cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
	571	return;
	572	}
	573	Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
	574	Float_t r1,r2,phi1,phi2,hit[3],hit2[3];
	575	//Float_t phi_0,kappa;
	576	for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
	577	{
	578	Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
	579	for(Int_t d1=0; d1<fData[rowindex1].npads; d1++)
	580	{
	581	pad1 = d1 + fData[rowindex1].minpad;
	582	for(Int_t j=0; j<10; j++)
	583	{
	584	time1 = j + fData[rowindex1].mintime;
	585	charge1 = fData[rowindex1].pads[d1][j];
	586	if(charge1==0) continue;
	587	AliL3Transform::Slice2Sector(0,i,sector,row);
	588	AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
	589	r1 = sqrt(hit[0]hit[0]+hit[1]hit[1])/2;
	590	phi1 = atan2(hit[1],hit[0]);
	591
	592	for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
	593	{
	594	Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
	595	for(Int_t d2=0; d2<fData[rowindex2].npads; d2++)
	596	{
	597	pad2 = d2 + fData[rowindex2].minpad;
	598	for(Int_t k=0; k<10; k++)
	599	{
	600	time2 = k + fData[rowindex2].mintime;
	601	charge2 = fData[rowindex2].pads[d2][k];
	602	if(charge2==0) continue;
	603	AliL3Transform::Slice2Sector(0,j,sector,row);
604	AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
605	r2 = sqrt(hit2[0]hit2[0]+hit2[1]hit2[1])/2;
606	phi2 = atan2(hit2[1],hit2[0]);
607	Float_t yc = (r2-(r1/cos(phi1))cos(phi2))/(sin(phi2)-tan(phi1)cos(phi2));
608	Float_t xc = r1/cos(phi1) - (r2tan(phi1)-r1sin(phi1)cos(phi2))/(sin(phi2)-tan(phi1)cos(phi2));
609	hist->Fill(xc,yc,charge1+charge2);
610	}
611	}
612	}
613	}
614	}
615	}
616	}
617
618
619	void AliL3HoughTest::FindAbsMaxima(TH3 *hist,Int_t zsearch,Float_t &max_x,Float_t &max_y,Float_t &max_z,Int_t &maxvalue)
620	{
621
622	TH1 *h1 = hist->Project3D("z");
623
624	TAxis *z = hist->GetZaxis();
625	Float_t zpeak[50];
626	Int_t n=0,i,zbin[50];
627	for(i=z->GetFirst()+1; i<z->GetLast()-1; i++)
628	{
629	int bin1 = h1->GetBin(i-1);
630	int bin2 = h1->GetBin(i);
631	int bin3 = h1->GetBin(i+1);
632	if(h1->GetBinContent(bin2) > h1->GetBinContent(bin1) && h1->GetBinContent(bin2) > h1->GetBinContent(bin3))
633	{
634	zbin[n]=bin2;
635	zpeak[n++]=h1->GetBinCenter(bin2);
636	}
637	}
638
639	Int_t zrange[2] = {z->GetFirst(),z->GetLast()};
640	z->SetRange(zbin[0]-zsearch,zbin[0]+zsearch);
641
642	TH1 *h2 = hist->Project3D("yx");
643	z->SetRange(zrange[0],zrange[1]);
644
645	TAxis *x = h2->GetXaxis();
646	TAxis *y = h2->GetYaxis();
647	maxvalue=0;
648	for(i=0; i<n; i++)
649	{
650	for(Int_t xbin=x->GetFirst(); xbin<=x->GetLast(); xbin++)
651	{
652	Float_t xvalue = x->GetBinCenter(xbin);
653	for(Int_t ybin=y->GetFirst(); ybin<=y->GetLast(); ybin++)
654	{
655	Float_t yvalue = y->GetBinCenter(ybin);
656	Int_t value = (Int_t)h2->GetBinContent(xbin,ybin);
657	if(value > maxvalue)
658	{
659	maxvalue=value;
660	max_x = xvalue;
661	max_y = yvalue;
662	max_z = zpeak[i];
663	}
664	}
665	}
666	}
667
668	delete h1;
669	delete h2;
670	}
671