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

//Author:        Anders Strand Vestbo
//Last Modified: 28.6.01

#include <math.h>
#include <TH2.h>
#include <time.h>

#include <TFile.h>
#include <TTree.h>

#include "AliTPCParam.h"
#include "AliSimDigits.h"

#include "AliL3Histogram.h"
#include "AliL3Logging.h"
#include "AliL3Defs.h"
#include "AliL3Transform.h"
#include "AliL3HoughTransformer.h"
#include "AliL3HoughTrack.h"
#include "AliL3TrackArray.h"
#include "AliL3DigitData.h"

ClassImp(AliL3HoughTransformer)

AliL3HoughTransformer::AliL3HoughTransformer()
{
  //Default constructor
  
}


AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
{
  //Constructor
  
  fTransform = new AliL3Transform();

  fEtaMin = etarange[0];
  fEtaMax = etarange[1];
  fSlice = slice; 
  fPatch = patch;
  fNumOfPadRows=NRowsSlice;
}

AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
{
  
  fTransform = new AliL3Transform();
  if(etarange)
    {
      //assumes you want to look at a given etaslice only
      fEtaMin = etarange[0];
      fEtaMax = etarange[1];
    }
  else
    {
      //transform in all etaslices
      fEtaMin = 0;
      fEtaMax = slice < 18 ? 0.9 : -0.9;
    }
  fNumEtaSegments = n_eta_segments;
  fSlice = slice;
  fPatch = patch;
  fNumOfPadRows = NRowsSlice;
  
  fNRowsInPatch = NRows[fPatch][1]-NRows[fPatch][0] + 1;
  fBinTableBounds = (fNRowsInPatch+1)*(MaxNPads+1);
  fNDigitRowData = 0;
  fDigitRowData = 0;
  fHistoPt = 0;
}


AliL3HoughTransformer::~AliL3HoughTransformer()
{
  //Destructor
  if(fBinTable)
    {
      for(Int_t i=0; i<fBinTableBounds; i++)
	delete [] fBinTable[i];
      delete [] fBinTable;
    }
  if(fEtaIndex)
    delete [] fEtaIndex;
  if(fTrackTable)
    {
      for(Int_t i=0; i<fNumEtaSegments; i++)
	delete [] fTrackTable[i];
      delete [] fTrackTable;
    }
  if(fTransform)
    delete fTransform;
  
  
}

void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
{
  fNDigitRowData = ndigits;
  fDigitRowData = ptr;
}

void AliL3HoughTransformer::InitTables()
{
  //Create LUT for the circle transform.
  //the actual transformation is done in TransformTables.

  AliL3Histogram *hist = fHistoPt;
  
  Int_t nbinsy = hist->GetNbinsY();
  
  fBinTable = new Int_t*[fBinTableBounds];
  Int_t index;

  Double_t etaslice = fEtaMax/fNumEtaSegments;
  
  Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);

  fEtaIndex = new Int_t[etabounds];
  for(Int_t i=0; i<etabounds; i++)
    fEtaIndex[i] = -1;
  
  fTrackTable = new Char_t*[fNumEtaSegments];
  for(Int_t i=0; i<fNumEtaSegments; i++)
    {
      fTrackTable[i] = new Char_t[fBinTableBounds];
      for(Int_t j=0; j<fBinTableBounds; j++)
	fTrackTable[i][j] = 0;
    }
  
  for(Int_t r=NRows[fPatch][0]; r<=NRows[fPatch][1]; r++)
    {
      Int_t prow = r - NRows[fPatch][0];
      
      for(Int_t p=0; p<fTransform->GetNPads(r); p++)
	{
	  Float_t xyz[3];
	  Int_t sector,row;
	  for(Int_t t=0; t<fTransform->GetNTimeBins(); t++)
	    {
	      fTransform->Slice2Sector(fSlice,r,sector,row);
	      fTransform->Raw2Local(xyz,sector,row,p,t);
	      Double_t eta = fTransform->GetEta(xyz);
	      if(eta < fEtaMin || eta > fEtaMax) continue;
	      Int_t ind = (prow<<17) + (p<<9) + t;
	      if(fEtaIndex[ind]>=0)
		printf("AliL3HoughTransformer::InitTables : Overlapping indexes in eta!!\n");
	      Int_t eta_index = (Int_t)(eta/etaslice);	      
	      if(eta_index < 0 || eta_index > fNumEtaSegments)
		continue;
	      fEtaIndex[ind] = eta_index;
	      
	    }
	  
	  Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
	  Double_t phi_pix = fTransform->GetPhi(xyz);
	  index = (prow<<8) + p;
	  fBinTable[index] = new Int_t[nbinsy+1];
	  
	  for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
	    {
	      Double_t phi0 = hist->GetBinCenterY(b);
	      Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	      Int_t bin = hist->FindBin(kappa,phi0);
	      if(fBinTable[index][b]!=0)
		printf("AliL3HoughTransformer::InitTables : Overlapping indexes %d %d b %d\n",fBinTable[index][b],index,b);
	      fBinTable[index][b] = bin;
	    }
	}

    }

}

