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

#include <math.h>

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

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

#include "AliL3Defs.h"
#include "AliL3Transform.h"
#include "AliL3HoughPixel.h"
#include "AliL3HoughTransformer.h"
#include "AliL3HoughTrack.h"
#include "AliL3TrackArray.h"

ClassImp(AliL3HoughTransformer)

AliL3HoughTransformer::AliL3HoughTransformer()
{
  //Default constructor
  fTransform = 0;
  fEtaMin = 0;
  fEtaMax = 0;
  fSlice = 0;
  fPatch = 0;
  fRowContainer = 0;
  fPhiRowContainer = 0;
  fVolume=0;
  fNumOfPadRows=0;
  fContainerBounds=0;
  fNDigits=0;
  fIndex = 0;
}


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;
}


AliL3HoughTransformer::~AliL3HoughTransformer()
{
  //Destructor
  if(fRowContainer)
    delete [] fRowContainer;
  if(fTransform)
    delete fTransform;
  if(fPhiRowContainer)
    delete [] fPhiRowContainer;
  if(fVolume)
    delete [] fVolume;
  if(fIndex)
    {
      for(Int_t i=0; i<fNDigits; i++)
	delete [] fIndex[i];
      delete [] fIndex;
    }
}

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::CountBins()
{
  
  Int_t middle_row = 87; //middle of the slice
  
  Double_t r_in_bundle = fTransform->Row2X(middle_row);
  //  Double_t phi_min = (fSlice*20 - 10)*ToRad;
  //Double_t phi_max = (fSlice*20 + 10)*ToRad;
  Double_t phi_min = -15*ToRad;
  Double_t phi_max = 15*ToRad;

  Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
  Double_t min_phi0 = 10000;
  Double_t max_phi0 = 0;
  Double_t min_kappa = 100000;
  Double_t max_kappa = 0;
  
  Int_t xbin = 60;
  Int_t ybin = 60;
  Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
  Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
  
  TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  
  for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
    {
      for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
	{
	  for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
	    {
	      Float_t xyz[3];
	      Int_t sector,row;
	      fTransform->Slice2Sector(fSlice,padrow,sector,row);
	      fTransform->Raw2Global(xyz,sector,row,pad,time);
	      Double_t eta = fTransform->GetEta(xyz);
	      if(eta < fEtaMin || eta > fEtaMax) continue;
	      fTransform->Global2Local(xyz,sector);
	      Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
	      Double_t phi_pix = fTransform->GetPhi(xyz);
	      
	      for(Int_t p=0; p<=fNPhiSegments; p++)
		{
		  Double_t phi_in_bundle = phi_min + p*phi_slice;
		  
		  Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
		  
		  Double_t phi0 = atan(tanPhi0);
		  //  if(phi0 < 0.55 || phi0 > 0.88) continue;
		  
		  //if(phi0 < 0) phi0 = phi0 +2*Pi;
		  //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
		  
		  Double_t angle = phi_pix - phi0;
		  Double_t kappa = 2*sin(angle)/r_pix;
		  histo->Fill(kappa,phi0,1);
		  if(phi0 < min_phi0)
		    min_phi0 = phi0;
		  if(phi0 > max_phi0)
		    max_phi0 = phi0;
		  if(kappa < min_kappa)
		    min_kappa = kappa;
		  if(kappa > max_kappa)
		    max_kappa = kappa;
		    		  
		}
	      
	    }
	  
	}
      
    }
  Int_t count=0,bi=0;
    
  Int_t xmin = histo->GetXaxis()->GetFirst();
  Int_t xmax = histo->GetXaxis()->GetLast();
  Int_t ymin = histo->GetYaxis()->GetFirst();
  Int_t ymax = histo->GetYaxis()->GetLast();


  for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
    {
      for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
	{
	  bi++;
	  Int_t bin = histo->GetBin(xbin,ybin);
	  if(histo->GetBinContent(bin)>0)
	    count++;
	}
    }


  printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
  printf("Phi, min %f max %f\n",min_phi0,max_phi0);
  printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
  histo->Draw("box");
}


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 firstly!!!!

  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)
	{
	  
	  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 *hist)
  {
  //Transformation is done with respect to local coordinates in slice.
  //Transform every pixel into whole phirange, using parametrisation:
  //kappa = 2*sin(phi-phi0)/R
  
  printf("Transforming 1 pixel only\n");
  
  AliL3Digits *pix1;
  Int_t sector,row;
  
  Int_t ymin = hist->GetYaxis()->GetFirst();
  Int_t ymax = hist->GetYaxis()->GetLast();
  
  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
  {
  
  for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
  {
  
  Float_t xyz[3];
  fTransform->Slice2Sector(fSlice,padrow,sector,row);
  fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
  
  Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  
  Double_t phi_pix = fTransform->GetPhi(xyz);
  Short_t signal = pix1->fCharge;
  
  for(Int_t p=ymin+1; p<=ymax; p++)
  {
  Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
  Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
  hist->Fill(kappa,phi0,signal);
  }
  
  }
  
  }
  
  
  }
