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

// @(#) $Id$

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

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3HoughTransformerGlobal.h"
#include "AliL3FileHandler.h"
#include "AliL3Transform.h"
#include "AliL3DigitData.h"
#include "AliL3Track.h"
#include "AliL3Histogram.h"
#include "AliL3TrackArray.h"
#include "AliL3HoughMaxFinder.h"
#include "AliL3HoughTrack.h"
#include "AliL3SpacePointData.h"

#include <TH2.h>

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

//_____________________________________________________________
// AliL3HoughTransformerGlobal
//
// Hough transformation class
//

ClassImp(AliL3HoughTransformerGlobal)

AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal()
{
  //default ctor
  fPadMin=0;
  fPadMax=0;
  fNActiveSlice=0;
  fMemHandler=0;
  fEvent=0;
  fPsi=0;
  fPtMin=0;
  fTracks=0;
  fPeakFinder=0;
  fSeedPadRow=-1;
}

AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal(Char_t *path,Int_t event)
{
  //normal ctor
  strcpy(fPath,path);
  fEvent=event;
  fMemHandler=0;
  fPadMin=0;
  fPadMax=0;
  fPsi=0;
  fPtMin=0;
  fNActiveSlice=0;
  fTracks = new AliL3TrackArray("AliL3HoughTrack");
  fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
}

AliL3HoughTransformerGlobal::~AliL3HoughTransformerGlobal()
{
  //dtor
  if(fPadMin)
    delete [] fPadMin;
  if(fPadMax)
    delete [] fPadMax;
  if(fTracks)
    delete fTracks;
  if(fPeakFinder)
    delete fPeakFinder;
  UnloadActiveSlices();
  
}

void AliL3HoughTransformerGlobal::CreateHistograms(Float_t /*ptmin*/,Int_t nxbin,Int_t nybin)
{
  //create histograms for HT
  if(fPsi==0)
    {
      cerr<<"AliL3HoughTransformerGlobal::CreateHistograms : Call DefineRegion first"<<endl;
      exit(5);
    }
  AliL3HoughTransformer::CreateHistograms(nxbin,fPtMin,nybin,-fPsi,fPsi);
  //AliL3HoughTransformer::CreateHistograms(fPtMin,ptmax,ptres,nybin,fPsi);
}

void AliL3HoughTransformerGlobal::TransformCircleC()
{
  //Hough Transform
  AliL3Seed *clusters = new AliL3Seed[1000];
  Int_t nclusters = LoadClusterSeeds(clusters);
  
  Int_t i=0,sector,row,slice;
  UInt_t dummy=0;
  Float_t xyz[3],rpe[3],kappa,psi;
  
  for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
    {
      if(sl < 0 || sl < 18 && GetSlice() >= 18)
	slice = sl + 18;
      else if(sl > 35 || sl > 17 && GetSlice() <= 17)
	slice = sl-18;
      else
	slice = sl;
      cout<<"Transforming in slice "<<slice<<endl;
      
      AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
      
      for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	{
	  if(padrow == fSeedPadRow) continue;
	  AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
	  
	  for(UInt_t j=0; j<rowPt->fNDigit; j++)
	    {
	      Int_t pad = digits[j].fPad;

	      if(i==1 && pad < fPadMin[padrow])
		continue;
	      if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
		continue;

	      Int_t time = digits[j].fTime;
	      Int_t charge = digits[j].fCharge;
	      
	      if(charge > GetUpperThreshold() || charge < GetLowerThreshold())
		continue;
	      
	      AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	      AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	      
	      Rotate(xyz,i-1-fNActiveSlice);
	      
	      AliL3Transform::XYZtoRPhiEta(rpe,xyz);
	      rpe[0] = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
	      
	      Int_t etaindex = GetEtaIndex(rpe[2]);
	      
	      AliL3Histogram *hist = GetHistogram(etaindex);
	      
	      for(Int_t l=0; l<nclusters; l++)
		{
		  if(clusters[l].fIndex != etaindex) continue;
		  psi = atan( (clusters[l].fRadius*sin(rpe[1])-rpe[0]*sin(clusters[l].fPhi))/
			      (clusters[l].fRadius*cos(rpe[1])-rpe[0]*cos(clusters[l].fPhi)) );
		  kappa = 2*sin(clusters[l].fPhi-psi)/clusters[l].fRadius;
		  if(fabs(kappa) < 1e-33)
  		    continue;
		  //hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
		  hist->Fill(kappa,psi,charge);
		}
	    }
	  AliL3MemHandler::UpdateRowPointer(rowPt);
	}
    }
  delete [] clusters;
}

