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

// @(#) $Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
//*-- Copyright &copy ALICE HLT Group

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3ClustFinderNew.h"
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
#include "AliL3SpacePointData.h"
#include "AliL3MemHandler.h"

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

/** \class AliL3ClustFinderNew
<pre>
//_____________________________________________________________
// AliL3ClustFinderNew
//
// The current cluster finder for HLT
// (Based on STAR L3)
// 
// The cluster finder is initialized with the Init function, 
// providing the slice and patch information to work on. 
// The input is a AliL3DigitRowData structure using the 
// Read function. The resulting space points will be in the
// array given by the SetOutputArray function.
// 
// There are several setters which control the behaviour:
//
// - SetXYError(Float_t):   set fixed error in XY direction
// - SetZError(Float_t):    set fixed error in Z  direction
//                            (used if errors are not calculated) 
// - SetDeconv(Bool_t):     switch on/off deconvolution
// - SetThreshold(UInt_t):  set charge threshold for cluster
// - SetMatchWidth(UInt_t): set the match distance in 
//                            time for sequences to be merged 
// - SetSTDOutput(Bool_t):  switch on/off output about found clusters   
// - SetCalcErr(Bool_t):    switch on/off calculation of 
//                          space point errors (or widths in raw system)
// - SetRawSP(Bool_t):      switch on/off convertion to raw system
//
//
// Example Usage:
//
// AliL3FileHandler *file = new AliL3FileHandler();
// file->SetAliInput(digitfile); //give some input file
// for(int slice=0; slice<=35; slice++){
//   for(int patch=0; pat<6; pat++){
//     file->Init(slice,patch);
//     UInt_t ndigits=0;
//     UInt_t maxclusters=100000;
//     UInt_t pointsize = maxclusters*sizeof(AliL3SpacePointData);
//     AliL3SpacePointData *points = (AliL3SpacePointData*)memory->Allocate(pointsize);
//     AliL3DigitRowData *digits = (AliL3DigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
//     AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
//     cf->SetMatchWidth(2);
//     cf->InitSlice(slice,patch,maxclusters);
//     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
//     cf->SetRawSP(kFALSE);       //Convert space points to local system
//     cf->SetThreshold(5);        //Threshold of cluster charge
//     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
//     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
//     cf->SetOutputArray(points); //Move the spacepoints to the array
//     cf->Read(ndigits,digits);   //give the data to the cf
//     cf->ProcessDigits();        //process the rows given by init
//     Int_t npoints = cf->GetNumberOfClusters();
//     AliL3MemHandler *out= new AliL3MemHandler();
//     out->SetBinaryOutput(fname);
//     out->Memory2Binary(npoints,points); //store the spacepoints
//     out->CloseBinaryOutput();
//     delete out;
//     file->free();
//     delete cf;
//   }
// }
</pre> 
*/

ClassImp(AliL3ClustFinderNew)

AliL3ClustFinderNew::AliL3ClustFinderNew()
{
  fMatch = 1;
  fThreshold = 10;
  fXYErr = 0.2;
  fZErr = 0.3;
  fDeconvPad = kTRUE;
  fDeconvTime = kTRUE;
  fStdout = kFALSE;
  fCalcerr = kTRUE;
  fRawSP = kFALSE;
  fFirstRow=0;
  fLastRow=0;
}

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

void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
{
  fNClusters = 0;
  fMaxNClusters = nmaxpoints;
  fCurrentSlice = slice;
  fCurrentPatch = patch;
  fFirstRow=AliL3Transform::GetFirstRow(patch);
  fLastRow=AliL3Transform::GetLastRow(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;
      //if(tempPt->fRow!=fCurrentRow) continue;
      ProcessRow(tempPt);
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
      tmp += size;
      tempPt = (AliL3DigitRowData*)tmp;
    }
  LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
    <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
    <<" patch "<<fCurrentPatch<<ENDLOG;
}