void AliL3HoughTransformer::TransformTables(AliL3Histogram **histos,AliL3Histogram **images)
{
  //Transform all the pixels while reading, and fill the corresponding histograms.
  //Transform is done using LUT created in InitTables.
  //fTrackTable : table telling whether a specific pixel is active (nonzero):
  //fTrackTable = 0  ->  no track
  //fTrackindex > 0  ->  track present
  //fTrackindex = -1  ->  track has been removed (already found)
  //fEtaIndex : table storing the etaindex -> used to find correct histogram to fill
  //fBinTable : table storing all the bins to fill for each nonzero pixel
    
  Int_t eta_index;
  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
  AliL3Histogram *hist;
  
  if(!tempPt)
    {
      LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformTables","data")<<
	"Zero datapointer"<<ENDLOG;
      return;
    }

  Int_t out_count=0,tot_count=0;
  Int_t index,ind;

  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      Int_t prow = i - NRows[fPatch][0];
      AliL3DigitData *bins = tempPt->fDigitData;
      for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
	{
	  index = (prow<<8) + bins[dig].fPad;
	  ind = (prow<<17) + (bins[dig].fPad<<9) + bins[dig].fTime;
	  eta_index = fEtaIndex[ind];
	  if(eta_index < 0) continue;  //pixel out of etarange
	  
	  if(fTrackTable[eta_index][index]<0) continue; //this pixel has already been removed. 
	  fTrackTable[eta_index][index]++; //this pixel is now marked as active (it is nonzero)
	  
	  tot_count++;
	  hist = histos[eta_index];
	  
	  if(images)
	    {
	      //Display the transformed images.
	      AliL3Histogram *image = images[eta_index];
	      Float_t xyz_local[3];
	      Int_t sector,row;
	      fTransform->Slice2Sector(fSlice,i,sector,row);
	      fTransform->Raw2Local(xyz_local,sector,row,bins[dig].fPad,bins[dig].fTime);
	      image->Fill(xyz_local[0],xyz_local[1],bins[dig].fCharge);
	    }
	  
	  if(!hist)
	    {
	      printf("Error getting histogram!\n");
	      continue;
	    }
	  for(Int_t p=hist->GetFirstYbin(); p<=hist->GetLastYbin(); p++)
	    hist->AddBinContent(fBinTable[index][p],bins[dig].fCharge);
	  
	}
      
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
    }  
  
}

void AliL3HoughTransformer::WriteTables()
{
  //Write the tables to asci files.
  
  AliL3Histogram *hist = fHistoPt;
  Char_t name[100];
  sprintf(name,"histogram_table_%d.txt",fPatch);
  FILE *file = fopen(name,"w");
  
  Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);
  for(Int_t i=0; i<etabounds; i++)
    {
      if(fEtaIndex[i]<0) continue;
      fprintf(file,"%d %d\n",i,fEtaIndex[i]);
    }
  fclose(file);
  
  sprintf(name,"bin_table_%d.txt",fPatch);
  FILE *file2 = fopen(name,"w");
  for(Int_t i=0; i<fBinTableBounds; i++)
    {
      if(!fBinTable[i]) continue;
      fprintf(file2,"%d ",i);
      for(Int_t j=hist->GetFirstYbin(); j<=hist->GetLastYbin(); j++)
	{
	  fprintf(file2,"%d ",fBinTable[i][j]);
	}
      fprintf(file2,"\n");
    }
  fclose(file2);
}

Double_t AliL3HoughTransformer::CpuTime()
{
  return (Double_t)(clock()) / CLOCKS_PER_SEC;
}