void AliL3HoughTransformerGlobal::TransformCircle()
{
  //Hough Transform  
  Int_t i=0,sector,row,slice;
  UInt_t dummy=0;
  Float_t xyz[3],rpe[3],kappa,psi;
  
  for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
    {
      if(sl < 0 || sl < 18 && GetSlice() >= 18)
	slice = sl + 18;
      else if(sl > 35 || sl > 17 && GetSlice() <= 17)
	slice = sl-18;
      else
	slice = sl;
      cout<<"Transforming in slice "<<slice<<endl;
      
      AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
      
      for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	{

	  AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
	  
	  for(UInt_t j=0; j<rowPt->fNDigit; j++)
	    {
	      Int_t pad = digits[j].fPad;
	      
	      if(i==1 && pad < fPadMin[padrow])
		continue;
	      if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
		continue;
	      
	      Int_t time = digits[j].fTime;
	      Int_t charge = digits[j].fCharge;
	      
	      if(charge > GetUpperThreshold() || charge < GetLowerThreshold())
		continue;
	      
	      AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	      AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	      
	      Rotate(xyz,i-1-fNActiveSlice);
	      
	      AliL3Transform::XYZtoRPhiEta(rpe,xyz);
	      rpe[0] = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
	      
	      Int_t etaindex = GetEtaIndex(rpe[2]);
	      
	      AliL3Histogram *hist = GetHistogram(etaindex);
	      
	      for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
		{
		  psi = hist->GetBinCenterY(b);
		  kappa = 2*sin(rpe[1] - psi)/rpe[0];
		  
		  if(fabs(kappa) < 1e-33)
		    continue;
		  
		  hist->Fill(kappa,psi,charge);
		  //hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
		}
	    }
	  AliL3MemHandler::UpdateRowPointer(rowPt);
	}
    }
}

void AliL3HoughTransformerGlobal::VerifyTracks(AliL3TrackArray *tracks,Int_t &index)
{
  //remove tracks which do not cross the seed padrow
  for(int i=index; i<tracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!tr) continue;

      Float_t angle=0;
      AliL3Transform::Local2GlobalAngle(&angle,GetSlice());
      if(!tr->CalculateReferencePoint(angle,AliL3Transform::Row2X(fSeedPadRow)))
	{
	  tracks->Remove(i);
	  continue;
	}
      Int_t sector,row;
      AliL3Transform::Slice2Sector(GetSlice(),fSeedPadRow,sector,row);
      float xyz[3] = {tr->GetPointX(),tr->GetPointY(),tr->GetPointZ()};
      AliL3Transform::Global2Raw(xyz,sector,row);
      if(xyz[1]<0 || xyz[1]>=AliL3Transform::GetNPads(fSeedPadRow) ||
	 xyz[2]<0 || xyz[2]>=AliL3Transform::GetNTimeBins())
	tracks->Remove(i);
    }
  tracks->Compress();
  index = tracks->GetNTracks();
}

void AliL3HoughTransformerGlobal::FindPeaks(AliL3Histogram *hist,Float_t eta)
{
  //Peak Finder  
  fPeakFinder->Reset();
  fPeakFinder->SetHistogram(hist);
  fPeakFinder->FindAbsMaxima();
  if(fPeakFinder->GetWeight(0) < 1000)
    return;
  
  AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks->NextTrack();
  track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
  track->SetEta(eta);
  track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
}

void AliL3HoughTransformerGlobal::Rotate(Float_t *xyz,Int_t relslice)
{
  //Rotate coordinates from one slice to the slice relative to it;
  //-1 means lower
  //1 means upper

  Float_t angle = (Float_t)relslice*AliL3Transform::Deg2Rad(20);
  Float_t x=xyz[0],y=xyz[1];
  xyz[0] = x*cos(angle) - y*sin(angle);
  xyz[1] = x*sin(angle) + y*cos(angle);
}


void AliL3HoughTransformerGlobal::DefineRegion(Float_t minpt,Float_t /*linephi*/,Int_t seedpadrow)
{
  //Setup the region to be included in the transform
  //This function should be called only once, since it is the same for all slices.
  
  fSeedPadRow=seedpadrow;
  fPtMin = minpt;
  
  fPadMin = new Int_t[AliL3Transform::GetNRows()];
  fPadMax = new Int_t[AliL3Transform::GetNRows()];
  
  //Phirange of data which should be included in the transform
  //Here we assume a min pt to define it; fPtMin
  //Based on pt (kappa) and the two points; origo and middle point, we can get psi angle
  
  //Calculate the upper point, the lower point is symmetrical around the x-axis
  //The track in this calculation has a positive charge, and the kappa sign is the opposite:  
  
  Float_t xyz[3],phi,phi2;
  
  phi = CalculateBorder(xyz,-1); //upward bending track
  phi2 = CalculateBorder(xyz,1); //downward bending track
  if(phi2 > phi)
    phi = phi2;
  
  cout<<"Phiangle "<<phi*180/AliL3Transform::Pi()<<" psi "<<fPsi*180/AliL3Transform::Pi()<<" nslices "<<fNActiveSlice<<endl;
  
  Float_t rotangle = fNActiveSlice*20;
  
  //Calculate the LUT for min/max pad for every padrow, and check which slices we need.
  Int_t pad,sector,row;
  for(Int_t i=0; i<AliL3Transform::GetNRows(); i++)
    {
      AliL3Transform::Slice2Sector(0,i,sector,row);
      
      //Lower boundary:
      pad = AliL3Transform::GetNPads(i)-1;
      while(pad >= 0)
	{
	  AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
	  if(AliL3Transform::GetPhi(xyz) > -1.*phi + AliL3Transform::Deg2Rad(rotangle))
	    fPadMin[i] = pad;
	  else
	    break;
	  pad--;
	}
      
      //Upper boundary
      pad = 0;
      while(pad < AliL3Transform::GetNPads(i))
	{
	  AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
	  if(AliL3Transform::GetPhi(xyz) < phi - AliL3Transform::Deg2Rad(rotangle))
	    fPadMax[i] = pad;
	  else
	    break;
	  pad++;
	}
    }
  
  cout<<"Padmax "<<fPadMax[155]<<endl;

}