void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
{

  UInt_t last_pad = 123456789;

  ClusterData *pad1[5000]; //2 lists for internal memory=2pads
  ClusterData *pad2[5000]; //2 lists for internal memory=2pads
  ClusterData clusterlist[10000]; //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 the lists:
	  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,seq_error=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 >= AliL3Transform::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;
	  seq_error += data[bin].fTime*data[bin].fTime*charge;

	  //Check where to stop:
	  if(bin >= tempPt->fNDigit - 1) //out of range
	    break;
	  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;
      Int_t pad_error = data[bin].fPad*pad_mean;

      //Compare with results on previous pad:
      for(UInt_t p=0; p<n_previous; p++)
	{
	  //dont merge sequences on the same pad twice
	  if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;

	  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:
	      local->fTotalCharge += seq_charge;
	      local->fPad += pad_mean;
	      local->fPad2 += pad_error;
	      local->fTime += seq_average;
	      local->fTime2 += seq_error;
	      local->fMean = seq_mean;
	      local->fFlags++; //means we have more than one pad 
	      local->fLastMergedPad = data[bin].fPad;

	      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->fPad2 = pad_error;
	  tmp->fTime = seq_average;
	  tmp->fTime2 = seq_error;
	  tmp->fMean = seq_mean;
	  tmp->fFlags = 0;  //flags for single pad clusters
	  tmp->fLastMergedPad = data[bin].fPad;

	  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 single 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 fpad2=fXYErr*fXYErr; //fixed given error
      Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
      Float_t ftime2=fZErr*fZErr;  //fixed given error

      if(fCalcerr) { //calc the errors, otherwice take the fixed error 
	Int_t patch = AliL3Transform::GetPatch(fCurrentRow);

	Float_t sy2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
	if(sy2 < 0) {
	    LOG(AliL3Log::kError,"AliL3ClustFinderNew::WriteClusters","Cluster width")
	      <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	    continue;
	} else {
	  if(!fRawSP){
	    fpad2 = (sy2 + 1./12)*AliL3Transform::GetPadPitchWidth(patch)*AliL3Transform::GetPadPitchWidth(patch);
	    if(sy2 != 0){
	      fpad2*=0.108; //constants are from offline studies
	      if(patch<2)
		fpad2*=2.07;
	    }
	  } else fpad2=sy2; //take the width not the error
	}
	Float_t sz2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
	if(sz2 < 0){
	  LOG(AliL3Log::kError,"AliL3ClustFinderNew::WriteClusters","Cluster width")
	    <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	  continue;
	} else {
	  if(!fRawSP){
	    ftime2 = (sz2 + 1./12)*AliL3Transform::GetZWidth()*AliL3Transform::GetZWidth();
	    if(sz2 != 0) {
	      ftime2 *= 0.169; //constants are from offline studies
	      if(patch<2)
		ftime2 *= 1.77;
	    }
	  } else ftime2=sz2; //take the width, not the error
	}
      }
      if(fStdout==kTRUE)
	cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
      
      if(!fRawSP){
	AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
	AliL3Transform::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 "<<fNClusters<<ENDLOG;
	    return;
	  }  
	
	fSpacePointData[counter].fX = xyz[0];
	fSpacePointData[counter].fY = xyz[1];
	fSpacePointData[counter].fZ = xyz[2];
	
      } else {
	fSpacePointData[counter].fX = fCurrentRow;
	fSpacePointData[counter].fY = fpad;
	fSpacePointData[counter].fZ = ftime;
      }
      
      fSpacePointData[counter].fCharge = list[j].fTotalCharge;
      fSpacePointData[counter].fPadRow = fCurrentRow;
      fSpacePointData[counter].fSigmaY2 = fpad2;
      fSpacePointData[counter].fSigmaZ2  = ftime2;

      Int_t patch=fCurrentPatch;
      if(patch==-1) patch=0; //never store negative patch number
      fSpacePointData[counter].fID = counter
	+((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
#ifdef do_mc
      Int_t trackID[3];
      GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);

      fSpacePointData[counter].fTrackID[0] = trackID[0];
      fSpacePointData[counter].fTrackID[1] = trackID[1];
      fSpacePointData[counter].fTrackID[2] = trackID[2];

      //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
#endif
      
      fNClusters++;
      counter++;
    }
}