/*
void AliL3HoughTransformer::InitTemplates(TH2F *hist)
{

  AliL3Digits *pixel;

  Int_t ymin = hist->GetYaxis()->GetFirst();
  Int_t ymax = hist->GetYaxis()->GetLast();
  Int_t nbinsy = hist->GetNbinsY();

  fIndex = new Int_t*[fNDigits];
  for(Int_t i=0; i<fNDigits; i++)
    fIndex[i] = new Int_t[nbinsy+1];
    
  Int_t sector,row;
  
  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
    {        
      
      for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
	{
	  Float_t xyz[3];
	  fTransform->Slice2Sector(fSlice,padrow,sector,row);
	  fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
	  
	  Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
	  
	  Double_t phi_pix = fTransform->GetPhi(xyz);
	  Int_t index = pixel->fIndex;
	  if(index >= fNDigits)
	    printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
	  for(Int_t p=ymin; p<=ymax; p++)
	    {
	      
	      Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	      Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	      //printf("kappa %f phi0 %f\n",kappa,phi0);
	      Int_t bin = hist->FindBin(kappa,phi0);
	      if(fIndex[index][p]!=0)
		printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
	      fIndex[index][p] = bin;
	    }
	}
    }
  
}


void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
{
  //Transformation is done with respect to local coordinates in slice.
  //Transform every pixel into whole phirange, using parametrisation:
  //kappa = 2*sin(phi-phi0)/R
  //Assumes you run InitTemplates first!!!!

  AliL3Digits *pix1;
    
  Int_t nbinsy = hist->GetNbinsY();

  Int_t totsignal=0,vol_index;
  if(fNumEtaSegments==1)
    eta_index=0; //only looking in one etaslice.

  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
    {
      vol_index = eta_index*fNumOfPadRows + padrow;
      //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
      for(pix1=(AliL3Digits*)fVolume[vol_index].first; pix1!=0; pix1=(AliL3Digits*)pix1->fNextVolumePixel)
	{
	  Float_t xyz[3];
	  fTransform->Raw2Global(xyz,2,padrow,pix1->fPad,pix1->fTime);
	  Double_t eta = fTransform->GetEta(xyz);
	  if(eta < fEtaMin || eta > fEtaMax)
	    printf("\n Eta OUT OF RANGE\n");

	  Short_t signal = pix1->fCharge;
	  Int_t index = pix1->fIndex;
	  if(index < 0) continue; //This pixel has been removed.
	  totsignal += signal;
	  for(Int_t p=0; p <= nbinsy; p++)
	    hist->AddBinContent(fIndex[index][p],signal);
	  		  
	}
      
    }
    
  printf("Total signal %d\n",totsignal);
}

void AliL3HoughTransformer::Transform2Circle(TH2F **histos,Int_t n_eta_segments,UInt_t ndigits,AliL3DigitRowData *ptr)
{
  //Transform all the pixels while reading them, and fill the corresponding histograms.
  //Everything is done in one go here. 

  Double_t etaslice = 0.9/n_eta_segments;
  Int_t eta_index;
  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)ptr;
  TH2F *hist;

  Int_t out_count=0,tot_count=0;
  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      printf("doing row %d\n",i);
      for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
	{
	  
	  AliL3DigitData *bins = tempPt->fDigitData;
	  Float_t xyz[3];
	  Int_t sector,row;
	  fTransform->Slice2Sector(fSlice,i,sector,row);
	  fTransform->Raw2Local(xyz,sector,row,bins[dig].fPad,bins[dig].fTime);
	  Double_t eta = fTransform->GetEta(xyz);
	  eta_index = (Int_t)(eta/etaslice);
	  if(eta_index < 0 || eta_index >= n_eta_segments)
	    {
	      //printf("Eta index out of range %d\n",eta_index);
	      out_count++;
	      continue;
	    }
	  tot_count++;
	  hist = histos[eta_index];
	  if(!hist)
	    {
	      printf("Error getting histogramm!\n");
	      return;
	    }
	  
	  Int_t ymin = hist->GetYaxis()->GetFirst();
	  Int_t ymax = hist->GetYaxis()->GetLast();
	  
	  Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
	  Double_t phi_pix = fTransform->GetPhi(xyz);
	  for(Int_t p=ymin; p<=ymax; p++)
	    {
	      Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	      Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	      //printf("kappa %f phi0 %f\n",kappa,phi0);
	      hist->Fill(kappa,phi0,bins[dig].fCharge);
	    }
	}
      
      
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
      
    }
  
  printf("There were %d pixels out of range and %d inside\n", out_count,tot_count);

}


void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
{

  for(Int_t i=0; i<tracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
     
      Double_t xmin = fTransform->Row2X(track->GetFirstRow());
      Double_t xmax = fTransform->Row2X(track->GetLastRow());
      
      Double_t a = -1./tan(track->GetPsiLine());
      Double_t b = track->GetDLine()/sin(track->GetPsiLine());
      
      Double_t ymin = a*xmin + b;
      Double_t ymax = a*xmax + b;
      
      Double_t middle_x = xmin + (xmax-xmin)/2;
      Double_t middle_y = ymin + (ymax-ymin)/2;
      
      Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
      Double_t phi = atan2(middle_y,middle_x);
      Double_t phi0 = 2*phi - track->GetPsiLine();
      Double_t kappa = 2*sin(phi-phi0)/r_middle;
      hist->Fill(kappa,phi0,track->GetWeight());
     
    }

  
}

void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
{
  //Transform every pixel into histogram, using parametrisation:
  //D = x*cos(psi) + y*sin(psi)

  //  printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);

  AliL3Digits *pix1;
  
  Int_t xmin = hist->GetXaxis()->GetFirst();
  Int_t xmax = hist->GetXaxis()->GetLast();
  
  Double_t x0 = fTransform->Row2X(ref_row);
  Double_t y0 = 0;

  Int_t sector,row;
  //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
  
  Double_t phi_min = -10*ToRad;
  Double_t phi_max = 10*ToRad;
  Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
  

  Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
  Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
    
  
  Int_t index;
  
  for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
    {
      for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
	{
	  
	  index = phi*fNumOfPadRows + padrow;
	  //printf("Looping index %d\n",index);
	  if(index > fContainerBounds || index < 0)
	    {
	      printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
	      return;
	    }
	  for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
	  //for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
	    {
	      //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
	      Float_t xyz[3];
	      fTransform->Slice2Sector(fSlice,padrow,sector,row);
	      fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
	      	  
    	      if(raw)
		raw->Fill(xyz[0],xyz[1],pix1->fCharge);

	      xyz[0] = xyz[0]-x0;
	      xyz[1] = xyz[1]-y0;
	      
	      //printf("calculating...");
	      for(Int_t d=xmin+1; d<=xmax; d++)
		{
		  Double_t psi = hist->GetXaxis()->GetBinCenter(d);
		  Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
		  
		  Short_t signal = pix1->fCharge;
		  hist->Fill(psi,D,signal);
		  //printf("Filling histo, psi %f D %f\n",psi,D);
		}
	      //printf("done\n");
	    }
	  //printf(" done\n");
	}
      
    }
    
}


void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
{

  //read data from rootfile. more or less obsolete code this.

  TFile *file = new TFile(rootfile);
  file->cd();

  AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
  TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");      
    
  AliSimDigits da, *digarr=&da;
  t->GetBranch("Segment")->SetAddress(&digarr);
  Stat_t num_of_entries=t->GetEntries();

  Int_t digit_counter=0;
  Float_t xyz[3];
  Double_t eta;
  
  Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
  printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
  
  if(fRowContainer)
    delete [] fRowContainer;
  fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
  memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));


  //fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
  //printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);

  
  //  if(fPhiRowContainer)
  //  delete [] fPhiRowContainer;
  //  fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
  //  memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
  
  fContainerBounds = (fNumEtaSegments+1)*(fNumOfPadRows+1);
  if(fVolume)
    delete [] fVolume;
  fVolume = new AliL3HoughContainer[fContainerBounds];
  memset(fVolume,0,fContainerBounds*sizeof(AliL3HoughContainer));
  
  //  Double_t phi_min = -10*ToRad;
  //  Double_t phi_max = 10*ToRad;
  //  Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
  
  Double_t eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
  Int_t index;
  digit_counter=0;

  //printf("\nLoading ALL pixels in slice\n\n");

  for (Int_t i=0; i<num_of_entries; i++) 
    { 
      t->GetEvent(i);
      Int_t sector; 
      Int_t row;    
      param->AdjustSectorRow(digarr->GetID(),sector,row);
      Int_t slice,padrow;
      fTransform->Sector2Slice(slice,padrow,sector,row);
      if(slice != fSlice) continue;
      if(padrow < NRows[fPatch][0]) continue;
      if(padrow > NRows[fPatch][1]) break;
      digarr->First();
      do {
	Int_t time=digarr->CurrentRow();
	Int_t pad=digarr->CurrentColumn();
	Short_t signal=digarr->CurrentDigit();
	if(time < param->GetMaxTBin()-1 && time > 0)
	  if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
	     && digarr->GetDigit(time+1,pad) <= param->GetZeroSup()) 
	    continue;
	
	
	fTransform->Raw2Global(xyz,sector,row,pad,time);
	eta = fTransform->GetEta(xyz);
	
	if(eta < fEtaMin || eta > fEtaMax) continue;
	fTransform->Global2Local(xyz,sector);
	
	
	//phi = fTransform->GetPhi(xyz);
	if(hist)
	  hist->Fill(xyz[0],xyz[1],signal);
	
	AliL3Digits *dig = new AliL3Digits;
	dig->fIndex = digit_counter;
	digit_counter++;
	dig->fCharge = signal;
	dig->fPad = pad;
	dig->fTime = time;
	
	if(fRowContainer[padrow].first == NULL)
	  fRowContainer[padrow].first = (void*)dig;
	else
	  ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
	fRowContainer[padrow].last = (void*)dig;
	
	//thisHit->etaIndex=(Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1);
	Int_t eta_index = (Int_t)((eta-fEtaMin)/eta_slice);
	index = eta_index*fNumOfPadRows + padrow;
	if(index > fContainerBounds || index < 0)
	  {
	    
	    //printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
	    continue;
	  }
	
	if(fVolume[index].first == NULL)
	  fVolume[index].first = (void*)dig;
	else
	  ((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
	fVolume[index].last = (void*)dig;
	
	
	//  Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
	//  index = phi_index*fNumOfPadRows + padrow;
	//  if(phi_index > fContainerBounds || phi_index < 0)
	// {
	//  printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
	//  continue;
	//  }
	  
	//  if(fPhiRowContainer[index].first == NULL)
	//  fPhiRowContainer[index].first = (void*)dig;
	//  else
	//  ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
	//  fPhiRowContainer[index].last=(void*)dig;
	
      }while (digarr->Next());
      
    }
  
  fNDigits = digit_counter;
  printf("digitcounter %d\n",digit_counter);
  printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
  file->Close();
  delete file;
  
}
*/
Commit	Line	Data
4cafa5fc	1	//Author: Anders Strand Vestbo
	2	//Last Modified: 28.6.01
	3