Float_t AliL3HoughTransformerGlobal::CalculateBorder(Float_t *xyz,Int_t charge)
{
  //Define Hough space borders
  Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(fSeedPadRow),2) + pow(AliL3Transform::GetMaxY(fSeedPadRow),2));

  Double_t kappa = -charge*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/fPtMin;

  //Calculate the psi angle of the track = emission angle with x-axis in origo
  if( fabs(lineradius*kappa/2) > 1)
    {
      cerr<<"AliL3TransformerGlobal::DefineRegion : Angle too big"<<endl;
      exit(5);
    }
  
  Int_t nslices=0;
  Float_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
  if(charge > 0)
    fPsi=psi;
  
  cout<<"Calculated psi-angle "<<fPsi<<endl;
  
  AliL3Track track;
  track.SetFirstPoint(0,0,0);
  track.SetPt(fPtMin);
  track.SetPsi(psi);
  track.SetCharge(charge);
  track.CalculateHelix();
  
  Int_t crossingrow=0;
  if(charge < 0)
    crossingrow = AliL3Transform::GetNRows()-1;
  
  Float_t rotangle;
  
 redefine:
  
  rotangle = nslices*20;
  
  if(nslices==0)//here we are in local slice, so we use the appropriate function. a mess of course...
    {
      while(!track.GetCrossingPoint(crossingrow,xyz))
	{
	  if(crossingrow==0)
	    {	  
	      cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point1 on row "<<crossingrow<<endl;
	      exit(5);
	    }
	  crossingrow--;
	}
    }
  else
    {
      while(!track.CalculateReferencePoint(AliL3Transform::Deg2Rad(rotangle),AliL3Transform::Row2X(crossingrow)))
	{
	  if(crossingrow==0)
	    {
	      cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point2 on row "<<crossingrow<<endl;
	      exit(5);
	    }
	  crossingrow--;
	}
      
      xyz[0] = track.GetPointX();
      xyz[1] = track.GetPointY();
      xyz[2] = track.GetPointZ();
    }
  
  Float_t phi = atan2(xyz[1],xyz[0]);
  
  //Rotate coordinates back to local coordinates:
  Rotate(xyz,-1*nslices);
  
  Int_t sector,row;
  AliL3Transform::Slice2Sector(0,crossingrow,sector,row);
  AliL3Transform::Local2Raw(xyz,sector,row);
  if(xyz[1] < 0)
    {
      if(crossingrow>0)
	{
	  cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Wrong pad, probably a deadzone... "<<xyz[1]<<endl;
	  exit(5);
	}
      else
	return 0;//here you only want the crossing point with the outer padrow
    }
  if(xyz[1] >= AliL3Transform::GetNPads(crossingrow)) //Here the range is even one more slice away
    {
      cerr<<"One more slice with pad "<<xyz[1]<<endl;
      nslices++;
      if(charge < 0)
	crossingrow = AliL3Transform::GetNRows()-1;
      goto redefine;
    }
  
  if(nslices > fNActiveSlice)
    fNActiveSlice = nslices;
  
  cout<<"Calculated phi "<<phi<<" and pad "<<xyz[1]<<" on row "<<crossingrow<<endl;
  return phi;
}

void AliL3HoughTransformerGlobal::LoadActiveSlices(Bool_t binary)
{
  //Load active slices
  if(fMemHandler)
    UnloadActiveSlices();
  fMemHandler = new AliL3MemHandler*[(fNActiveSlice*2 + 1)];
  Char_t filename[1024];
  UInt_t dummy;
  Int_t i=0,slice;
  for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
    {
      if(sl < 0 || sl < 18 && GetSlice() >= 18)
	slice = sl + 18;
      else if(sl > 35 || sl > 17 && GetSlice() <= 17)
	slice = sl-18;
      else
	slice = sl;
      cout<<"Loading slice "<<slice<<endl;
      fMemHandler[i] = new AliL3FileHandler();
      fMemHandler[i]->Init(slice,-1);
      if(binary)
	{
	  sprintf(filename,"%s/binaries/digits_%d_%d_%d.raw",fPath,fEvent,slice,-1);
	  fMemHandler[i]->SetBinaryInput(filename);
	  fMemHandler[i]->CompBinary2Memory(dummy);
	  fMemHandler[i++]->CloseBinaryInput();
	}
      else
	{
	  sprintf(filename,"%s/digitfile.root",fPath);
	  fMemHandler[i]->SetAliInput(filename);
	  fMemHandler[i]->AliAltroDigits2Memory(dummy,fEvent);
	  fMemHandler[i++]->CloseAliInput();
	}
    }
  
}

void AliL3HoughTransformerGlobal::UnloadActiveSlices()
{
  //Unload active slices
  if(!fMemHandler)
    return;
  for(Int_t i=0; i<=fNActiveSlice*2; i++)
    {
      if(!fMemHandler[i]) continue;
      delete fMemHandler[i];
    }
  delete [] fMemHandler;
  fMemHandler=0;
}