#ifdef do_mc
void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
{
  AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
  
  trackID[0]=trackID[1]=trackID[2]=-2;
  //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
  for(Int_t i=fFirstRow; i<=fLastRow; i++){
    if(rowPt->fRow < (UInt_t)fCurrentRow){
      AliL3MemHandler::UpdateRowPointer(rowPt);
      continue;
    }
    AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
    for(UInt_t j=0; j<rowPt->fNDigit; j++){
      Int_t cpad = digPt[j].fPad;
      Int_t ctime = digPt[j].fTime;
      if(cpad != pad) continue;
      if(ctime != time) continue;

      trackID[0] = digPt[j].fTrackID[0];
      trackID[1] = digPt[j].fTrackID[1];
      trackID[2] = digPt[j].fTrackID[2];
      
      //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
      break;
    }
    break;
  }
}
#endif
Commit	Line	Data
02f030e3	1	// @(#) $Id$
4ec96209	2
02f030e3	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
02f030e3	4	//*-- Copyright &copy ALICE HLT Group
46fbeb8a	5
118c26c3	6	#include "AliL3StandardIncludes.h"
118c26c3	7
46fbeb8a	8	#include "AliL3Logging.h"
	9	#include "AliL3ClustFinderNew.h"
	10	#include "AliL3DigitData.h"
	11	#include "AliL3Transform.h"
	12	#include "AliL3SpacePointData.h"
10f815d9	13	#include "AliL3MemHandler.h"
46fbeb8a	14
118c26c3	15	#if GCCVERSION == 3
	16	using namespace std;
	17	#endif
	18
b62fa469	19	/** \class AliL3ClustFinderNew
3e87ef69	20	<pre>
b661165c	21	//_____________________________________________________________
	22	// AliL3ClustFinderNew
	23	//
	24	// The current cluster finder for HLT
3e87ef69	25	// (Based on STAR L3)
	26	//
	27	// The cluster finder is initialized with the Init function,
	28	// providing the slice and patch information to work on.
	29	// The input is a AliL3DigitRowData structure using the
	30	// Read function. The resulting space points will be in the
	31	// array given by the SetOutputArray function.
	32	//
	33	// There are several setters which control the behaviour:
	34	//
	35	// - SetXYError(Float_t): set fixed error in XY direction
	36	// - SetZError(Float_t): set fixed error in Z direction
	37	// (used if errors are not calculated)
	38	// - SetDeconv(Bool_t): switch on/off deconvolution
	39	// - SetThreshold(UInt_t): set charge threshold for cluster
	40	// - SetMatchWidth(UInt_t): set the match distance in
	41	// time for sequences to be merged
	42	// - SetSTDOutput(Bool_t): switch on/off output about found clusters
	43	// - SetCalcErr(Bool_t): switch on/off calculation of
	44	// space point errors (or widths in raw system)
	45	// - SetRawSP(Bool_t): switch on/off convertion to raw system
	46	//
	47	//
	48	// Example Usage:
	49	//
	50	// AliL3FileHandler *file = new AliL3FileHandler();
	51	// file->SetAliInput(digitfile); //give some input file
	52	// for(int slice=0; slice<=35; slice++){
	53	// for(int patch=0; pat<6; pat++){
	54	// file->Init(slice,patch);
	55	// UInt_t ndigits=0;
	56	// UInt_t maxclusters=100000;
	57	// UInt_t pointsize = maxclusters*sizeof(AliL3SpacePointData);
	58	// AliL3SpacePointData points = (AliL3SpacePointData)memory->Allocate(pointsize);
	59	// AliL3DigitRowData digits = (AliL3DigitRowData)file->AliAltroDigits2Memory(ndigits,event);
	60	// AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
	61	// cf->SetMatchWidth(2);
	62	// cf->InitSlice(slice,patch,maxclusters);
	63	// cf->SetSTDOutput(kTRUE); //Some output to standard IO
	64	// cf->SetRawSP(kFALSE); //Convert space points to local system
	65	// cf->SetThreshold(5); //Threshold of cluster charge
	66	// cf->SetDeconv(kTRUE); //Deconv in pad and time direction
	67	// cf->SetCalcErr(kTRUE); //Calculate the errors of the spacepoints
	68	// cf->SetOutputArray(points); //Move the spacepoints to the array
	69	// cf->Read(ndigits,digits); //give the data to the cf
	70	// cf->ProcessDigits(); //process the rows given by init
	71	// Int_t npoints = cf->GetNumberOfClusters();
	72	// AliL3MemHandler *out= new AliL3MemHandler();
	73	// out->SetBinaryOutput(fname);
	74	// out->Memory2Binary(npoints,points); //store the spacepoints
	75	// out->CloseBinaryOutput();
	76	// delete out;
	77	// file->free();
	78	// delete cf;
	79	// }
	80	// }
	81	</pre>
	82	*/
