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

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

#include <math.h>
#include <TH1.h>
#include <TFile.h>

#include "AliL3MemHandler.h"
#include "GetGoodParticles.h"
#include "AliL3TrackArray.h"
#include "AliL3Logging.h"
#include "AliL3HoughEval.h"
#include "AliL3HoughTransformer.h"
#include "AliL3DigitData.h"
#include "AliL3HoughTrack.h"
#include "AliL3Transform.h"
#include "AliL3Histogram.h"
#include "AliL3Histogram1D.h"
#include "AliL3Defs.h"

ClassImp(AliL3HoughEval)

AliL3HoughEval::AliL3HoughEval()
{
  
  fTransform = new AliL3Transform();
  fRemoveFoundTracks = kFALSE;
  fNumOfPadsToLook = 1;
  fNumOfRowsToMiss = 1;
  fEtaHistos=0;
  
}


AliL3HoughEval::~AliL3HoughEval()
{
  fHoughTransformer = 0;
  if(fTransform)
    delete fTransform;
  if(fRowPointers)
    delete [] fRowPointers;
}

void AliL3HoughEval::InitTransformer(AliL3HoughTransformer *transformer)
{
  fHoughTransformer = transformer;
  fSlice = fHoughTransformer->GetSlice();
  fPatch = fHoughTransformer->GetPatch();
  fNrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
  fNEtaSegments = fHoughTransformer->GetNEtaSegments();
  fEtaMin = fHoughTransformer->GetEtaMin();
  fEtaMax = fHoughTransformer->GetEtaMax();
  GenerateLUT();
}

void AliL3HoughEval::GenerateLUT()
{
  //Generate a Look-up table, to limit the access to raw data
  
  fRowPointers = new AliL3DigitRowData*[fNrows];

  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
  if(!tempPt)
    printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
  
  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
    {
      Int_t prow = i - NRows[fPatch][0];
      fRowPointers[prow] = tempPt;
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }
  
}

Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
{
  //Look at rawdata along the road specified by the track candidates.
  //If track is good, return true, if not return false.
  
  Int_t sector,row;
  
  Int_t nrow=0,npixs=0;
  Float_t xyz[3];
  
  Int_t total_charge=0;//total charge along the road
  Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
    {
      Int_t prow = padrow - NRows[fPatch][0];
      
      if(!track->GetCrossingPoint(padrow,xyz))  
	{
	  printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	  continue;
	}
      
      fTransform->Slice2Sector(fSlice,padrow,sector,row);
      fTransform->Local2Raw(xyz,sector,row);
      npixs=0;
      
      //Get the timebins for this pad
      AliL3DigitRowData *tempPt = fRowPointers[prow];
      if(!tempPt) 
	{
	  printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
	  continue;
	}
      
      //Look at both sides of the pad:
      for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
	{
	  AliL3DigitData *digPt = tempPt->fDigitData;
	  for(UInt_t j=0; j<tempPt->fNDigit; j++)
	    {
	      if(digPt->fCharge <= fHoughTransformer->GetThreshold()) continue;
	      UChar_t pad = digPt[j].fPad;
	      if(pad < p) continue;
	      if(pad > p) break;
	      UShort_t time = digPt[j].fTime;
	      Double_t eta = fTransform->GetEta(padrow,pad,time);
	      Int_t pixel_index = (Int_t)(eta/etaslice);
	      if(pixel_index > eta_index) continue;
	      if(pixel_index != eta_index) break;
	      total_charge += digPt[j].fCharge;
	      if(remove)
		digPt[j].fCharge = 0; //Erease the track from image
	      npixs++;
	    }
	}
            
      if(npixs > 1)//At least 2 digits on this padrow
	{
	  nrow++;
	}	  
    }
  if(remove)
    return kTRUE;
  
  if(nrow >= fNrows - fNumOfRowsToMiss)//this was a good track
    {
      Double_t eta_track = (Double_t)eta_index*etaslice;
      track->SetEtaIndex(eta_index);
      track->SetWeight(total_charge,kTRUE);
      track->SetEta(eta_track);
      if(fRemoveFoundTracks)
	LookInsideRoad(track,eta_index,kTRUE);
      return kTRUE;
    }
  else
    return kFALSE;
}