Int_t AliL3HoughTransformerGlobal::LoadClusterSeeds(AliL3Seed *seeds)
{
  //Load cluster seeds
  Char_t filename[1024];
  UInt_t npoints;
  sprintf(filename,"%s/hough/points_%d_%d_%d.raw",fPath,fEvent,GetSlice(),-1);
  //sprintf(filename,"../tracker/points_%d_%d_%d.raw",fEvent,GetSlice(),-1);
  AliL3MemHandler mem;
  if(!mem.SetBinaryInput(filename))
    {
      cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Cannot open file "<<filename<<endl;
      exit(5);
    }
  AliL3SpacePointData *points = (AliL3SpacePointData*)mem.Allocate();
  mem.Binary2Memory(npoints,points);
  mem.CloseBinaryInput();
  
  Int_t counter=0;
  for(UInt_t i=0; i<npoints; i++)
    {
      Int_t lrow = (Int_t)points[i].fPadRow;
      if(lrow < fSeedPadRow) continue;
      if(lrow > fSeedPadRow) break;
      if(fabs(points[i].fY) > AliL3Transform::GetMaxY(lrow))
	{
	  cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Seeds seems not to be in local coordinates: "
	      <<points[i].fY<<endl;
	  exit(5);
	}
      Float_t xyz[3] = {AliL3Transform::Row2X(lrow),points[i].fY,points[i].fZ};
      Double_t eta = AliL3Transform::GetEta(xyz);
      seeds[counter].fPhi = atan2(xyz[1],xyz[0]);
      seeds[counter].fRadius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
      seeds[counter].fIndex = GetEtaIndex(eta);
      seeds[counter++].fEta = eta;
    }
  cout<<"Loaded "<<counter<<" cluster seeds on slice "<<GetSlice()<<" padrow "<<fSeedPadRow<<endl;
  mem.Free();
  return counter;
}

void AliL3HoughTransformerGlobal::DisplayActiveRegion(TH2F *hist,Int_t etaindex)
{
  //Fill the active region in a histogram
  
  Int_t i=0,sector,row,slice;
  UInt_t dummy=0;
  Float_t xyz[3];
  for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
    {
      if(sl < 0 || sl < 18 && GetSlice() >= 18)
	slice = sl + 18;
      else if(sl > 35 || sl > 17 && GetSlice() <= 17)
	slice = sl-18;
      else
	slice = sl;
      cout<<"Displaying slice "<<slice<<endl;
      AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
      for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	{
	  AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
	  for(UInt_t j=0; j<rowPt->fNDigit; j++)
	    {
	      Int_t pad = digits[j].fPad;
	      if(i==1 && pad < fPadMin[padrow])
 		continue;
	      if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
		continue;
	      Int_t time = digits[j].fTime;
	      AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	      AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	     
	      Rotate(xyz,i-1-fNActiveSlice);
	      
	      Double_t eta = AliL3Transform::GetEta(xyz);
	      if(GetEtaIndex(eta) == etaindex)
		hist->Fill(xyz[0],xyz[1]);
	    }
	  AliL3MemHandler::UpdateRowPointer(rowPt);
	}
    }
}
Commit	Line	Data
3e87ef69	1	// @(#) $Id$
	2
	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	4	//*-- Copyright &copy ALICE HLT Group
	5
	6	#include "AliL3StandardIncludes.h"
	7
	8	#include "AliL3Logging.h"
	9	#include "AliL3HoughTransformerGlobal.h"
	10	#include "AliL3FileHandler.h"
	11	#include "AliL3Transform.h"
	12	#include "AliL3DigitData.h"
	13	#include "AliL3Track.h"
	14	#include "AliL3Histogram.h"
	15	#include "AliL3TrackArray.h"
	16	#include "AliL3HoughMaxFinder.h"
	17	#include "AliL3HoughTrack.h"
	18	#include "AliL3SpacePointData.h"
	19
	20	#include <TH2.h>
	21
0bd0c1ef	22	#if __GNUC__ == 3
3e87ef69	23	using namespace std;
	24	#endif
	25
	26	//_____________________________________________________________
	27	// AliL3HoughTransformerGlobal
	28	//
	29	// Hough transformation class
	30	//
	31
	32	ClassImp(AliL3HoughTransformerGlobal)
	33
	34	AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal()
	35	{
c62b480b	36	//default ctor
3e87ef69	37	fPadMin=0;
	38	fPadMax=0;
	39	fNActiveSlice=0;
	40	fMemHandler=0;
	41	fEvent=0;
	42	fPsi=0;
	43	fPtMin=0;
	44	fTracks=0;
	45	fPeakFinder=0;
	46	fSeedPadRow=-1;
	47	}
	48
	49	AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal(Char_t *path,Int_t event)
	50	{
c62b480b	51	//normal ctor
3e87ef69	52	strcpy(fPath,path);
	53	fEvent=event;
	54	fMemHandler=0;
	55	fPadMin=0;
	56	fPadMax=0;
	57	fPsi=0;
	58	fPtMin=0;
	59	fNActiveSlice=0;
	60	fTracks = new AliL3TrackArray("AliL3HoughTrack");
	61	fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
	62	}
	63
	64	AliL3HoughTransformerGlobal::~AliL3HoughTransformerGlobal()
	65	{
c62b480b	66	//dtor
3e87ef69	67	if(fPadMin)
	68	delete [] fPadMin;
	69	if(fPadMax)
	70	delete [] fPadMax;
	71	if(fTracks)
	72	delete fTracks;
	73	if(fPeakFinder)
	74	delete fPeakFinder;
	75	UnloadActiveSlices();
	76
	77	}
	78