b661165c	83
46fbeb8a	84	ClassImp(AliL3ClustFinderNew)
	85
	86	AliL3ClustFinderNew::AliL3ClustFinderNew()
	87	{
3e87ef69	88	fMatch = 1;
77f08cd3	89	fThreshold = 10;
	90	fXYErr = 0.2;
	91	fZErr = 0.3;
99a6d5c9	92	fDeconvPad = kTRUE;
2bbf57d8	93	fDeconvTime = kTRUE;
3e87ef69	94	fStdout = kFALSE;
	95	fCalcerr = kTRUE;
	96	fRawSP = kFALSE;
	97	fFirstRow=0;
	98	fLastRow=0;
46fbeb8a	99	}
46fbeb8a	100
46fbeb8a	101	AliL3ClustFinderNew::~AliL3ClustFinderNew()
46fbeb8a	102	{
46fbeb8a	103	}
	104
	105	void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
	106	{
	107	fNClusters = 0;
	108	fMaxNClusters = nmaxpoints;
	109	fCurrentSlice = slice;
	110	fCurrentPatch = patch;
	111	fFirstRow = firstrow;
	112	fLastRow = lastrow;
	113	}
	114
	115	void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
	116	{
	117	fNClusters = 0;
	118	fMaxNClusters = nmaxpoints;
	119	fCurrentSlice = slice;
	120	fCurrentPatch = patch;
3e87ef69	121	fFirstRow=AliL3Transform::GetFirstRow(patch);
3e87ef69	122	fLastRow=AliL3Transform::GetLastRow(patch);
46fbeb8a	123	}
	124
	125	void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
	126	{
46fbeb8a	127	fSpacePointData = pt;
	128	}
	129