4de874d1	4	#include <math.h>
4cafa5fc	5	#include <TH2.h>
9f33a1db	6	#include <time.h>
4de874d1	7
	8	#include <TFile.h>
	9	#include <TTree.h>
4de874d1	10
	11	#include "AliTPCParam.h"
	12	#include "AliSimDigits.h"
	13
4cafa5fc	14	#include "AliL3Histogram.h"
4cafa5fc	15	#include "AliL3Logging.h"
4de874d1	16	#include "AliL3Defs.h"
4de874d1	17	#include "AliL3Transform.h"
4de874d1	18	#include "AliL3HoughTransformer.h"
	19	#include "AliL3HoughTrack.h"
	20	#include "AliL3TrackArray.h"
4cafa5fc	21	#include "AliL3DigitData.h"
4de874d1	22
	23	ClassImp(AliL3HoughTransformer)
	24
	25	AliL3HoughTransformer::AliL3HoughTransformer()
	26	{
	27	//Default constructor
4cafa5fc	28
4de874d1	29	}
	30
	31
52a2a604	32	AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
4de874d1	33	{
	34	//Constructor
	35
	36	fTransform = new AliL3Transform();
4de874d1	37
	38	fEtaMin = etarange[0];
	39	fEtaMax = etarange[1];
	40	fSlice = slice;
	41	fPatch = patch;
4de874d1	42	fNumOfPadRows=NRowsSlice;
	43	}
	44