dd7d3870	79	void AliL3HoughTransformerGlobal::CreateHistograms(Float_t /ptmin/,Int_t nxbin,Int_t nybin)
3e87ef69	80	{
c62b480b	81	//create histograms for HT
3e87ef69	82	if(fPsi==0)
	83	{
	84	cerr<<"AliL3HoughTransformerGlobal::CreateHistograms : Call DefineRegion first"<<endl;
	85	exit(5);
	86	}
b2a02bce	87	AliL3HoughTransformer::CreateHistograms(nxbin,fPtMin,nybin,-fPsi,fPsi);
b2a02bce	88	//AliL3HoughTransformer::CreateHistograms(fPtMin,ptmax,ptres,nybin,fPsi);
3e87ef69	89	}
	90
	91	void AliL3HoughTransformerGlobal::TransformCircleC()
	92	{
c62b480b	93	//Hough Transform
c62b480b	94	AliL3Seed *clusters = new AliL3Seed[1000];
3e87ef69	95	Int_t nclusters = LoadClusterSeeds(clusters);
	96
	97	Int_t i=0,sector,row,slice;
	98	UInt_t dummy=0;
	99	Float_t xyz[3],rpe[3],kappa,psi;
	100
	101	for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
	102	{
	103	if(sl < 0 \|\| sl < 18 && GetSlice() >= 18)
	104	slice = sl + 18;
	105	else if(sl > 35 \|\| sl > 17 && GetSlice() <= 17)
	106	slice = sl-18;
	107	else
	108	slice = sl;
	109	cout<<"Transforming in slice "<<slice<<endl;
	110
	111	AliL3DigitRowData rowPt = (AliL3DigitRowData)fMemHandler[i++]->GetDataPointer(dummy);
	112
	113	for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	114	{
	115	if(padrow == fSeedPadRow) continue;
	116	AliL3DigitData digits = (AliL3DigitData)rowPt->fDigitData;
	117
	118	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	119	{
	120	Int_t pad = digits[j].fPad;
b2a02bce	121
3e87ef69	122	if(i==1 && pad < fPadMin[padrow])
	123	continue;
	124	if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
	125	continue;
b2a02bce	126
3e87ef69	127	Int_t time = digits[j].fTime;
	128	Int_t charge = digits[j].fCharge;
	129
	130	if(charge > GetUpperThreshold() \|\| charge < GetLowerThreshold())
	131	continue;
	132
	133	AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	134	AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	135
	136	Rotate(xyz,i-1-fNActiveSlice);
	137
	138	AliL3Transform::XYZtoRPhiEta(rpe,xyz);
	139	rpe[0] = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	140
c62b480b	141	Int_t etaindex = GetEtaIndex(rpe[2]);
3e87ef69	142
c62b480b	143	AliL3Histogram *hist = GetHistogram(etaindex);
3e87ef69	144
	145	for(Int_t l=0; l<nclusters; l++)
	146	{
c62b480b	147	if(clusters[l].fIndex != etaindex) continue;
3e87ef69	148	psi = atan( (clusters[l].fRadiussin(rpe[1])-rpe[0]sin(clusters[l].fPhi))/
	149	(clusters[l].fRadiuscos(rpe[1])-rpe[0]cos(clusters[l].fPhi)) );
	150	kappa = 2*sin(clusters[l].fPhi-psi)/clusters[l].fRadius;
	151	if(fabs(kappa) < 1e-33)
	152	continue;
b2a02bce	153	//hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
b2a02bce	154	hist->Fill(kappa,psi,charge);
3e87ef69	155	}
	156	}
	157	AliL3MemHandler::UpdateRowPointer(rowPt);
	158	}
	159	}
	160	delete [] clusters;
	161	}
	162
	163	void AliL3HoughTransformerGlobal::TransformCircle()
	164	{
c62b480b	165	//Hough Transform
3e87ef69	166	Int_t i=0,sector,row,slice;
	167	UInt_t dummy=0;
	168	Float_t xyz[3],rpe[3],kappa,psi;
	169
	170	for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
	171	{
	172	if(sl < 0 \|\| sl < 18 && GetSlice() >= 18)
	173	slice = sl + 18;
	174	else if(sl > 35 \|\| sl > 17 && GetSlice() <= 17)
	175	slice = sl-18;
	176	else
	177	slice = sl;
	178	cout<<"Transforming in slice "<<slice<<endl;
	179
	180	AliL3DigitRowData rowPt = (AliL3DigitRowData)fMemHandler[i++]->GetDataPointer(dummy);
	181
	182	for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	183	{
	184
	185	AliL3DigitData digits = (AliL3DigitData)rowPt->fDigitData;
	186
	187	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	188	{
	189	Int_t pad = digits[j].fPad;
	190
	191	if(i==1 && pad < fPadMin[padrow])
	192	continue;
	193	if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
	194	continue;
	195
	196	Int_t time = digits[j].fTime;
	197	Int_t charge = digits[j].fCharge;
	198
	199	if(charge > GetUpperThreshold() \|\| charge < GetLowerThreshold())
	200	continue;
	201
	202	AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	203	AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	204
	205	Rotate(xyz,i-1-fNActiveSlice);
	206
	207	AliL3Transform::XYZtoRPhiEta(rpe,xyz);
	208	rpe[0] = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	209
c62b480b	210	Int_t etaindex = GetEtaIndex(rpe[2]);
3e87ef69	211
c62b480b	212	AliL3Histogram *hist = GetHistogram(etaindex);
3e87ef69	213
	214	for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
	215	{
	216	psi = hist->GetBinCenterY(b);
	217	kappa = 2*sin(rpe[1] - psi)/rpe[0];
	218
	219	if(fabs(kappa) < 1e-33)
	220	continue;
	221
b2a02bce	222	hist->Fill(kappa,psi,charge);
b2a02bce	223	//hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
3e87ef69	224	}
	225	}
	226	AliL3MemHandler::UpdateRowPointer(rowPt);
	227	}
	228	}
	229	}
	230
	231	void AliL3HoughTransformerGlobal::VerifyTracks(AliL3TrackArray *tracks,Int_t &index)
	232	{
c62b480b	233	//remove tracks which do not cross the seed padrow
3e87ef69	234	for(int i=index; i<tracks->GetNTracks(); i++)
	235	{
	236	AliL3HoughTrack tr = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	237	if(!tr) continue;
	238
	239	Float_t angle=0;
	240	AliL3Transform::Local2GlobalAngle(&angle,GetSlice());
	241	if(!tr->CalculateReferencePoint(angle,AliL3Transform::Row2X(fSeedPadRow)))
	242	{
	243	tracks->Remove(i);
	244	continue;
	245	}
	246	Int_t sector,row;
	247	AliL3Transform::Slice2Sector(GetSlice(),fSeedPadRow,sector,row);
	248	float xyz[3] = {tr->GetPointX(),tr->GetPointY(),tr->GetPointZ()};
	249	AliL3Transform::Global2Raw(xyz,sector,row);
	250	if(xyz[1]<0 \|\| xyz[1]>=AliL3Transform::GetNPads(fSeedPadRow) \|\|
	251	xyz[2]<0 \|\| xyz[2]>=AliL3Transform::GetNTimeBins())
	252	tracks->Remove(i);
	253	}
	254	tracks->Compress();
	255	index = tracks->GetNTracks();
	256	}
	257
	258	void AliL3HoughTransformerGlobal::FindPeaks(AliL3Histogram *hist,Float_t eta)
	259	{
c62b480b	260	//Peak Finder
3e87ef69	261	fPeakFinder->Reset();
	262	fPeakFinder->SetHistogram(hist);
	263	fPeakFinder->FindAbsMaxima();
	264	if(fPeakFinder->GetWeight(0) < 1000)
	265	return;
	266
	267	AliL3HoughTrack track = (AliL3HoughTrack)fTracks->NextTrack();
	268	track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
	269	track->SetEta(eta);
	270	track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
	271	}
	272