46fbeb8a	130	void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
	131	{
	132	fNDigitRowData = ndigits;
	133	fDigitRowData = ptr;
46fbeb8a	134	}
	135
	136	void AliL3ClustFinderNew::ProcessDigits()
	137	{
	138	//Loop over rows, and call processrow
	139
	140
	141	AliL3DigitRowData tempPt = (AliL3DigitRowData)fDigitRowData;
	142
	143	for(Int_t i=fFirstRow; i<=fLastRow; i++)
	144	{
	145	fCurrentRow = i;
b2a02bce	146	//if(tempPt->fRow!=fCurrentRow) continue;
46fbeb8a	147	ProcessRow(tempPt);
	148	Byte_t tmp = (Byte_t)tempPt;
	149	Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
	150	tmp += size;
	151	tempPt = (AliL3DigitRowData*)tmp;
	152	}
355debe1	153	LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
	154	<<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
	155	<<" patch "<<fCurrentPatch<<ENDLOG;
46fbeb8a	156	}
46fbeb8a	157
46fbeb8a	158	void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
	159	{
	160
	161	UInt_t last_pad = 123456789;
	162
3e87ef69	163	ClusterData *pad1[5000]; //2 lists for internal memory=2pads
	164	ClusterData *pad2[5000]; //2 lists for internal memory=2pads
	165	ClusterData clusterlist[10000]; //Clusterlist
46fbeb8a	166
77f08cd3	167	ClusterData **currentPt; //List of pointers to the current pad
46fbeb8a	168	ClusterData **previousPt; //List of pointers to the previous pad
	169	currentPt = pad2;
	170	previousPt = pad1;
	171	UInt_t n_previous=0,n_current=0,n_total=0;
	172
	173	//Loop over sequences of this row:
	174	for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
	175	{
46fbeb8a	176	AliL3DigitData *data = tempPt->fDigitData;
	177	if(data[bin].fPad != last_pad)
	178	{
	179	//This is a new pad
46fbeb8a	180
3e87ef69	181	//Switch the lists:
46fbeb8a	182	if(currentPt == pad2)
	183	{
	184	currentPt = pad1;
	185	previousPt = pad2;
46fbeb8a	186	}
	187	else
	188	{
	189	currentPt = pad2;
	190	previousPt = pad1;
	191	}
	192	n_previous = n_current;
	193	n_current = 0;
99a6d5c9	194	if(bin[data].fPad != last_pad+1)
	195	{
	196	//this happens if there is a pad with no signal.
	197	n_previous = n_current = 0;
	198	}
	199	last_pad = data[bin].fPad;
46fbeb8a	200	}
	201
	202	Bool_t new_cluster = kTRUE;
c3dd27a3	203	UInt_t seq_charge=0,seq_average=0,seq_error=0;
99a6d5c9	204	UInt_t last_charge=0,last_was_falling=0;
99a6d5c9	205	Int_t new_bin=-1;
77f08cd3	206
99a6d5c9	207	if(fDeconvTime)
	208	{
	209	redo: //This is a goto.
	210	if(new_bin > -1)
	211	{
	212	bin = new_bin;
	213	new_bin = -1;
	214	}
	215
	216	last_charge=0;
	217	last_was_falling = 0;
	218	}
	219
46fbeb8a	220	while(1) //Loop over current sequence
46fbeb8a	221	{
494fad94	222	if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
adb0e845	223	{
	224	LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
	225	<<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
	226	break;
	227	}
	228
46fbeb8a	229	//Get the current ADC-value
46fbeb8a	230	UInt_t charge = data[bin].fCharge;
99a6d5c9	231
	232	if(fDeconvTime)
	233	{
	234	//Check if the last pixel in the sequence is smaller than this
	235	if(charge > last_charge)
	236	{
	237	if(last_was_falling)
	238	{
	239	new_bin = bin;
	240	break;
	241	}
	242	}
	243	else last_was_falling = 1; //last pixel was larger than this
	244	last_charge = charge;
	245	}
	246
46fbeb8a	247	//Sum the total charge of this sequence
46fbeb8a	248	seq_charge += charge;
46fbeb8a	249	seq_average += data[bin].fTime*charge;
c3dd27a3	250	seq_error += data[bin].fTimedata[bin].fTimecharge;
c3dd27a3	251
46fbeb8a	252	//Check where to stop:
3e87ef69	253	if(bin >= tempPt->fNDigit - 1) //out of range
3e87ef69	254	break;
46fbeb8a	255	if(data[bin+1].fPad != data[bin].fPad) //new pad
3e87ef69	256	break;
46fbeb8a	257	if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
	258	break;
	259
	260	bin++;
	261	}//end loop over sequence
	262
	263	//Calculate mean of sequence:
	264	Int_t seq_mean=0;
	265	if(seq_charge)
	266	seq_mean = seq_average/seq_charge;
	267	else
99a6d5c9	268	{
adb0e845	269	LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
adb0e845	270	<<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9	271	seq_mean = 1;
	272	seq_charge = 1;
	273	}
46fbeb8a	274
	275	//Calculate mean in pad direction:
	276	Int_t pad_mean = seq_charge*data[bin].fPad;
c3dd27a3	277	Int_t pad_error = data[bin].fPad*pad_mean;
3e87ef69	278
46fbeb8a	279	//Compare with results on previous pad:
	280	for(UInt_t p=0; p<n_previous; p++)
	281	{
b62fa469	282	//dont merge sequences on the same pad twice
	283	if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
	284
46fbeb8a	285	Int_t difference = seq_mean - previousPt[p]->fMean;
b62fa469	286	if(difference < -fMatch) break;
46fbeb8a	287
c3dd27a3	288	if(difference <= fMatch) //There is a match here!!
46fbeb8a	289	{
99a6d5c9	290	ClusterData *local = previousPt[p];
b62fa469	291
99a6d5c9	292	if(fDeconvPad)
	293	{
	294	if(seq_charge > local->fLastCharge)
	295	{
c3dd27a3	296	if(local->fChargeFalling) //The previous pad was falling
99a6d5c9	297	{
c3dd27a3	298	break; //create a new cluster
99a6d5c9	299	}
	300	}
	301	else
	302	local->fChargeFalling = 1;
	303	local->fLastCharge = seq_charge;
	304	}
	305
46fbeb8a	306	//Don't create a new cluster, because we found a match
	307	new_cluster = kFALSE;
	308
	309	//Update cluster on current pad with the matching one:
99a6d5c9	310	local->fTotalCharge += seq_charge;
99a6d5c9	311	local->fPad += pad_mean;
c3dd27a3	312	local->fPad2 += pad_error;
99a6d5c9	313	local->fTime += seq_average;
c3dd27a3	314	local->fTime2 += seq_error;
99a6d5c9	315	local->fMean = seq_mean;
99a6d5c9	316	local->fFlags++; //means we have more than one pad
b62fa469	317	local->fLastMergedPad = data[bin].fPad;
b62fa469	318
99a6d5c9	319	currentPt[n_current] = local;
46fbeb8a	320	n_current++;
46fbeb8a	321
46fbeb8a	322	break;
c3dd27a3	323	} //Checking for match at previous pad
c3dd27a3	324	} //Loop over results on previous pad.
46fbeb8a	325
	326	if(new_cluster)
	327	{
	328	//Start a new cluster. Add it to the clusterlist, and update
	329	//the list of pointers to clusters in current pad.
99a6d5c9	330	//current pad will be previous pad on next pad.
46fbeb8a	331
	332	//Add to the clusterlist:
	333	ClusterData *tmp = &clusterlist[n_total];
	334	tmp->fTotalCharge = seq_charge;
	335	tmp->fPad = pad_mean;
c3dd27a3	336	tmp->fPad2 = pad_error;
46fbeb8a	337	tmp->fTime = seq_average;
c3dd27a3	338	tmp->fTime2 = seq_error;
46fbeb8a	339	tmp->fMean = seq_mean;
3e87ef69	340	tmp->fFlags = 0; //flags for single pad clusters
b62fa469	341	tmp->fLastMergedPad = data[bin].fPad;
3e87ef69	342
99a6d5c9	343	if(fDeconvPad)
	344	{
	345	tmp->fChargeFalling = 0;
	346	tmp->fLastCharge = seq_charge;
	347	}
	348
46fbeb8a	349	//Update list of pointers to previous pad:
	350	currentPt[n_current] = &clusterlist[n_total];
	351	n_total++;
	352	n_current++;
	353	}
b62fa469	354
99a6d5c9	355	if(fDeconvTime)
99a6d5c9	356	if(new_bin >= 0) goto redo;
46fbeb8a	357	}//Loop over digits on this padrow
	358
	359	WriteClusters(n_total,clusterlist);
46fbeb8a	360	}
	361
	362	void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
	363	{
	364	Int_t thisrow,thissector;
	365	UInt_t counter = fNClusters;
	366
	367	for(int j=0; j<n_clusters; j++)
	368	{
3e87ef69	369	if(!list[j].fFlags) continue; //discard single pad clusters
46fbeb8a	370	if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
d3b44290	371
d3b44290	372	Float_t xyz[3];
c3dd27a3	373	Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
3e87ef69	374	Float_t fpad2=fXYErr*fXYErr; //fixed given error
c3dd27a3	375	Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
3e87ef69	376	Float_t ftime2=fZErr*fZErr; //fixed given error
	377
	378	if(fCalcerr) { //calc the errors, otherwice take the fixed error
	379	Int_t patch = AliL3Transform::GetPatch(fCurrentRow);
	380
	381	Float_t sy2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
	382	if(sy2 < 0) {
	383	LOG(AliL3Log::kError,"AliL3ClustFinderNew::WriteClusters","Cluster width")
	384	<<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	385	continue;
	386	} else {
	387	if(!fRawSP){
	388	fpad2 = (sy2 + 1./12)AliL3Transform::GetPadPitchWidth(patch)AliL3Transform::GetPadPitchWidth(patch);
	389	if(sy2 != 0){
	390	fpad2*=0.108; //constants are from offline studies
	391	if(patch<2)
	392	fpad2*=2.07;
	393	}
	394	} else fpad2=sy2; //take the width not the error
	395	}
	396	Float_t sz2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
	397	if(sz2 < 0){
	398	LOG(AliL3Log::kError,"AliL3ClustFinderNew::WriteClusters","Cluster width")
	399	<<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	400	continue;
	401	} else {
	402	if(!fRawSP){
	403	ftime2 = (sz2 + 1./12)AliL3Transform::GetZWidth()AliL3Transform::GetZWidth();
	404	if(sz2 != 0) {
	405	ftime2 *= 0.169; //constants are from offline studies
	406	if(patch<2)
	407	ftime2 *= 1.77;
	408	}
	409	} else ftime2=sz2; //take the width, not the error
	410	}
c3dd27a3	411	}
3e87ef69	412	if(fStdout==kTRUE)
b62fa469	413	cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
3e87ef69	414
	415	if(!fRawSP){
	416	AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
	417	AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
240d63be	418
3e87ef69	419	if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
	420	<<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
	421	if(fNClusters >= fMaxNClusters)
	422	{
	423	LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
	424	<<AliL3Log::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
	425	return;
	426	}
240d63be	427
3e87ef69	428	fSpacePointData[counter].fX = xyz[0];
	429	fSpacePointData[counter].fY = xyz[1];
	430	fSpacePointData[counter].fZ = xyz[2];
240d63be	431
3e87ef69	432	} else {
3e87ef69	433	fSpacePointData[counter].fX = fCurrentRow;
	434	fSpacePointData[counter].fY = fpad;
	435	fSpacePointData[counter].fZ = ftime;
	436	}
240d63be	437
240d63be	438	fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a	439	fSpacePointData[counter].fPadRow = fCurrentRow;
3e87ef69	440	fSpacePointData[counter].fSigmaY2 = fpad2;
3e87ef69	441	fSpacePointData[counter].fSigmaZ2 = ftime2;
02f030e3	442
	443	Int_t patch=fCurrentPatch;
	444	if(patch==-1) patch=0; //never store negative patch number
46fbeb8a	445	fSpacePointData[counter].fID = counter
02f030e3	446	+((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
10f815d9	447	#ifdef do_mc
	448	Int_t trackID[3];
	449	GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
3e87ef69	450
10f815d9	451	fSpacePointData[counter].fTrackID[0] = trackID[0];
	452	fSpacePointData[counter].fTrackID[1] = trackID[1];
	453	fSpacePointData[counter].fTrackID[2] = trackID[2];
9953e05f	454
3e87ef69	455	//cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
	456	#endif
	457
46fbeb8a	458	fNClusters++;
46fbeb8a	459	counter++;
46fbeb8a	460	}
46fbeb8a	461	}
10f815d9	462
	463	#ifdef do_mc
	464	void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
	465	{
	466	AliL3DigitRowData rowPt = (AliL3DigitRowData)fDigitRowData;
	467
355debe1	468	trackID[0]=trackID[1]=trackID[2]=-2;
10f815d9	469	//cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
3e87ef69	470	for(Int_t i=fFirstRow; i<=fLastRow; i++){
	471	if(rowPt->fRow < (UInt_t)fCurrentRow){
	472	AliL3MemHandler::UpdateRowPointer(rowPt);
	473	continue;
	474	}
	475	AliL3DigitData digPt = (AliL3DigitData)rowPt->fDigitData;
	476	for(UInt_t j=0; j<rowPt->fNDigit; j++){
	477	Int_t cpad = digPt[j].fPad;
	478	Int_t ctime = digPt[j].fTime;
	479	if(cpad != pad) continue;
	480	if(ctime != time) continue;
	481
	482	trackID[0] = digPt[j].fTrackID[0];
	483	trackID[1] = digPt[j].fTrackID[1];
	484	trackID[2] = digPt[j].fTrackID[2];
	485
	486	//cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
10f815d9	487	break;
10f815d9	488	}
3e87ef69	489	break;
3e87ef69	490	}
10f815d9	491	}
10f815d9	492	#endif