[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinderNew.cxx

//$Id$

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

#include <iostream.h>
#include "AliL3Logging.h"
#include "AliL3ClustFinderNew.h"
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
#include "AliL3SpacePointData.h"

//_____________________________________________________________
// AliL3ClustFinderNew
//
// The current cluster finder for HLT
// Based on STAR L3

ClassImp(AliL3ClustFinderNew)

AliL3ClustFinderNew::AliL3ClustFinderNew()
{
  fMatch = 4;
  fThreshold =10;
  fDeconvPad = kTRUE;
  fDeconvTime = kTRUE;
}

AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
{
  fTransform = transform;
  fMatch = 4;
  fThreshold =10;
  fDeconvTime = kTRUE;
  fDeconvPad = kTRUE;
}

AliL3ClustFinderNew::~AliL3ClustFinderNew()
{


}

void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
{
  fNClusters = 0;
  fMaxNClusters = nmaxpoints;
  fCurrentSlice = slice;
  fCurrentPatch = patch;
  fFirstRow = firstrow;
  fLastRow = lastrow;
  fDeconvTime = kTRUE;
  fDeconvPad = kTRUE;
}

void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
{
  fNClusters = 0;
  fMaxNClusters = nmaxpoints;
  fCurrentSlice = slice;
  fCurrentPatch = patch;
}

void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
{
  
  fSpacePointData = pt;
}


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

void AliL3ClustFinderNew::ProcessDigits()
{
  //Loop over rows, and call processrow
  
  
  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
  
  for(Int_t i=fFirstRow; i<=fLastRow; i++)
    {
      fCurrentRow = i;
      ProcessRow(tempPt);
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
    }
}


void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
{

  UInt_t last_pad = 123456789;

  ClusterData *pad1[2500]; //2 lists for internal memory=2pads
  ClusterData *pad2[2500]; //2 lists for internal memory=2pads
  ClusterData clusterlist[5000]; //Clusterlist

  ClusterData **currentPt; //List of pointers to the current pad
  ClusterData **previousPt; //List of pointers to the previous pad
  currentPt = pad2;
  previousPt = pad1;
  UInt_t n_previous=0,n_current=0,n_total=0;

  //Loop over sequences of this row:
  for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
    {
      
      AliL3DigitData *data = tempPt->fDigitData;
      if(data[bin].fPad != last_pad)
	{
	  //This is a new pad
	  
	  //Switch:
	  if(currentPt == pad2)
	    {
	      currentPt = pad1;
	      previousPt = pad2;
	      
	    }
	  else 
	    {
	      currentPt = pad2;
	      previousPt = pad1;
	    }
	  n_previous = n_current;
	  n_current = 0;
	  if(bin[data].fPad != last_pad+1)
	    {
	      //this happens if there is a pad with no signal.
	      n_previous = n_current = 0;
	    }
	  last_pad = data[bin].fPad;
	}

      Bool_t new_cluster = kTRUE;

      UInt_t seq_charge=0,seq_average=0;
      
      UInt_t last_charge=0,last_was_falling=0;
      Int_t new_bin=-1;
      if(fDeconvTime)
	{
	redo: //This is a goto.
	  if(new_bin > -1)
	    {
	      bin = new_bin;
	      new_bin = -1;
	    }
	  
	  last_charge=0;
	  last_was_falling = 0;
	}
      

      while(1) //Loop over current sequence
	{
	  if(data[bin].fTime >= fTransform->GetNTimeBins())
	    {
	      LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
		<<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
	      break;
	    }

	  //Get the current ADC-value
	  UInt_t charge = data[bin].fCharge;
	  
	  if(fDeconvTime)
	    {
	      //Check if the last pixel in the sequence is smaller than this
	      if(charge > last_charge)
		{
		  if(last_was_falling)
		    {
		      new_bin = bin;
		      break;
		    }
		}
	      else last_was_falling = 1; //last pixel was larger than this
	      last_charge = charge;
	    }
	  
	  //Sum the total charge of this sequence
	  seq_charge += charge;
	  seq_average += data[bin].fTime*charge;
	  
	  //Check where to stop:
	  if(data[bin+1].fPad != data[bin].fPad) //new pad
	    break; 
	  if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
	    break;
	  
	  bin++;
	}//end loop over sequence
      
      //Calculate mean of sequence:
      Int_t seq_mean=0;
      if(seq_charge)
	seq_mean = seq_average/seq_charge;
      else
	{
	  LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
	    <<"Error in data given to the cluster finder"<<ENDLOG;
	  seq_mean = 1;
	  seq_charge = 1;
	}
      
      //Calculate mean in pad direction:
      Int_t pad_mean = seq_charge*data[bin].fPad;

      //Compare with results on previous pad:
      for(UInt_t p=0; p<n_previous; p++)
	{
	  Int_t difference = seq_mean - previousPt[p]->fMean;
	  if(difference < -fMatch)
	    break;

	  if(difference <= fMatch)//There is a match here!!
	    {
	      
	      ClusterData *local = previousPt[p];
	      if(fDeconvPad)
		{
		  if(seq_charge > local->fLastCharge)
		    {
		      if(local->fChargeFalling)//The previous pad was falling
			{			
			  break;//create a new cluster
			}		    
		    }
		  else
		    local->fChargeFalling = 1;
		  local->fLastCharge = seq_charge;
		}
	      
	      //Don't create a new cluster, because we found a match
	      new_cluster = kFALSE;
	      
	      //Update cluster on current pad with the matching one:
	      //currentPt[n_current] = previousPt[p];
	      
	      local->fTotalCharge += seq_charge;
	      local->fPad += pad_mean;
	      local->fTime += seq_average;
	      local->fMean = seq_mean;
	      local->fFlags++; //means we have more than one pad 
	      
	      currentPt[n_current] = local;
	      n_current++;
	      
	      break;
	    }//Checking for match at previous pad
	}//Loop over results on previous pad.
      
      if(new_cluster)
	{
	  //Start a new cluster. Add it to the clusterlist, and update
	  //the list of pointers to clusters in current pad.
	  //current pad will be previous pad on next pad.

	  //Add to the clusterlist:
	  ClusterData *tmp = &clusterlist[n_total];
	  tmp->fTotalCharge = seq_charge;
	  tmp->fPad = pad_mean;
	  tmp->fTime = seq_average;
	  tmp->fMean = seq_mean;
	  tmp->fFlags = 0;  //flags for 1 pad clusters
	  if(fDeconvPad)
	    {
	      tmp->fChargeFalling = 0;
	      tmp->fLastCharge = seq_charge;
	    }

	  //Update list of pointers to previous pad:
	  currentPt[n_current] = &clusterlist[n_total];
	  n_total++;
	  n_current++;
	}
      if(fDeconvTime)
	if(new_bin >= 0) goto redo;
      
    }//Loop over digits on this padrow
  
  WriteClusters(n_total,clusterlist);


}

void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
{
  Int_t thisrow,thissector;
  UInt_t counter = fNClusters;
  
  for(int j=0; j<n_clusters; j++)
    {
      if(!list[j].fFlags) continue; //discard 1 pad clusters
      if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
      Float_t xyz[3];
      
      Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
      Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
      //printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
      fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
      fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
      if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
	<<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
      if(fNClusters >= fMaxNClusters)
	{
	  LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
	    <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
	  return;
	  }  
      fSpacePointData[counter].fCharge = list[j].fTotalCharge;
      fSpacePointData[counter].fX = xyz[0];
      fSpacePointData[counter].fY = xyz[1];
      fSpacePointData[counter].fZ = xyz[2];
      fSpacePointData[counter].fPadRow = fCurrentRow;
      fSpacePointData[counter].fXYErr = fXYErr;
      fSpacePointData[counter].fZErr = fZErr;
      fSpacePointData[counter].fID = counter
	+((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
      
      fNClusters++;
      counter++;

    }

}
Commit	Line	Data
4ec96209	1	//$Id$
4ec96209	2
b661165c	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
b661165c	4	//*-- Copyright &copy ASV
46fbeb8a	5
	6	#include <iostream.h>
	7	#include "AliL3Logging.h"
	8	#include "AliL3ClustFinderNew.h"
	9	#include "AliL3DigitData.h"
	10	#include "AliL3Transform.h"
	11	#include "AliL3SpacePointData.h"
	12
b661165c	13	//_____________________________________________________________
	14	// AliL3ClustFinderNew
	15	//
	16	// The current cluster finder for HLT
	17	// Based on STAR L3
	18
46fbeb8a	19	ClassImp(AliL3ClustFinderNew)
	20
	21	AliL3ClustFinderNew::AliL3ClustFinderNew()
	22	{
	23	fMatch = 4;
	24	fThreshold =10;
99a6d5c9	25	fDeconvPad = kTRUE;
2bbf57d8	26	fDeconvTime = kTRUE;
46fbeb8a	27	}
	28
	29	AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
	30	{
	31	fTransform = transform;
	32	fMatch = 4;
	33	fThreshold =10;
99a6d5c9	34	fDeconvTime = kTRUE;
99a6d5c9	35	fDeconvPad = kTRUE;
46fbeb8a	36	}
	37
	38	AliL3ClustFinderNew::~AliL3ClustFinderNew()
	39	{
	40
	41
	42	}
	43
	44	void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
	45	{
	46	fNClusters = 0;
	47	fMaxNClusters = nmaxpoints;
	48	fCurrentSlice = slice;
	49	fCurrentPatch = patch;
	50	fFirstRow = firstrow;
	51	fLastRow = lastrow;
99a6d5c9	52	fDeconvTime = kTRUE;
99a6d5c9	53	fDeconvPad = kTRUE;
46fbeb8a	54	}
	55
	56	void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
	57	{
	58	fNClusters = 0;
	59	fMaxNClusters = nmaxpoints;
	60	fCurrentSlice = slice;
	61	fCurrentPatch = patch;
	62	}
	63
	64	void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
	65	{
	66
	67	fSpacePointData = pt;
	68	}
	69
	70
	71	void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
	72	{
	73	fNDigitRowData = ndigits;
	74	fDigitRowData = ptr;
	75
	76	}
	77
	78	void AliL3ClustFinderNew::ProcessDigits()
	79	{
	80	//Loop over rows, and call processrow
	81
	82
	83	AliL3DigitRowData tempPt = (AliL3DigitRowData)fDigitRowData;
	84
	85	for(Int_t i=fFirstRow; i<=fLastRow; i++)
	86	{
	87	fCurrentRow = i;
	88	ProcessRow(tempPt);
	89	Byte_t tmp = (Byte_t)tempPt;
	90	Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
	91	tmp += size;
	92	tempPt = (AliL3DigitRowData*)tmp;
	93	}
	94	}
	95
	96
	97
	98	void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
	99	{
	100
	101	UInt_t last_pad = 123456789;
	102
	103	ClusterData *pad1[2500]; //2 lists for internal memory=2pads
	104	ClusterData *pad2[2500]; //2 lists for internal memory=2pads
	105	ClusterData clusterlist[5000]; //Clusterlist
	106
	107	ClusterData **currentPt; //List of pointers to the current pad
	108	ClusterData **previousPt; //List of pointers to the previous pad
	109	currentPt = pad2;
	110	previousPt = pad1;
	111	UInt_t n_previous=0,n_current=0,n_total=0;
	112
	113	//Loop over sequences of this row:
	114	for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
	115	{
	116
	117	AliL3DigitData *data = tempPt->fDigitData;
118	if(data[bin].fPad != last_pad)
119	{
120	//This is a new pad
46fbeb8a	121
	122	//Switch:
	123	if(currentPt == pad2)
	124	{
	125	currentPt = pad1;
	126	previousPt = pad2;
	127
	128	}
	129	else
	130	{
	131	currentPt = pad2;
	132	previousPt = pad1;
	133	}
	134	n_previous = n_current;
	135	n_current = 0;
99a6d5c9	136	if(bin[data].fPad != last_pad+1)
	137	{
	138	//this happens if there is a pad with no signal.
	139	n_previous = n_current = 0;
	140	}
	141	last_pad = data[bin].fPad;
46fbeb8a	142	}
	143
	144	Bool_t new_cluster = kTRUE;
	145
	146	UInt_t seq_charge=0,seq_average=0;
	147
99a6d5c9	148	UInt_t last_charge=0,last_was_falling=0;
	149	Int_t new_bin=-1;
	150	if(fDeconvTime)
	151	{
	152	redo: //This is a goto.
	153	if(new_bin > -1)
	154	{
	155	bin = new_bin;
	156	new_bin = -1;
	157	}
	158
	159	last_charge=0;
	160	last_was_falling = 0;
	161	}
	162
	163
46fbeb8a	164	while(1) //Loop over current sequence
46fbeb8a	165	{
adb0e845	166	if(data[bin].fTime >= fTransform->GetNTimeBins())
	167	{
	168	LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
	169	<<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
	170	break;
	171	}
	172
46fbeb8a	173	//Get the current ADC-value
46fbeb8a	174	UInt_t charge = data[bin].fCharge;
99a6d5c9	175
	176	if(fDeconvTime)
	177	{
	178	//Check if the last pixel in the sequence is smaller than this
	179	if(charge > last_charge)
	180	{
	181	if(last_was_falling)
	182	{
	183	new_bin = bin;
	184	break;
	185	}
	186	}
	187	else last_was_falling = 1; //last pixel was larger than this
	188	last_charge = charge;
	189	}
	190
46fbeb8a	191	//Sum the total charge of this sequence
46fbeb8a	192	seq_charge += charge;
46fbeb8a	193	seq_average += data[bin].fTime*charge;
	194
	195	//Check where to stop:
	196	if(data[bin+1].fPad != data[bin].fPad) //new pad
	197	break;
	198	if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
	199	break;
	200
	201	bin++;
	202	}//end loop over sequence
	203
	204	//Calculate mean of sequence:
	205	Int_t seq_mean=0;
	206	if(seq_charge)
	207	seq_mean = seq_average/seq_charge;
	208	else
99a6d5c9	209	{
adb0e845	210	LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
adb0e845	211	<<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9	212	seq_mean = 1;
	213	seq_charge = 1;
	214	}
46fbeb8a	215
	216	//Calculate mean in pad direction:
	217	Int_t pad_mean = seq_charge*data[bin].fPad;
	218
	219	//Compare with results on previous pad:
	220	for(UInt_t p=0; p<n_previous; p++)
	221	{
	222	Int_t difference = seq_mean - previousPt[p]->fMean;
	223	if(difference < -fMatch)
	224	break;
	225
	226	if(difference <= fMatch)//There is a match here!!
	227	{
99a6d5c9	228
	229	ClusterData *local = previousPt[p];
	230	if(fDeconvPad)
	231	{
	232	if(seq_charge > local->fLastCharge)
	233	{
	234	if(local->fChargeFalling)//The previous pad was falling
	235	{
	236	break;//create a new cluster
	237	}
	238	}
	239	else
	240	local->fChargeFalling = 1;
	241	local->fLastCharge = seq_charge;
	242	}
	243
46fbeb8a	244	//Don't create a new cluster, because we found a match
	245	new_cluster = kFALSE;
	246
	247	//Update cluster on current pad with the matching one:
99a6d5c9	248	//currentPt[n_current] = previousPt[p];
	249
	250	local->fTotalCharge += seq_charge;
	251	local->fPad += pad_mean;
	252	local->fTime += seq_average;
	253	local->fMean = seq_mean;
	254	local->fFlags++; //means we have more than one pad
	255
	256	currentPt[n_current] = local;
46fbeb8a	257	n_current++;
46fbeb8a	258
46fbeb8a	259	break;
	260	}//Checking for match at previous pad
	261	}//Loop over results on previous pad.
	262
	263	if(new_cluster)
	264	{
	265	//Start a new cluster. Add it to the clusterlist, and update
	266	//the list of pointers to clusters in current pad.
99a6d5c9	267	//current pad will be previous pad on next pad.
46fbeb8a	268
	269	//Add to the clusterlist:
	270	ClusterData *tmp = &clusterlist[n_total];
	271	tmp->fTotalCharge = seq_charge;
	272	tmp->fPad = pad_mean;
	273	tmp->fTime = seq_average;
	274	tmp->fMean = seq_mean;
	275	tmp->fFlags = 0; //flags for 1 pad clusters
99a6d5c9	276	if(fDeconvPad)
	277	{
	278	tmp->fChargeFalling = 0;
	279	tmp->fLastCharge = seq_charge;
	280	}
	281
46fbeb8a	282	//Update list of pointers to previous pad:
	283	currentPt[n_current] = &clusterlist[n_total];
	284	n_total++;
	285	n_current++;
	286	}
99a6d5c9	287	if(fDeconvTime)
	288	if(new_bin >= 0) goto redo;
	289
46fbeb8a	290	}//Loop over digits on this padrow
	291
	292	WriteClusters(n_total,clusterlist);
	293
	294
	295	}
	296
	297	void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
	298	{
	299	Int_t thisrow,thissector;
	300	UInt_t counter = fNClusters;
	301
	302	for(int j=0; j<n_clusters; j++)
	303	{
	304	if(!list[j].fFlags) continue; //discard 1 pad clusters
	305	if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
	306	Float_t xyz[3];
	307
	308	Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
	309	Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
	310	//printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
	311	fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
	312	fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
	313	if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
	314	<<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
	315	if(fNClusters >= fMaxNClusters)
	316	{
	317	LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
	318	<<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
	319	return;
	320	}
95a00d93	321	fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a	322	fSpacePointData[counter].fX = xyz[0];
	323	fSpacePointData[counter].fY = xyz[1];
	324	fSpacePointData[counter].fZ = xyz[2];
	325	fSpacePointData[counter].fPadRow = fCurrentRow;
	326	fSpacePointData[counter].fXYErr = fXYErr;
	327	fSpacePointData[counter].fZErr = fZErr;
	328	fSpacePointData[counter].fID = counter
	329	+((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
	330
46fbeb8a	331	fNClusters++;
	332	counter++;
	333
	334	}
	335
	336	}