c62b480b	273	void AliL3HoughTransformerGlobal::Rotate(Float_t *xyz,Int_t relslice)
3e87ef69	274	{
	275	//Rotate coordinates from one slice to the slice relative to it;
	276	//-1 means lower
	277	//1 means upper
	278
c62b480b	279	Float_t angle = (Float_t)relslice*AliL3Transform::Deg2Rad(20);
3e87ef69	280	Float_t x=xyz[0],y=xyz[1];
	281	xyz[0] = xcos(angle) - ysin(angle);
	282	xyz[1] = xsin(angle) + ycos(angle);
	283	}
	284
	285
dd7d3870	286	void AliL3HoughTransformerGlobal::DefineRegion(Float_t minpt,Float_t /linephi/,Int_t seedpadrow)
3e87ef69	287	{
	288	//Setup the region to be included in the transform
	289	//This function should be called only once, since it is the same for all slices.
	290
	291	fSeedPadRow=seedpadrow;
	292	fPtMin = minpt;
	293
	294	fPadMin = new Int_t[AliL3Transform::GetNRows()];
	295	fPadMax = new Int_t[AliL3Transform::GetNRows()];
	296
	297	//Phirange of data which should be included in the transform
	298	//Here we assume a min pt to define it; fPtMin
	299	//Based on pt (kappa) and the two points; origo and middle point, we can get psi angle
	300
	301	//Calculate the upper point, the lower point is symmetrical around the x-axis
	302	//The track in this calculation has a positive charge, and the kappa sign is the opposite:
	303
	304	Float_t xyz[3],phi,phi2;
	305
	306	phi = CalculateBorder(xyz,-1); //upward bending track
	307	phi2 = CalculateBorder(xyz,1); //downward bending track
	308	if(phi2 > phi)
	309	phi = phi2;
	310
	311	cout<<"Phiangle "<<phi180/AliL3Transform::Pi()<<" psi "<<fPsi180/AliL3Transform::Pi()<<" nslices "<<fNActiveSlice<<endl;
	312
	313	Float_t rotangle = fNActiveSlice*20;
	314
	315	//Calculate the LUT for min/max pad for every padrow, and check which slices we need.
	316	Int_t pad,sector,row;
	317	for(Int_t i=0; i<AliL3Transform::GetNRows(); i++)
	318	{
	319	AliL3Transform::Slice2Sector(0,i,sector,row);
	320
	321	//Lower boundary:
	322	pad = AliL3Transform::GetNPads(i)-1;
	323	while(pad >= 0)
	324	{
	325	AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
	326	if(AliL3Transform::GetPhi(xyz) > -1.*phi + AliL3Transform::Deg2Rad(rotangle))
	327	fPadMin[i] = pad;
	328	else
	329	break;
	330	pad--;
	331	}
	332
	333	//Upper boundary
	334	pad = 0;
	335	while(pad < AliL3Transform::GetNPads(i))
	336	{
	337	AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
	338	if(AliL3Transform::GetPhi(xyz) < phi - AliL3Transform::Deg2Rad(rotangle))
	339	fPadMax[i] = pad;
	340	else
	341	break;
	342	pad++;
	343	}
	344	}
	345
	346	cout<<"Padmax "<<fPadMax[155]<<endl;
	347
	348	}
	349
	350	Float_t AliL3HoughTransformerGlobal::CalculateBorder(Float_t *xyz,Int_t charge)
