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