*/

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)
{

  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
4de874d1	1	#include <math.h>
	2
	3	#include <TFile.h>
	4	#include <TTree.h>
	5	#include <TH2.h>
	6
	7	#include "AliTPCParam.h"
	8	#include "AliSimDigits.h"
	9
	10	#include "AliL3Defs.h"
	11	#include "AliL3Transform.h"
	12	#include "AliL3HoughPixel.h"
	13	#include "AliL3HoughTransformer.h"
	14	#include "AliL3HoughTrack.h"
	15	#include "AliL3TrackArray.h"
	16
	17	ClassImp(AliL3HoughTransformer)
	18
	19	AliL3HoughTransformer::AliL3HoughTransformer()
	20	{
	21	//Default constructor
	22	fTransform = 0;
	23	fEtaMin = 0;
	24	fEtaMax = 0;
	25	fSlice = 0;
	26	fPatch = 0;
	27	fRowContainer = 0;
	28	fPhiRowContainer = 0;
52a2a604	29	fVolume=0;
4de874d1	30	fNumOfPadRows=0;
	31	fContainerBounds=0;
	32	fNDigits=0;
	33	fIndex = 0;
	34	}
	35
	36
52a2a604	37	AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
4de874d1	38	{
	39	//Constructor
	40
	41	fTransform = new AliL3Transform();
4de874d1	42
	43	fEtaMin = etarange[0];
	44	fEtaMax = etarange[1];
	45	fSlice = slice;
	46	fPatch = patch;
4de874d1	47	fNumOfPadRows=NRowsSlice;
	48	}
	49
52a2a604	50	AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
	51	{
	52
	53	fTransform = new AliL3Transform();
	54	if(etarange)
	55	{
	56	//assumes you want to look at a given etaslice only
	57	fEtaMin = etarange[0];
	58	fEtaMax = etarange[1];
	59	}
	60	else
	61	{
	62	//transform in all etaslices
	63	fEtaMin = 0;
	64	fEtaMax = slice < 18 ? 0.9 : -0.9;
	65	}
	66	fNumEtaSegments = n_eta_segments;
	67	fSlice = slice;
	68	fPatch = patch;
	69	fNumOfPadRows = NRowsSlice;
	70	}
	71
4de874d1	72
	73	AliL3HoughTransformer::~AliL3HoughTransformer()
	74	{
	75	//Destructor
	76	if(fRowContainer)
	77	delete [] fRowContainer;
	78	if(fTransform)
	79	delete fTransform;
	80	if(fPhiRowContainer)
	81	delete [] fPhiRowContainer;
52a2a604	82	if(fVolume)
52a2a604	83	delete [] fVolume;
4de874d1	84	if(fIndex)
	85	{
	86	for(Int_t i=0; i<fNDigits; i++)
	87	delete [] fIndex[i];
	88	delete [] fIndex;
	89	}
	90	}
	91
	92	void AliL3HoughTransformer::InitTemplates(TH2F *hist)
	93	{
	94
	95	AliL3Digits *pixel;
	96
	97	Int_t ymin = hist->GetYaxis()->GetFirst();
	98	Int_t ymax = hist->GetYaxis()->GetLast();
	99	Int_t nbinsy = hist->GetNbinsY();
	100
	101	fIndex = new Int_t*[fNDigits];
	102	for(Int_t i=0; i<fNDigits; i++)
	103	fIndex[i] = new Int_t[nbinsy+1];
	104
	105	Int_t sector,row;
52a2a604	106
4de874d1	107	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
52a2a604	108	{
4de874d1	109
	110	for(pixel=(AliL3Digits)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits)pixel->nextRowPixel)
	111	{
	112	Float_t xyz[3];
	113	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	114	fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
	115
	116	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	117
	118	Double_t phi_pix = fTransform->GetPhi(xyz);
4de874d1	119	Int_t index = pixel->fIndex;
	120	if(index >= fNDigits)
	121	printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
	122	for(Int_t p=ymin; p<=ymax; p++)
	123	{
	124
	125	Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	126	Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
52a2a604	127	//printf("kappa %f phi0 %f\n",kappa,phi0);
4de874d1	128	Int_t bin = hist->FindBin(kappa,phi0);
	129	if(fIndex[index][p]!=0)
	130	printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
	131	fIndex[index][p] = bin;
	132	}
	133	}
	134	}
	135
	136	}
	137
	138
	139	void AliL3HoughTransformer::CountBins()
	140	{
	141
	142	Int_t middle_row = 87; //middle of the slice
	143
	144	Double_t r_in_bundle = fTransform->Row2X(middle_row);
	145	// Double_t phi_min = (fSlice20 - 10)ToRad;
	146	//Double_t phi_max = (fSlice20 + 10)ToRad;
	147	Double_t phi_min = -15*ToRad;
	148	Double_t phi_max = 15*ToRad;
	149
	150	Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
	151	Double_t min_phi0 = 10000;
	152	Double_t max_phi0 = 0;
	153	Double_t min_kappa = 100000;
	154	Double_t max_kappa = 0;
	155
	156	Int_t xbin = 60;
	157	Int_t ybin = 60;
	158	Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
	159	Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
	160
	161	TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
	162
	163	for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	164	{
	165	for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
	166	{
	167	for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
	168	{
	169	Float_t xyz[3];
	170	Int_t sector,row;
	171	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	172	fTransform->Raw2Global(xyz,sector,row,pad,time);
	173	Double_t eta = fTransform->GetEta(xyz);
	174	if(eta < fEtaMin \|\| eta > fEtaMax) continue;
	175	fTransform->Global2Local(xyz,sector);
	176	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	177	Double_t phi_pix = fTransform->GetPhi(xyz);
	178
	179	for(Int_t p=0; p<=fNPhiSegments; p++)
	180	{
	181	Double_t phi_in_bundle = phi_min + p*phi_slice;
	182
	183	Double_t tanPhi0 = (r_pixsin(phi_in_bundle)-r_in_bundlesin(phi_pix))/(r_pixcos(phi_in_bundle)-r_in_bundlecos(phi_pix));
	184
	185	Double_t phi0 = atan(tanPhi0);
	186	// if(phi0 < 0.55 \|\| phi0 > 0.88) continue;
	187
	188	//if(phi0 < 0) phi0 = phi0 +2*Pi;
	189	//Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
	190
	191	Double_t angle = phi_pix - phi0;
192	Double_t kappa = 2*sin(angle)/r_pix;
193	histo->Fill(kappa,phi0,1);
194	if(phi0 < min_phi0)
195	min_phi0 = phi0;
196	if(phi0 > max_phi0)
197	max_phi0 = phi0;
198	if(kappa < min_kappa)
199	min_kappa = kappa;
200	if(kappa > max_kappa)
201	max_kappa = kappa;
202
203	}
204
205	}
206
207	}
208
209	}
210	Int_t count=0,bi=0;
211
212	Int_t xmin = histo->GetXaxis()->GetFirst();
213	Int_t xmax = histo->GetXaxis()->GetLast();
214	Int_t ymin = histo->GetYaxis()->GetFirst();
215	Int_t ymax = histo->GetYaxis()->GetLast();
216
217
218	for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
219	{
220	for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
221	{
222	bi++;
223	Int_t bin = histo->GetBin(xbin,ybin);
224	if(histo->GetBinContent(bin)>0)
225	count++;
226	}
227	}
228
229
230	printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
231	printf("Phi, min %f max %f\n",min_phi0,max_phi0);
232	printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
233	histo->Draw("box");
234	}
235
236
52a2a604	237	void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
4de874d1	238	{
	239	//Transformation is done with respect to local coordinates in slice.
	240	//Transform every pixel into whole phirange, using parametrisation:
	241	//kappa = 2*sin(phi-phi0)/R
e4c21048	242	//Assumes you run InitTemplates firstly!!!!
4de874d1	243
4de874d1	244	AliL3Digits *pix1;
e4c21048	245
4de874d1	246	Int_t nbinsy = hist->GetNbinsY();
4de874d1	247
52a2a604	248	Int_t totsignal=0,vol_index;
	249	if(fNumEtaSegments==1)
	250	eta_index=0; //only looking in one etaslice.
	251
4de874d1	252	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
4de874d1	253	{
52a2a604	254	vol_index = eta_index*fNumOfPadRows + padrow;
	255	//for(pix1=(AliL3Digits)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits)pix1->nextRowPixel)
	256	for(pix1=(AliL3Digits)fVolume[vol_index].first; pix1!=0; pix1=(AliL3Digits)pix1->fNextVolumePixel)
4de874d1	257	{
	258
	259	Short_t signal = pix1->fCharge;
	260	Int_t index = pix1->fIndex;
52a2a604	261	if(index < 0) continue; //This pixel has been removed.
52a2a604	262	totsignal += signal;
4de874d1	263	for(Int_t p=0; p <= nbinsy; p++)
	264	hist->AddBinContent(fIndex[index][p],signal);
	265
	266	}
	267
	268	}
	269
52a2a604	270	printf("Total signal %d\n",totsignal);
4de874d1	271	}
	272
	273	/*
	274	void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
	275	{
	276	//Transformation is done with respect to local coordinates in slice.
	277	//Transform every pixel into whole phirange, using parametrisation:
	278	//kappa = 2*sin(phi-phi0)/R
	279
	280	printf("Transforming 1 pixel only\n");
	281
	282	AliL3Digits *pix1;
	283	Int_t sector,row;
	284
	285	Int_t ymin = hist->GetYaxis()->GetFirst();
	286	Int_t ymax = hist->GetYaxis()->GetLast();
	287
	288	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	289	{
	290
	291	for(pix1=(AliL3Digits)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits)pix1->nextRowPixel)
	292	{
	293
	294	Float_t xyz[3];
	295	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	296	fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
	297
	298	Double_t r_pix = sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	299
	300	Double_t phi_pix = fTransform->GetPhi(xyz);
	301	Short_t signal = pix1->fCharge;
	302
	303	for(Int_t p=ymin+1; p<=ymax; p++)
	304	{
	305	Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
	306	Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
	307	hist->Fill(kappa,phi0,signal);
	308	}
	309
	310	}
	311
	312	}
	313
	314
	315	}
	316	*/
	317
	318	void AliL3HoughTransformer::TransformLines2Circle(TH2F hist,AliL3TrackArray tracks)
	319	{
	320
	321	for(Int_t i=0; i<tracks->GetNTracks(); i++)
	322	{
	323	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	324	if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
	325
	326	Double_t xmin = fTransform->Row2X(track->GetFirstRow());
	327	Double_t xmax = fTransform->Row2X(track->GetLastRow());
	328
	329	Double_t a = -1./tan(track->GetPsiLine());
	330	Double_t b = track->GetDLine()/sin(track->GetPsiLine());
	331
	332	Double_t ymin = a*xmin + b;
	333	Double_t ymax = a*xmax + b;
	334
335	Double_t middle_x = xmin + (xmax-xmin)/2;
336	Double_t middle_y = ymin + (ymax-ymin)/2;
337
338	Double_t r_middle = sqrt(middle_xmiddle_x + middle_ymiddle_y);
339	Double_t phi = atan2(middle_y,middle_x);
340	Double_t phi0 = 2*phi - track->GetPsiLine();
341	Double_t kappa = 2*sin(phi-phi0)/r_middle;
342	hist->Fill(kappa,phi0,track->GetWeight());
343
344	}
345
346
347	}
348
349	void AliL3HoughTransformer::Transform2Line(TH2F hist,Int_t ref_row,Int_t rowrange,Double_t phirange,TH2F raw)
350	{
351	//Transform every pixel into histogram, using parametrisation:
352	//D = xcos(psi) + ysin(psi)
353
354	// printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
355
356	AliL3Digits *pix1;
357
358	Int_t xmin = hist->GetXaxis()->GetFirst();
359	Int_t xmax = hist->GetXaxis()->GetLast();
360
361	Double_t x0 = fTransform->Row2X(ref_row);
362	Double_t y0 = 0;
363
364	Int_t sector,row;
365	//for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
366
367	Double_t phi_min = -10*ToRad;
368	Double_t phi_max = 10*ToRad;
369	Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
370
371
372	Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
373	Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
374
375
376	Int_t index;
377
378	for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
379	{
380	for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
381	{
382
383	index = phi*fNumOfPadRows + padrow;
384	//printf("Looping index %d\n",index);
385	if(index > fContainerBounds \|\| index < 0)
386	{
387	printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
388	return;
389	}
52a2a604	390	for(pix1=(AliL3Digits)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits)pix1->nextRowPixel)
52a2a604	391	//for(pix1=(AliL3Digits)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits)pix1->nextPhiRowPixel)
4de874d1	392	{
	393	//printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
	394	Float_t xyz[3];
	395	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	396	fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
	397
	398	if(raw)
	399	raw->Fill(xyz[0],xyz[1],pix1->fCharge);
	400
	401	xyz[0] = xyz[0]-x0;
	402	xyz[1] = xyz[1]-y0;
	403
	404	//printf("calculating...");
	405	for(Int_t d=xmin+1; d<=xmax; d++)
	406	{
	407	Double_t psi = hist->GetXaxis()->GetBinCenter(d);
	408	Double_t D = xyz[0]cos(psi) + xyz[1]sin(psi);
	409
	410	Short_t signal = pix1->fCharge;
	411	hist->Fill(psi,D,signal);
	412	//printf("Filling histo, psi %f D %f\n",psi,D);
	413	}
	414	//printf("done\n");
	415	}
	416	//printf(" done\n");
	417	}
	418
	419	}
	420
	421	}
	422
	423	void AliL3HoughTransformer::GetPixels(Char_t rootfile,TH2F hist)
	424	{
	425
	426	TFile *file = new TFile(rootfile);
	427	file->cd();
	428
	429	AliTPCParam param=(AliTPCParam )file->Get("75x40_100x60");
	430	TTree t=(TTree)file->Get("TreeD_75x40_100x60");
	431
	432	AliSimDigits da, *digarr=&da;
	433	t->GetBranch("Segment")->SetAddress(&digarr);
	434	Stat_t num_of_entries=t->GetEntries();
	435
	436	Int_t digit_counter=0;
	437	Float_t xyz[3];
52a2a604	438	Double_t eta;
4de874d1	439
	440	Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
	441	printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
	442
	443	if(fRowContainer)
	444	delete [] fRowContainer;
	445	fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
	446	memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
	447
	448
52a2a604	449	//fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
52a2a604	450	//printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
4de874d1	451
52a2a604	452	/*
	453	if(fPhiRowContainer)
	454	delete [] fPhiRowContainer;
	455	fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
	456	memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
	457	*/
	458	fContainerBounds = (fNumEtaSegments+1)*(fNumOfPadRows+1);
	459	if(fVolume)
	460	delete [] fVolume;
	461	fVolume = new AliL3HoughContainer[fContainerBounds];
	462	memset(fVolume,0,fContainerBounds*sizeof(AliL3HoughContainer));
	463	/*
	464	Double_t phi_min = -10*ToRad;
	465	Double_t phi_max = 10*ToRad;
	466	Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
	467	*/
	468	Double_t eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
4de874d1	469	Int_t index;
	470	digit_counter=0;
	471
52a2a604	472	printf("\nLoading ALL pixels in slice\n\n");
52a2a604	473
4de874d1	474	for (Int_t i=0; i<num_of_entries; i++)
	475	{
	476	t->GetEvent(i);
	477	Int_t sector;
	478	Int_t row;
	479	param->AdjustSectorRow(digarr->GetID(),sector,row);
	480	Int_t slice,padrow;
	481	fTransform->Sector2Slice(slice,padrow,sector,row);
	482	if(slice != fSlice) continue;
	483	if(padrow < NRows[fPatch][0]) continue;
	484	if(padrow > NRows[fPatch][1]) break;
	485	digarr->First();
	486	do {
	487	Int_t time=digarr->CurrentRow();
	488	Int_t pad=digarr->CurrentColumn();
	489	Short_t signal=digarr->CurrentDigit();
	490	if(time < param->GetMaxTBin()-1 && time > 0)
	491	if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
	492	&& digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
	493	continue;
	494
	495
	496	fTransform->Raw2Global(xyz,sector,row,pad,time);
	497	eta = fTransform->GetEta(xyz);
52a2a604	498
4de874d1	499	if(eta < fEtaMin \|\| eta > fEtaMax) continue;
	500	fTransform->Global2Local(xyz,sector);
	501
	502
52a2a604	503	//phi = fTransform->GetPhi(xyz);
4de874d1	504	if(hist)
	505	hist->Fill(xyz[0],xyz[1],signal);
	506
	507	AliL3Digits *dig = new AliL3Digits;
	508	dig->fIndex = digit_counter;
	509	digit_counter++;
	510	dig->fCharge = signal;
	511	dig->fPad = pad;
	512	dig->fTime = time;
	513
	514	if(fRowContainer[padrow].first == NULL)
	515	fRowContainer[padrow].first = (void*)dig;
	516	else
	517	((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
	518	fRowContainer[padrow].last = (void*)dig;
	519
52a2a604	520	//thisHit->etaIndex=(Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1);
	521	Int_t eta_index = (Int_t)((eta-fEtaMin)/eta_slice);
	522	index = eta_index*fNumOfPadRows + padrow;
	523	if(index > fContainerBounds \|\| index < 0)
4de874d1	524	{
52a2a604	525
52a2a604	526	//printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
4de874d1	527	continue;
4de874d1	528	}
52a2a604	529
	530	if(fVolume[index].first == NULL)
	531	fVolume[index].first = (void*)dig;
4de874d1	532	else
52a2a604	533	((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
52a2a604	534	fVolume[index].last = (void*)dig;
4de874d1	535
52a2a604	536	/*
	537	Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
	538	index = phi_index*fNumOfPadRows + padrow;
	539	if(phi_index > fContainerBounds \|\| phi_index < 0)
	540	{
	541	printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
	542	continue;
	543	}
	544
	545	if(fPhiRowContainer[index].first == NULL)
	546	fPhiRowContainer[index].first = (void*)dig;
	547	else
	548	((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
	549	fPhiRowContainer[index].last=(void*)dig;
	550	*/
4de874d1	551	}while (digarr->Next());
	552
	553	}
	554
	555	fNDigits = digit_counter;
	556	printf("digitcounter %d\n",digit_counter);
	557	printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
	558	file->Close();
	559	delete file;
	560
	561	}