351	{
c62b480b	352	//Define Hough space borders
3e87ef69	353	Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(fSeedPadRow),2) + pow(AliL3Transform::GetMaxY(fSeedPadRow),2));
	354
	355	Double_t kappa = -chargeAliL3Transform::GetBField()AliL3Transform::GetBFact()/fPtMin;
	356
	357	//Calculate the psi angle of the track = emission angle with x-axis in origo
	358	if( fabs(lineradius*kappa/2) > 1)
	359	{
	360	cerr<<"AliL3TransformerGlobal::DefineRegion : Angle too big"<<endl;
	361	exit(5);
	362	}
	363
	364	Int_t nslices=0;
	365	Float_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
	366	if(charge > 0)
	367	fPsi=psi;
	368
b2a02bce	369	cout<<"Calculated psi-angle "<<fPsi<<endl;
b2a02bce	370
3e87ef69	371	AliL3Track track;
	372	track.SetFirstPoint(0,0,0);
	373	track.SetPt(fPtMin);
	374	track.SetPsi(psi);
	375	track.SetCharge(charge);
	376	track.CalculateHelix();
	377
	378	Int_t crossingrow=0;
	379	if(charge < 0)
	380	crossingrow = AliL3Transform::GetNRows()-1;
	381
	382	Float_t rotangle;
	383
	384	redefine:
	385
	386	rotangle = nslices*20;
	387
	388	if(nslices==0)//here we are in local slice, so we use the appropriate function. a mess of course...
	389	{
	390	while(!track.GetCrossingPoint(crossingrow,xyz))
	391	{
	392	if(crossingrow==0)
	393	{
	394	cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point1 on row "<<crossingrow<<endl;
	395	exit(5);
	396	}
	397	crossingrow--;
	398	}
	399	}
	400	else
	401	{
	402	while(!track.CalculateReferencePoint(AliL3Transform::Deg2Rad(rotangle),AliL3Transform::Row2X(crossingrow)))
	403	{
	404	if(crossingrow==0)
	405	{
	406	cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point2 on row "<<crossingrow<<endl;
	407	exit(5);
	408	}
	409	crossingrow--;
	410	}
	411
	412	xyz[0] = track.GetPointX();
	413	xyz[1] = track.GetPointY();
	414	xyz[2] = track.GetPointZ();
	415	}
	416
	417	Float_t phi = atan2(xyz[1],xyz[0]);
	418
	419	//Rotate coordinates back to local coordinates:
	420	Rotate(xyz,-1*nslices);
	421
	422	Int_t sector,row;
	423	AliL3Transform::Slice2Sector(0,crossingrow,sector,row);
	424	AliL3Transform::Local2Raw(xyz,sector,row);
	425	if(xyz[1] < 0)
	426	{
	427	if(crossingrow>0)
	428	{
	429	cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Wrong pad, probably a deadzone... "<<xyz[1]<<endl;
	430	exit(5);
	431	}
	432	else
	433	return 0;//here you only want the crossing point with the outer padrow
	434	}