52a2a604	45	AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
	46	{
	47
	48	fTransform = new AliL3Transform();
	49	if(etarange)
	50	{
	51	//assumes you want to look at a given etaslice only
	52	fEtaMin = etarange[0];
	53	fEtaMax = etarange[1];
	54	}
	55	else
	56	{
	57	//transform in all etaslices
	58	fEtaMin = 0;
	59	fEtaMax = slice < 18 ? 0.9 : -0.9;
	60	}
	61	fNumEtaSegments = n_eta_segments;
	62	fSlice = slice;
	63	fPatch = patch;
	64	fNumOfPadRows = NRowsSlice;
4cafa5fc	65
	66	fNRowsInPatch = NRows[fPatch][1]-NRows[fPatch][0] + 1;
	67	fBinTableBounds = (fNRowsInPatch+1)*(MaxNPads+1);
	68	fNDigitRowData = 0;
	69	fDigitRowData = 0;
	70	fHistoPt = 0;
52a2a604	71	}
52a2a604	72
4de874d1	73
	74	AliL3HoughTransformer::~AliL3HoughTransformer()
	75	{
	76	//Destructor
4cafa5fc	77	if(fBinTable)
	78	{
	79	for(Int_t i=0; i<fBinTableBounds; i++)
	80	delete [] fBinTable[i];
	81	delete [] fBinTable;
	82	}
	83	if(fEtaIndex)
	84	delete [] fEtaIndex;
	85	if(fTrackTable)
	86	{
	87	for(Int_t i=0; i<fNumEtaSegments; i++)
	88	delete [] fTrackTable[i];
	89	delete [] fTrackTable;
	90	}
4de874d1	91	if(fTransform)
4de874d1	92	delete fTransform;
4cafa5fc	93
	94
	95	}
	96
	97	void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
	98	{
	99	fNDigitRowData = ndigits;
	100	fDigitRowData = ptr;
	101	}
	102
	103	void AliL3HoughTransformer::InitTables()
	104	{
	105	//Create LUT for the circle transform.
	106	//the actual transformation is done in TransformTables.
	107
	108	AliL3Histogram *hist = fHistoPt;
	109
	110	Int_t nbinsy = hist->GetNbinsY();
	111
	112	fBinTable = new Int_t*[fBinTableBounds];
	113	Int_t index;
	114
	115	Double_t etaslice = fEtaMax/fNumEtaSegments;
	116
	117	Int_t etabounds = (MaxNPads+1)(fNRowsInPatch+1)(MaxNTimeBins+1);
	118
	119	fEtaIndex = new Int_t[etabounds];
	120	for(Int_t i=0; i<etabounds; i++)
	121	fEtaIndex[i] = -1;
	122
9f33a1db	123	fTrackTable = new Char_t*[fNumEtaSegments];
4cafa5fc	124	for(Int_t i=0; i<fNumEtaSegments; i++)
4cafa5fc	125	{
9f33a1db	126	fTrackTable[i] = new Char_t[fBinTableBounds];
4cafa5fc	127	for(Int_t j=0; j<fBinTableBounds; j++)
	128	fTrackTable[i][j] = 0;
	129	}
	130
	131	for(Int_t r=NRows[fPatch][0]; r<=NRows[fPatch][1]; r++)
	132	{
	133	Int_t prow = r - NRows[fPatch][0];
	134
	135	for(Int_t p=0; p<fTransform->GetNPads(r); p++)
	136	{
	137	Float_t xyz[3];
	138	Int_t sector,row;
	139	for(Int_t t=0; t<fTransform->GetNTimeBins(); t++)
	140	{
	141	fTransform->Slice2Sector(fSlice,r,sector,row);
	142	fTransform->Raw2Local(xyz,sector,row,p,t);
	143	Double_t eta = fTransform->GetEta(xyz);
	144	if(eta < fEtaMin \|\| eta > fEtaMax) continue;
	145	Int_t ind = (prow<<17) + (p<<9) + t;
	146	if(fEtaIndex[ind]>=0)
	147	printf("AliL3HoughTransformer::InitTables : Overlapping indexes in eta!!\n");
	148	Int_t eta_index = (Int_t)(eta/etaslice);
	149	if(eta_index < 0 \|\| eta_index > fNumEtaSegments)
	150	continue;
	151	fEtaIndex[ind] = eta_index;
	152
	153	}
	154
	155	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	156	Double_t phi_pix = fTransform->GetPhi(xyz);
	157	index = (prow<<8) + p;
	158	fBinTable[index] = new Int_t[nbinsy+1];
	159
	160	for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
	161	{
	162	Double_t phi0 = hist->GetBinCenterY(b);
	163	Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	164	Int_t bin = hist->FindBin(kappa,phi0);
	165	if(fBinTable[index][b]!=0)
	166	printf("AliL3HoughTransformer::InitTables : Overlapping indexes %d %d b %d\n",fBinTable[index][b],index,b);
	167	fBinTable[index][b] = bin;
	168	}
	169	}
	170
	171	}
	172
	173	}
	174
	175	void AliL3HoughTransformer::TransformTables(AliL3Histogram histos,AliL3Histogram images)
	176	{
	177	//Transform all the pixels while reading, and fill the corresponding histograms.
	178	//Transform is done using LUT created in InitTables.
	179	//fTrackTable : table telling whether a specific pixel is active (nonzero):
	180	//fTrackTable = 0 -> no track
9f33a1db	181	//fTrackindex > 0 -> track present
9f33a1db	182	//fTrackindex = -1 -> track has been removed (already found)
4cafa5fc	183	//fEtaIndex : table storing the etaindex -> used to find correct histogram to fill
	184	//fBinTable : table storing all the bins to fill for each nonzero pixel
	185
	186	Int_t eta_index;
	187	AliL3DigitRowData tempPt = (AliL3DigitRowData)fDigitRowData;
	188	AliL3Histogram *hist;
	189
	190	if(!tempPt)
4de874d1	191	{
4cafa5fc	192	LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformTables","data")<<
	193	"Zero datapointer"<<ENDLOG;
	194	return;
4de874d1	195	}
4cafa5fc	196
	197	Int_t out_count=0,tot_count=0;
	198	Int_t index,ind;
9f33a1db	199
4cafa5fc	200	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
	201	{
	202	Int_t prow = i - NRows[fPatch][0];
	203	AliL3DigitData *bins = tempPt->fDigitData;
	204	for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
	205	{
	206	index = (prow<<8) + bins[dig].fPad;
	207	ind = (prow<<17) + (bins[dig].fPad<<9) + bins[dig].fTime;
	208	eta_index = fEtaIndex[ind];
	209	if(eta_index < 0) continue; //pixel out of etarange
	210
9f33a1db	211	if(fTrackTable[eta_index][index]<0) continue; //this pixel has already been removed.
9f33a1db	212	fTrackTable[eta_index][index]++; //this pixel is now marked as active (it is nonzero)
4cafa5fc	213
	214	tot_count++;
	215	hist = histos[eta_index];
	216
	217	if(images)
	218	{
	219	//Display the transformed images.
	220	AliL3Histogram *image = images[eta_index];
	221	Float_t xyz_local[3];
	222	Int_t sector,row;
	223	fTransform->Slice2Sector(fSlice,i,sector,row);
	224	fTransform->Raw2Local(xyz_local,sector,row,bins[dig].fPad,bins[dig].fTime);
	225	image->Fill(xyz_local[0],xyz_local[1],bins[dig].fCharge);
	226	}
	227
	228	if(!hist)
	229	{
	230	printf("Error getting histogram!\n");
	231	continue;
	232	}
	233	for(Int_t p=hist->GetFirstYbin(); p<=hist->GetLastYbin(); p++)
	234	hist->AddBinContent(fBinTable[index][p],bins[dig].fCharge);
	235
	236	}
	237
	238	Byte_t tmp = (Byte_t)tempPt;
	239	Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
	240	tmp += size;
	241	tempPt = (AliL3DigitRowData*)tmp;
	242	}
	243
	244	}
	245
	246	void AliL3HoughTransformer::WriteTables()
	247	{
	248	//Write the tables to asci files.
	249
	250	AliL3Histogram *hist = fHistoPt;
	251	Char_t name[100];
	252	sprintf(name,"histogram_table_%d.txt",fPatch);
	253	FILE *file = fopen(name,"w");
	254
	255	Int_t etabounds = (MaxNPads+1)(fNRowsInPatch+1)(MaxNTimeBins+1);
	256	for(Int_t i=0; i<etabounds; i++)
	257	{
	258	if(fEtaIndex[i]<0) continue;
	259	fprintf(file,"%d %d\n",i,fEtaIndex[i]);
	260	}
	261	fclose(file);
	262
	263	sprintf(name,"bin_table_%d.txt",fPatch);
	264	FILE *file2 = fopen(name,"w");
	265	for(Int_t i=0; i<fBinTableBounds; i++)
	266	{
	267	if(!fBinTable[i]) continue;
	268	fprintf(file2,"%d ",i);
	269	for(Int_t j=hist->GetFirstYbin(); j<=hist->GetLastYbin(); j++)
	270	{
	271	fprintf(file2,"%d ",fBinTable[i][j]);
	272	}
	273	fprintf(file2,"\n");
	274	}
	275	fclose(file2);
9f33a1db	276	}
	277
	278	Double_t AliL3HoughTransformer::CpuTime()
	279	{
	280	return (Double_t)(clock()) / CLOCKS_PER_SEC;
4de874d1	281	}
4de874d1	282
4cafa5fc	283	/*
4de874d1	284	void AliL3HoughTransformer::InitTemplates(TH2F *hist)
	285	{
	286
	287	AliL3Digits *pixel;
	288
	289	Int_t ymin = hist->GetYaxis()->GetFirst();
	290	Int_t ymax = hist->GetYaxis()->GetLast();
	291	Int_t nbinsy = hist->GetNbinsY();
	292
	293	fIndex = new Int_t*[fNDigits];
	294	for(Int_t i=0; i<fNDigits; i++)
	295	fIndex[i] = new Int_t[nbinsy+1];
	296
	297	Int_t sector,row;
52a2a604	298
4de874d1	299	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
52a2a604	300	{
4de874d1	301
	302	for(pixel=(AliL3Digits)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits)pixel->nextRowPixel)
	303	{
	304	Float_t xyz[3];
	305	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	306	fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
	307
	308	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	309
	310	Double_t phi_pix = fTransform->GetPhi(xyz);
4de874d1	311	Int_t index = pixel->fIndex;
	312	if(index >= fNDigits)
	313	printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
	314	for(Int_t p=ymin; p<=ymax; p++)
	315	{
	316
	317	Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	318	Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
52a2a604	319	//printf("kappa %f phi0 %f\n",kappa,phi0);
4de874d1	320	Int_t bin = hist->FindBin(kappa,phi0);
	321	if(fIndex[index][p]!=0)
	322	printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
	323	fIndex[index][p] = bin;
	324	}
	325	}
	326	}
	327
	328	}
	329
	330
52a2a604	331	void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
4de874d1	332	{
	333	//Transformation is done with respect to local coordinates in slice.
	334	//Transform every pixel into whole phirange, using parametrisation:
	335	//kappa = 2*sin(phi-phi0)/R
4cafa5fc	336	//Assumes you run InitTemplates first!!!!
4de874d1	337
4de874d1	338	AliL3Digits *pix1;
e4c21048	339
4de874d1	340	Int_t nbinsy = hist->GetNbinsY();
4de874d1	341
52a2a604	342	Int_t totsignal=0,vol_index;
	343	if(fNumEtaSegments==1)
	344	eta_index=0; //only looking in one etaslice.
	345
4de874d1	346	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
4de874d1	347	{
52a2a604	348	vol_index = eta_index*fNumOfPadRows + padrow;
	349	//for(pix1=(AliL3Digits)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits)pix1->nextRowPixel)
	350	for(pix1=(AliL3Digits)fVolume[vol_index].first; pix1!=0; pix1=(AliL3Digits)pix1->fNextVolumePixel)
4de874d1	351	{
4cafa5fc	352	Float_t xyz[3];
	353	fTransform->Raw2Global(xyz,2,padrow,pix1->fPad,pix1->fTime);
	354	Double_t eta = fTransform->GetEta(xyz);
	355	if(eta < fEtaMin \|\| eta > fEtaMax)
	356	printf("\n Eta OUT OF RANGE\n");
	357
4de874d1	358	Short_t signal = pix1->fCharge;
4de874d1	359	Int_t index = pix1->fIndex;
52a2a604	360	if(index < 0) continue; //This pixel has been removed.
52a2a604	361	totsignal += signal;
4de874d1	362	for(Int_t p=0; p <= nbinsy; p++)
	363	hist->AddBinContent(fIndex[index][p],signal);
	364
	365	}
	366
	367	}
	368
52a2a604	369	printf("Total signal %d\n",totsignal);
4de874d1	370	}
4de874d1	371
4cafa5fc	372	void AliL3HoughTransformer::Transform2Circle(TH2F *histos,Int_t n_eta_segments,UInt_t ndigits,AliL3DigitRowData ptr)
	373	{
	374	//Transform all the pixels while reading them, and fill the corresponding histograms.
	375	//Everything is done in one go here.
	376
	377	Double_t etaslice = 0.9/n_eta_segments;
	378	Int_t eta_index;
	379	AliL3DigitRowData tempPt = (AliL3DigitRowData)ptr;
	380	TH2F *hist;
	381
	382	Int_t out_count=0,tot_count=0;
	383	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
	384	{
	385	printf("doing row %d\n",i);
	386	for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
	387	{
	388
	389	AliL3DigitData *bins = tempPt->fDigitData;
	390	Float_t xyz[3];
	391	Int_t sector,row;
	392	fTransform->Slice2Sector(fSlice,i,sector,row);
	393	fTransform->Raw2Local(xyz,sector,row,bins[dig].fPad,bins[dig].fTime);
	394	Double_t eta = fTransform->GetEta(xyz);
	395	eta_index = (Int_t)(eta/etaslice);
	396	if(eta_index < 0 \|\| eta_index >= n_eta_segments)
	397	{
	398	//printf("Eta index out of range %d\n",eta_index);
	399	out_count++;
	400	continue;
	401	}
	402	tot_count++;
	403	hist = histos[eta_index];
	404	if(!hist)
	405	{
	406	printf("Error getting histogramm!\n");
	407	return;
	408	}
	409
	410	Int_t ymin = hist->GetYaxis()->GetFirst();
	411	Int_t ymax = hist->GetYaxis()->GetLast();
	412
	413	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	414	Double_t phi_pix = fTransform->GetPhi(xyz);
	415	for(Int_t p=ymin; p<=ymax; p++)
	416	{
	417	Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	418	Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	419	//printf("kappa %f phi0 %f\n",kappa,phi0);
	420	hist->Fill(kappa,phi0,bins[dig].fCharge);
	421	}
	422	}
	423
	424
	425	Byte_t tmp = (Byte_t)tempPt;
	426	Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
	427	tmp += size;
	428	tempPt = (AliL3DigitRowData*)tmp;
	429
	430	}
4de874d1	431
4cafa5fc	432	printf("There were %d pixels out of range and %d inside\n", out_count,tot_count);
	433
	434	}
	435
4de874d1	436
	437	void AliL3HoughTransformer::TransformLines2Circle(TH2F hist,AliL3TrackArray tracks)
	438	{
	439
	440	for(Int_t i=0; i<tracks->GetNTracks(); i++)
	441	{
	442	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	443	if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
	444
	445	Double_t xmin = fTransform->Row2X(track->GetFirstRow());
	446	Double_t xmax = fTransform->Row2X(track->GetLastRow());
	447
	448	Double_t a = -1./tan(track->GetPsiLine());
	449	Double_t b = track->GetDLine()/sin(track->GetPsiLine());
	450
	451	Double_t ymin = a*xmin + b;
	452	Double_t ymax = a*xmax + b;
	453
	454	Double_t middle_x = xmin + (xmax-xmin)/2;
	455	Double_t middle_y = ymin + (ymax-ymin)/2;
	456
	457	Double_t r_middle = sqrt(middle_xmiddle_x + middle_ymiddle_y);
	458	Double_t phi = atan2(middle_y,middle_x);
	459	Double_t phi0 = 2*phi - track->GetPsiLine();
	460	Double_t kappa = 2*sin(phi-phi0)/r_middle;
	461	hist->Fill(kappa,phi0,track->GetWeight());
	462
	463	}
	464
	465
	466	}
	467
	468	void AliL3HoughTransformer::Transform2Line(TH2F hist,Int_t ref_row,Int_t rowrange,Double_t phirange,TH2F raw)
	469	{
	470	//Transform every pixel into histogram, using parametrisation:
	471	//D = xcos(psi) + ysin(psi)
	472
	473	// printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
	474
	475	AliL3Digits *pix1;
	476
	477	Int_t xmin = hist->GetXaxis()->GetFirst();
	478	Int_t xmax = hist->GetXaxis()->GetLast();
	479
	480	Double_t x0 = fTransform->Row2X(ref_row);
	481	Double_t y0 = 0;
	482
	483	Int_t sector,row;
	484	//for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	485
	486	Double_t phi_min = -10*ToRad;
	487	Double_t phi_max = 10*ToRad;
	488	Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
	489
	490
	491	Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
	492	Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
	493
	494
	495	Int_t index;
	496
	497	for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
	498	{
	499	for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
500	{
501
502	index = phi*fNumOfPadRows + padrow;
503	//printf("Looping index %d\n",index);
504	if(index > fContainerBounds \|\| index < 0)
505	{
506	printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
507	return;
508	}
52a2a604	509	for(pix1=(AliL3Digits)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits)pix1->nextRowPixel)
52a2a604	510	//for(pix1=(AliL3Digits)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits)pix1->nextPhiRowPixel)
4de874d1	511	{
	512	//printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
	513	Float_t xyz[3];
	514	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	515	fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
	516
	517	if(raw)
	518	raw->Fill(xyz[0],xyz[1],pix1->fCharge);
	519
	520	xyz[0] = xyz[0]-x0;
	521	xyz[1] = xyz[1]-y0;
	522
	523	//printf("calculating...");
	524	for(Int_t d=xmin+1; d<=xmax; d++)
	525	{
	526	Double_t psi = hist->GetXaxis()->GetBinCenter(d);
	527	Double_t D = xyz[0]cos(psi) + xyz[1]sin(psi);
	528
	529	Short_t signal = pix1->fCharge;
	530	hist->Fill(psi,D,signal);
	531	//printf("Filling histo, psi %f D %f\n",psi,D);
	532	}
	533	//printf("done\n");
	534	}
	535	//printf(" done\n");
	536	}
	537
	538	}
	539
	540	}
	541
4cafa5fc	542
4de874d1	543	void AliL3HoughTransformer::GetPixels(Char_t rootfile,TH2F hist)
	544	{
	545
4cafa5fc	546	//read data from rootfile. more or less obsolete code this.
4cafa5fc	547
4de874d1	548	TFile *file = new TFile(rootfile);
	549	file->cd();
	550
	551	AliTPCParam param=(AliTPCParam )file->Get("75x40_100x60");
	552	TTree t=(TTree)file->Get("TreeD_75x40_100x60");
	553
	554	AliSimDigits da, *digarr=&da;
	555	t->GetBranch("Segment")->SetAddress(&digarr);
	556	Stat_t num_of_entries=t->GetEntries();
	557
	558	Int_t digit_counter=0;
	559	Float_t xyz[3];
52a2a604	560	Double_t eta;
4de874d1	561
	562	Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
	563	printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
	564
	565	if(fRowContainer)
	566	delete [] fRowContainer;
	567	fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
	568	memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
	569
	570
52a2a604	571	//fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
52a2a604	572	//printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
4de874d1	573
4cafa5fc	574
	575	// if(fPhiRowContainer)
	576	// delete [] fPhiRowContainer;
	577	// fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
	578	// memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
	579
52a2a604	580	fContainerBounds = (fNumEtaSegments+1)*(fNumOfPadRows+1);
	581	if(fVolume)
	582	delete [] fVolume;
	583	fVolume = new AliL3HoughContainer[fContainerBounds];
	584	memset(fVolume,0,fContainerBounds*sizeof(AliL3HoughContainer));
4cafa5fc	585
	586	// Double_t phi_min = -10*ToRad;
	587	// Double_t phi_max = 10*ToRad;
	588	// Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
	589
52a2a604	590	Double_t eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
4de874d1	591	Int_t index;
	592	digit_counter=0;
	593
4cafa5fc	594	//printf("\nLoading ALL pixels in slice\n\n");
52a2a604	595
4de874d1	596	for (Int_t i=0; i<num_of_entries; i++)
	597	{
	598	t->GetEvent(i);
	599	Int_t sector;
	600	Int_t row;
	601	param->AdjustSectorRow(digarr->GetID(),sector,row);
	602	Int_t slice,padrow;
	603	fTransform->Sector2Slice(slice,padrow,sector,row);
	604	if(slice != fSlice) continue;
	605	if(padrow < NRows[fPatch][0]) continue;
	606	if(padrow > NRows[fPatch][1]) break;
	607	digarr->First();
	608	do {
	609	Int_t time=digarr->CurrentRow();
	610	Int_t pad=digarr->CurrentColumn();
	611	Short_t signal=digarr->CurrentDigit();
	612	if(time < param->GetMaxTBin()-1 && time > 0)
	613	if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
	614	&& digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
	615	continue;
	616
	617
	618	fTransform->Raw2Global(xyz,sector,row,pad,time);
	619	eta = fTransform->GetEta(xyz);
52a2a604	620
4de874d1	621	if(eta < fEtaMin \|\| eta > fEtaMax) continue;
	622	fTransform->Global2Local(xyz,sector);
	623
	624
52a2a604	625	//phi = fTransform->GetPhi(xyz);
4de874d1	626	if(hist)
	627	hist->Fill(xyz[0],xyz[1],signal);
	628
	629	AliL3Digits *dig = new AliL3Digits;
	630	dig->fIndex = digit_counter;
	631	digit_counter++;
	632	dig->fCharge = signal;
	633	dig->fPad = pad;
	634	dig->fTime = time;
	635
	636	if(fRowContainer[padrow].first == NULL)
	637	fRowContainer[padrow].first = (void*)dig;
	638	else
	639	((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
	640	fRowContainer[padrow].last = (void*)dig;
	641
52a2a604	642	//thisHit->etaIndex=(Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1);
	643	Int_t eta_index = (Int_t)((eta-fEtaMin)/eta_slice);
	644	index = eta_index*fNumOfPadRows + padrow;
	645	if(index > fContainerBounds \|\| index < 0)
4de874d1	646	{
52a2a604	647
52a2a604	648	//printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
4de874d1	649	continue;
4de874d1	650	}
52a2a604	651
	652	if(fVolume[index].first == NULL)
	653	fVolume[index].first = (void*)dig;
4de874d1	654	else
52a2a604	655	((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
52a2a604	656	fVolume[index].last = (void*)dig;
4de874d1	657
4cafa5fc	658
	659	// Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
	660	// index = phi_index*fNumOfPadRows + padrow;
	661	// if(phi_index > fContainerBounds \|\| phi_index < 0)
	662	// {
	663	// printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
	664	// continue;
	665	// }
52a2a604	666
4cafa5fc	667	// if(fPhiRowContainer[index].first == NULL)
	668	// fPhiRowContainer[index].first = (void*)dig;
	669	// else
	670	// ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
	671	// fPhiRowContainer[index].last=(void*)dig;
	672
4de874d1	673	}while (digarr->Next());
	674
	675	}
	676
	677	fNDigits = digit_counter;
	678	printf("digitcounter %d\n",digit_counter);
	679	printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
	680	file->Close();
	681	delete file;
	682
	683	}
4cafa5fc	684	*/