void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
{
  
  Int_t sector,row;
  Float_t xyz[3];
  
  Int_t ntracks = tracks->GetNTracks();
  fEtaHistos = new AliL3Histogram1D*[ntracks];
  
  Char_t hname[100];
  for(Int_t i=0; i<ntracks; i++)
    {
      sprintf(hname,"etahist_%d",i);
      fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
    }
  Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
  
  for(Int_t ntr=0; ntr<ntracks; ntr++)
    {
      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
      if(!track) continue;
      for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	{
	  Int_t prow = padrow - NRows[fPatch][0];
	  
	  if(!track->GetCrossingPoint(padrow,xyz))  
	    {
	      printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	      continue;
	    }
	  
	  fTransform->Slice2Sector(fSlice,padrow,sector,row);
	  fTransform->Local2Raw(xyz,sector,row);
	  
	  //Get the timebins for this pad
	  AliL3DigitRowData *tempPt = fRowPointers[prow];
	  if(!tempPt) 
	    {
	      printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
	      continue;
	    }
	  
	  //Look at both sides of the pad:
	  for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
	    {
	      AliL3DigitData *digPt = tempPt->fDigitData;
	      for(UInt_t j=0; j<tempPt->fNDigit; j++)
		{
		  UChar_t pad = digPt[j].fPad;
		  
		  if(pad < p) continue;
		  if(pad > p) break;
		  UShort_t time = digPt[j].fTime;
		  Double_t eta = fTransform->GetEta(padrow,pad,time);
		  Int_t pixel_index = (Int_t)(eta/etaslice);
		  if(pixel_index > track->GetEtaIndex()+1) continue;
		  if(pixel_index < track->GetEtaIndex()-1) break;
		  fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
		}
	    }
	}
    }
  
  for(Int_t i=0; i<ntracks; i++)
    {
      AliL3Histogram1D *hist = fEtaHistos[i];
      Int_t max_bin = hist->GetMaximumBin();
      Double_t max_value = hist->GetBinContent(max_bin);
      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!track) continue;
      if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
	{
	  track->SetWeight((Int_t)max_value,kTRUE); 
	  track->SetEta(hist->GetBinCenter(max_bin));
	  track->SetNHits(track->GetWeight());
	}
      else
	{
	  track->SetWeight(0);
	  tracks->Remove(i); //remove this track, because it was not a peak
	}    
    }
  tracks->Compress();
  
  //for(Int_t i=0; i<ntracks; i++)
  //delete fEtaHistos[i];