435	if(xyz[1] >= AliL3Transform::GetNPads(crossingrow)) //Here the range is even one more slice away
436	{
437	cerr<<"One more slice with pad "<<xyz[1]<<endl;
438	nslices++;
439	if(charge < 0)
440	crossingrow = AliL3Transform::GetNRows()-1;
441	goto redefine;
442	}
443
444	if(nslices > fNActiveSlice)
445	fNActiveSlice = nslices;
446
447	cout<<"Calculated phi "<<phi<<" and pad "<<xyz[1]<<" on row "<<crossingrow<<endl;
448	return phi;
449	}
450
451	void AliL3HoughTransformerGlobal::LoadActiveSlices(Bool_t binary)
452	{
c62b480b	453	//Load active slices
3e87ef69	454	if(fMemHandler)
	455	UnloadActiveSlices();
	456	fMemHandler = new AliL3MemHandler[(fNActiveSlice2 + 1)];
	457	Char_t filename[1024];
	458	UInt_t dummy;
	459	Int_t i=0,slice;
	460	for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
	461	{
	462	if(sl < 0 \|\| sl < 18 && GetSlice() >= 18)
	463	slice = sl + 18;
	464	else if(sl > 35 \|\| sl > 17 && GetSlice() <= 17)
	465	slice = sl-18;
	466	else
	467	slice = sl;
	468	cout<<"Loading slice "<<slice<<endl;
	469	fMemHandler[i] = new AliL3FileHandler();
	470	fMemHandler[i]->Init(slice,-1);
	471	if(binary)
	472	{
	473	sprintf(filename,"%s/binaries/digits_%d_%d_%d.raw",fPath,fEvent,slice,-1);
	474	fMemHandler[i]->SetBinaryInput(filename);
	475	fMemHandler[i]->CompBinary2Memory(dummy);
	476	fMemHandler[i++]->CloseBinaryInput();
	477	}
	478	else
	479	{
b2a02bce	480	sprintf(filename,"%s/digitfile.root",fPath);
b2a02bce	481	fMemHandler[i]->SetAliInput(filename);
3e87ef69	482	fMemHandler[i]->AliAltroDigits2Memory(dummy,fEvent);
	483	fMemHandler[i++]->CloseAliInput();
	484	}
	485	}
	486
	487	}
	488
	489	void AliL3HoughTransformerGlobal::UnloadActiveSlices()
	490	{
c62b480b	491	//Unload active slices
3e87ef69	492	if(!fMemHandler)
	493	return;
	494	for(Int_t i=0; i<=fNActiveSlice*2; i++)
	495	{
	496	if(!fMemHandler[i]) continue;
	497	delete fMemHandler[i];
	498	}
	499	delete [] fMemHandler;
	500	fMemHandler=0;
	501	}
	502
c62b480b	503	Int_t AliL3HoughTransformerGlobal::LoadClusterSeeds(AliL3Seed *seeds)
3e87ef69	504	{
c62b480b	505	//Load cluster seeds
3e87ef69	506	Char_t filename[1024];
	507	UInt_t npoints;
	508	sprintf(filename,"%s/hough/points_%d_%d_%d.raw",fPath,fEvent,GetSlice(),-1);
	509	//sprintf(filename,"../tracker/points_%d_%d_%d.raw",fEvent,GetSlice(),-1);
	510	AliL3MemHandler mem;
	511	if(!mem.SetBinaryInput(filename))
	512	{
	513	cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Cannot open file "<<filename<<endl;
	514	exit(5);
	515	}
	516	AliL3SpacePointData points = (AliL3SpacePointData)mem.Allocate();
	517	mem.Binary2Memory(npoints,points);
	518	mem.CloseBinaryInput();
	519
	520	Int_t counter=0;
	521	for(UInt_t i=0; i<npoints; i++)
	522	{
	523	Int_t lrow = (Int_t)points[i].fPadRow;
	524	if(lrow < fSeedPadRow) continue;
	525	if(lrow > fSeedPadRow) break;
	526	if(fabs(points[i].fY) > AliL3Transform::GetMaxY(lrow))
	527	{
	528	cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Seeds seems not to be in local coordinates: "
	529	<<points[i].fY<<endl;
	530	exit(5);
	531	}
	532	Float_t xyz[3] = {AliL3Transform::Row2X(lrow),points[i].fY,points[i].fZ};
	533	Double_t eta = AliL3Transform::GetEta(xyz);
	534	seeds[counter].fPhi = atan2(xyz[1],xyz[0]);
	535	seeds[counter].fRadius = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	536	seeds[counter].fIndex = GetEtaIndex(eta);
	537	seeds[counter++].fEta = eta;
	538	}
	539	cout<<"Loaded "<<counter<<" cluster seeds on slice "<<GetSlice()<<" padrow "<<fSeedPadRow<<endl;
	540	mem.Free();
	541	return counter;
	542	}
	543
c62b480b	544	void AliL3HoughTransformerGlobal::DisplayActiveRegion(TH2F *hist,Int_t etaindex)
3e87ef69	545	{
	546	//Fill the active region in a histogram
	547
	548	Int_t i=0,sector,row,slice;
	549	UInt_t dummy=0;
	550	Float_t xyz[3];
	551	for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
	552	{
	553	if(sl < 0 \|\| sl < 18 && GetSlice() >= 18)
	554	slice = sl + 18;
	555	else if(sl > 35 \|\| sl > 17 && GetSlice() <= 17)
	556	slice = sl-18;
	557	else
	558	slice = sl;
	559	cout<<"Displaying slice "<<slice<<endl;
	560	AliL3DigitRowData rowPt = (AliL3DigitRowData)fMemHandler[i++]->GetDataPointer(dummy);
	561	for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
	562	{
	563	AliL3DigitData digits = (AliL3DigitData)rowPt->fDigitData;
	564	for(UInt_t j=0; j<rowPt->fNDigit; j++)
	565	{
	566	Int_t pad = digits[j].fPad;
	567	if(i==1 && pad < fPadMin[padrow])
	568	continue;
	569	if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
	570	continue;
	571	Int_t time = digits[j].fTime;
	572	AliL3Transform::Slice2Sector(slice,padrow,sector,row);
	573	AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
	574
	575	Rotate(xyz,i-1-fNActiveSlice);
	576
	577	Double_t eta = AliL3Transform::GetEta(xyz);
c62b480b	578	if(GetEtaIndex(eta) == etaindex)
3e87ef69	579	hist->Fill(xyz[0],xyz[1]);
	580	}
	581	AliL3MemHandler::UpdateRowPointer(rowPt);
	582	}
	583	}
	584	}