Commit	Line	Data
	1	//Author: Anders Strand Vestbo
	2	//Last Modified: 28.6.01
	3
	4	#include <math.h>
	5	#include <TH1.h>
	6	#include <TFile.h>
	7
	8	#include "AliL3MemHandler.h"
	9	#include "GetGoodParticles.h"
	10	#include "AliL3TrackArray.h"
	11	#include "AliL3Logging.h"
	12	#include "AliL3HoughEval.h"
	13	#include "AliL3HoughTransformer.h"
	14	#include "AliL3DigitData.h"
	15	#include "AliL3HoughTrack.h"
	16	#include "AliL3Transform.h"
	17	#include "AliL3Histogram.h"
	18	#include "AliL3Histogram1D.h"
	19	#include "AliL3Defs.h"
	20
	21	ClassImp(AliL3HoughEval)
	22
	23	AliL3HoughEval::AliL3HoughEval()
	24	{
	25
	26	fTransform = new AliL3Transform();
	27	fRemoveFoundTracks = kFALSE;
	28	fNumOfPadsToLook = 1;
	29	fNumOfRowsToMiss = 1;
	30	fEtaHistos=0;
	31
	32	}
	33
	34
	35	AliL3HoughEval::~AliL3HoughEval()
	36	{
	37	fHoughTransformer = 0;
	38	if(fTransform)
	39	delete fTransform;
	40	if(fRowPointers)
	41	delete [] fRowPointers;
	42	}
	43
	44	void AliL3HoughEval::InitTransformer(AliL3HoughTransformer *transformer)
	45	{
	46	fHoughTransformer = transformer;
	47	fSlice = fHoughTransformer->GetSlice();
	48	fPatch = fHoughTransformer->GetPatch();
	49	fNrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
	50	fNEtaSegments = fHoughTransformer->GetNEtaSegments();
	51	fEtaMin = fHoughTransformer->GetEtaMin();
	52	fEtaMax = fHoughTransformer->GetEtaMax();
	53	GenerateLUT();
	54	}
	55
	56	void AliL3HoughEval::GenerateLUT()
	57	{
	58	//Generate a Look-up table, to limit the access to raw data
	59
	60	fRowPointers = new AliL3DigitRowData*[fNrows];
	61
	62	AliL3DigitRowData tempPt = (AliL3DigitRowData)fHoughTransformer->GetDataPointer();
	63	if(!tempPt)
	64	printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
	65
	66	for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
	67	{
	68	Int_t prow = i - NRows[fPatch][0];
	69	fRowPointers[prow] = tempPt;
	70	AliL3MemHandler::UpdateRowPointer(tempPt);
	71	}
	72
	73	}
	74
	75	Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
	76	{
	77	//Look at rawdata along the road specified by the track candidates.
	78	//If track is good, return true, if not return false.
	79
	80	Int_t sector,row;
	81
	82	Int_t nrow=0,npixs=0;
	83	Float_t xyz[3];
	84
	85	Int_t total_charge=0;//total charge along the road
	86	Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
	87	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	88	{
	89	Int_t prow = padrow - NRows[fPatch][0];
	90
	91	if(!track->GetCrossingPoint(padrow,xyz))
	92	{
	93	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	94	continue;
	95	}
	96
	97	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	98	fTransform->Local2Raw(xyz,sector,row);
	99	npixs=0;
	100
	101	//Get the timebins for this pad
	102	AliL3DigitRowData *tempPt = fRowPointers[prow];
	103	if(!tempPt)
	104	{
	105	printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
	106	continue;
	107	}
	108
	109	//Look at both sides of the pad:
	110	for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
	111	{
	112	AliL3DigitData *digPt = tempPt->fDigitData;
	113	for(UInt_t j=0; j<tempPt->fNDigit; j++)
	114	{
	115	if(digPt->fCharge <= fHoughTransformer->GetThreshold()) continue;
	116	UChar_t pad = digPt[j].fPad;
	117	if(pad < p) continue;
	118	if(pad > p) break;
	119	UShort_t time = digPt[j].fTime;
	120	Double_t eta = fTransform->GetEta(padrow,pad,time);
	121	Int_t pixel_index = (Int_t)(eta/etaslice);
	122	if(pixel_index > eta_index) continue;
	123	if(pixel_index != eta_index) break;
	124	total_charge += digPt[j].fCharge;
	125	if(remove)
	126	digPt[j].fCharge = 0; //Erease the track from image
	127	npixs++;
	128	}
	129	}
	130
	131	if(npixs > 1)//At least 2 digits on this padrow
	132	{
	133	nrow++;
	134	}
	135	}
	136	if(remove)
	137	return kTRUE;
	138
	139	if(nrow >= fNrows - fNumOfRowsToMiss)//this was a good track
	140	{
	141	Double_t eta_track = (Double_t)eta_index*etaslice;
	142	track->SetEtaIndex(eta_index);
	143	track->SetWeight(total_charge,kTRUE);
	144	track->SetEta(eta_track);
	145	if(fRemoveFoundTracks)
	146	LookInsideRoad(track,eta_index,kTRUE);
	147	return kTRUE;
	148	}
	149	else
	150	return kFALSE;
	151	}
	152
	153	void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
	154	{
	155
	156	Int_t sector,row;
	157	Float_t xyz[3];
	158
	159	Int_t ntracks = tracks->GetNTracks();
	160	fEtaHistos = new AliL3Histogram1D*[ntracks];
	161
	162	Char_t hname[100];
	163	for(Int_t i=0; i<ntracks; i++)
	164	{
	165	sprintf(hname,"etahist_%d",i);
	166	fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
	167	}
	168	Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
	169
	170	for(Int_t ntr=0; ntr<ntracks; ntr++)
	171	{
	172	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(ntr);
	173	if(!track) continue;
	174	for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
	175	{
	176	Int_t prow = padrow - NRows[fPatch][0];
	177
	178	if(!track->GetCrossingPoint(padrow,xyz))
	179	{
	180	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	181	continue;
	182	}
	183
	184	fTransform->Slice2Sector(fSlice,padrow,sector,row);
	185	fTransform->Local2Raw(xyz,sector,row);
	186
	187	//Get the timebins for this pad
	188	AliL3DigitRowData *tempPt = fRowPointers[prow];
	189	if(!tempPt)
	190	{
	191	printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
	192	continue;
	193	}
	194
	195	//Look at both sides of the pad:
	196	for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
	197	{
	198	AliL3DigitData *digPt = tempPt->fDigitData;
	199	for(UInt_t j=0; j<tempPt->fNDigit; j++)
	200	{
	201	UChar_t pad = digPt[j].fPad;
	202
	203	if(pad < p) continue;
	204	if(pad > p) break;
	205	UShort_t time = digPt[j].fTime;
	206	Double_t eta = fTransform->GetEta(padrow,pad,time);
	207	Int_t pixel_index = (Int_t)(eta/etaslice);
	208	if(pixel_index > track->GetEtaIndex()+1) continue;
	209	if(pixel_index < track->GetEtaIndex()-1) break;
	210	fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
	211	}
	212	}
	213	}
	214	}
	215
	216	for(Int_t i=0; i<ntracks; i++)
	217	{
	218	AliL3Histogram1D *hist = fEtaHistos[i];
	219	Int_t max_bin = hist->GetMaximumBin();
	220	Double_t max_value = hist->GetBinContent(max_bin);
	221	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	222	if(!track) continue;
	223	if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
	224	{
	225	track->SetWeight((Int_t)max_value,kTRUE);
	226	track->SetEta(hist->GetBinCenter(max_bin));
	227	track->SetNHits(track->GetWeight());
	228	}
	229	else
	230	{
	231	track->SetWeight(0);
	232	tracks->Remove(i); //remove this track, because it was not a peak
	233	}
	234	}
	235	tracks->Compress();
	236
	237	//for(Int_t i=0; i<ntracks; i++)
	238	//delete fEtaHistos[i];
	239