[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClustFinderNew.cxx

// @(#) $Id$

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

#include "AliHLTTPCStandardIncludes.h"
#include "AliHLTTPCRootTypes.h"
#include "AliHLTTPCLogging.h"
#include "AliHLTTPCClustFinderNew.h"
#include "AliHLTTPCDigitData.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCMemHandler.h"

#if __GNUC__ >= 3
using namespace std;
#endif

/** \class AliHLTTPCClustFinderNew
<pre>
//_____________________________________________________________
// AliHLTTPCClustFinderNew
//
// 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 AliHLTTPCDigitRowData 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:
//
// AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
// 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(AliHLTTPCSpacePointData);
//     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
//     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
//     AliHLTTPCClustFinderNew *cf = new AliHLTTPCClustFinderNew();
//     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();
//     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
//     out->SetBinaryOutput(fname);
//     out->Memory2Binary(npoints,points); //store the spacepoints
//     out->CloseBinaryOutput();
//     delete out;
//     file->free();
//     delete cf;
//   }
// }
</pre> 
*/

ClassImp(AliHLTTPCClustFinderNew)

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

AliHLTTPCClustFinderNew::~AliHLTTPCClustFinderNew()
{
  //destructor
  ;
}

void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
{
  //init slice
  fNClusters = 0;
  fMaxNClusters = nmaxpoints;
  fCurrentSlice = slice;
  fCurrentPatch = patch;
  fFirstRow = firstrow;
  fLastRow = lastrow;
}

void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
{
  //init slice
  fNClusters = 0;
  fMaxNClusters = nmaxpoints;
  fCurrentSlice = slice;
  fCurrentPatch = patch;
  fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
  fLastRow=AliHLTTPCTransform::GetLastRow(patch);
}

void AliHLTTPCClustFinderNew::SetOutputArray(AliHLTTPCSpacePointData *pt)
{
  //set pointer to output
  fSpacePointData = pt;
}

void AliHLTTPCClustFinderNew::Read(UInt_t ndigits,AliHLTTPCDigitRowData *ptr)
{
  //set input pointer
  fNDigitRowData = ndigits;
  fDigitRowData = ptr;
}

void AliHLTTPCClustFinderNew::ProcessDigits()
{
  //Loop over rows, and call processrow
  AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fDigitRowData;
  fNClusters = 0; 
  for(Int_t i=fFirstRow; i<=fLastRow; i++)
    {
      fCurrentRow = i;
      if((Int_t)tempPt->fRow!=fCurrentRow){
	LOG(AliHLTTPCLog::kWarning,"AliHLTTPCClustFinderNew::ProcessDigits","Digits")
	  <<"Row number should match! "<<tempPt->fRow<<" "<<fCurrentRow<<ENDLOG;
	continue;
      }
      ProcessRow(tempPt);
      Byte_t *tmp = (Byte_t*)tempPt;
      Int_t size = sizeof(AliHLTTPCDigitRowData) + tempPt->fNDigit*sizeof(AliHLTTPCDigitData);
      tmp += size;
      tempPt = (AliHLTTPCDigitRowData*)tmp;
    }
  LOG(AliHLTTPCLog::kInformational,"AliHLTTPCClustFinderNew::ProcessDigits","Space points")
    <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
    <<" patch "<<fCurrentPatch<<ENDLOG;
}

void AliHLTTPCClustFinderNew::ProcessRow(AliHLTTPCDigitRowData *tempPt)
{
  //process row
  UInt_t lastpad = 123456789;

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

  AliClusterData **currentPt;  //List of pointers to the current pad
  AliClusterData **previousPt; //List of pointers to the previous pad
  currentPt = pad2;
  previousPt = pad1;
  UInt_t nprevious=0,ncurrent=0,ntotal=0;

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

      Bool_t newcluster = kTRUE;
      UInt_t seqcharge=0,seqaverage=0,seqerror=0;
      UInt_t lastcharge=0,lastwas_falling=0;
      Int_t newbin=-1;

      if(fDeconvTime)
	{
	redo: //This is a goto.
	  if(newbin > -1)
	    {
	      bin = newbin;
	      newbin = -1;
	    }
	  
	  lastcharge=0;
	  lastwas_falling = 0;
	}
      
      while(1) //Loop over current sequence
	{
	  if(data[bin].fTime >= AliHLTTPCTransform::GetNTimeBins())
	    {
	      LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::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 > lastcharge)
		{
		  if(lastwas_falling)
		    {
		      newbin = bin;
		      break;
		    }
		}
	      else lastwas_falling = 1; //last pixel was larger than this
	      lastcharge = charge;
	    }
	  
	  //Sum the total charge of this sequence
	  seqcharge += charge;
	  seqaverage += data[bin].fTime*charge;
	  seqerror += 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 seqmean=0;
      if(seqcharge)
	seqmean = seqaverage/seqcharge;
      else
	{
	  LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Data")
	    <<"Error in data given to the cluster finder"<<ENDLOG;
	  seqmean = 1;
	  seqcharge = 1;
	}
      
      //Calculate mean in pad direction:
      Int_t padmean = seqcharge*data[bin].fPad;
      Int_t paderror = data[bin].fPad*padmean;

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

	  Int_t difference = seqmean - previousPt[p]->fMean;
	  if(difference < -fMatch) break;

	  if(difference <= fMatch) //There is a match here!!
	    {
	      AliClusterData *local = previousPt[p];

	      if(fDeconvPad)
		{
		  if(seqcharge > local->fLastCharge)
		    {
		      if(local->fChargeFalling) //The previous pad was falling
			{			
			  break; //create a new cluster
			}		    
		    }
		  else
		    local->fChargeFalling = 1;
		  local->fLastCharge = seqcharge;
		}
	      
	      //Don't create a new cluster, because we found a match
	      newcluster = kFALSE;
	      
	      //Update cluster on current pad with the matching one:
	      local->fTotalCharge += seqcharge;
	      local->fPad += padmean;
	      local->fPad2 += paderror;
	      local->fTime += seqaverage;
	      local->fTime2 += seqerror;
	      local->fMean = seqmean;
	      local->fFlags++; //means we have more than one pad 
	      local->fLastMergedPad = data[bin].fPad;

	      currentPt[ncurrent] = local;
	      ncurrent++;
	      
	      break;
	    } //Checking for match at previous pad
	} //Loop over results on previous pad.
      
      if(newcluster)
	{
	  //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:
	  AliClusterData *tmp = &clusterlist[ntotal];
	  tmp->fTotalCharge = seqcharge;
	  tmp->fPad = padmean;
	  tmp->fPad2 = paderror;
	  tmp->fTime = seqaverage;
	  tmp->fTime2 = seqerror;
	  tmp->fMean = seqmean;
	  tmp->fFlags = 0;  //flags for single pad clusters
	  tmp->fLastMergedPad = data[bin].fPad;

	  if(fDeconvPad)
	    {
	      tmp->fChargeFalling = 0;
	      tmp->fLastCharge = seqcharge;
	    }

	  //Update list of pointers to previous pad:
	  currentPt[ncurrent] = &clusterlist[ntotal];
	  ntotal++;
	  ncurrent++;
	}

      if(fDeconvTime)
	if(newbin >= 0) goto redo;
    }//Loop over digits on this padrow
  
  WriteClusters(ntotal,clusterlist);
}

void AliHLTTPCClustFinderNew::WriteClusters(Int_t nclusters,AliClusterData *list)
{
  //write cluster to output pointer
  Int_t thisrow,thissector;
  UInt_t counter = fNClusters;
  
  for(int j=0; j<nclusters; 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 / list[j].fTotalCharge;
      Float_t fpad2=fXYErr*fXYErr; //fixed given error
      Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
      Float_t ftime2=fZErr*fZErr;  //fixed given error

      if(fCalcerr) { //calc the errors, otherwice take the fixed error 
	Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
	UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
	Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
	sy2/=q2;
	if(sy2 < 0) {
	    LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
	      <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	    continue;
	} else {
	  if(!fRawSP){
	    fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::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=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
	sz2/=q2;
	if(sz2 < 0){
	  LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
	    <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	  continue;
	} else {
	  if(!fRawSP){
	    ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::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){
	AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
	AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
	
	if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
	  <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
	if(fNClusters >= fMaxNClusters)
	  {
	    LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
	      <<AliHLTTPCLog::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 AliHLTTPCClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
{
  //get mc id
  AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)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){
      AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
      continue;
    }
    AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)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
a6c02c85	1	// @(#) $Id$
	2
	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
	4	//*-- Copyright &copy ALICE HLT Group
	5
	6	#include "AliHLTTPCStandardIncludes.h"
	7	#include "AliHLTTPCRootTypes.h"
	8	#include "AliHLTTPCLogging.h"
	9	#include "AliHLTTPCClustFinderNew.h"
	10	#include "AliHLTTPCDigitData.h"
	11	#include "AliHLTTPCTransform.h"
	12	#include "AliHLTTPCSpacePointData.h"
	13	#include "AliHLTTPCMemHandler.h"
	14
	15	#if __GNUC__ >= 3
	16	using namespace std;
	17	#endif
	18
	19	/** \class AliHLTTPCClustFinderNew
	20	<pre>
	21	//_____________________________________________________________
	22	// AliHLTTPCClustFinderNew
	23	//
	24	// The current cluster finder for HLT
	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 AliHLTTPCDigitRowData 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	// AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
	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(AliHLTTPCSpacePointData);
	58	// AliHLTTPCSpacePointData points = (AliHLTTPCSpacePointData)memory->Allocate(pointsize);
	59	// AliHLTTPCDigitRowData digits = (AliHLTTPCDigitRowData)file->AliAltroDigits2Memory(ndigits,event);
	60	// AliHLTTPCClustFinderNew *cf = new AliHLTTPCClustFinderNew();
	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	// AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
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	*/
83
84	ClassImp(AliHLTTPCClustFinderNew)
85
86	AliHLTTPCClustFinderNew::AliHLTTPCClustFinderNew()
87	{
88	//constructor
89	fMatch = 1;
90	fThreshold = 10;
91	fXYErr = 0.2;
92	fZErr = 0.3;
93	fDeconvPad = kTRUE;
94	fDeconvTime = kTRUE;
95	fStdout = kFALSE;
96	fCalcerr = kTRUE;
97	fRawSP = kFALSE;
98	fFirstRow=0;
99	fLastRow=0;
100	}
101
102	AliHLTTPCClustFinderNew::~AliHLTTPCClustFinderNew()
103	{
104	//destructor
105	;
106	}
107
108	void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
109	{
110	//init slice
111	fNClusters = 0;
112	fMaxNClusters = nmaxpoints;
113	fCurrentSlice = slice;
114	fCurrentPatch = patch;
115	fFirstRow = firstrow;
116	fLastRow = lastrow;
117	}
118
119	void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
120	{
121	//init slice
122	fNClusters = 0;
123	fMaxNClusters = nmaxpoints;
124	fCurrentSlice = slice;
125	fCurrentPatch = patch;
126	fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
127	fLastRow=AliHLTTPCTransform::GetLastRow(patch);
128	}
129
130	void AliHLTTPCClustFinderNew::SetOutputArray(AliHLTTPCSpacePointData *pt)
131	{
132	//set pointer to output
133	fSpacePointData = pt;
134	}
135
136	void AliHLTTPCClustFinderNew::Read(UInt_t ndigits,AliHLTTPCDigitRowData *ptr)
137	{
138	//set input pointer
139	fNDigitRowData = ndigits;
140	fDigitRowData = ptr;
141	}
142
143	void AliHLTTPCClustFinderNew::ProcessDigits()
144	{
145	//Loop over rows, and call processrow
146	AliHLTTPCDigitRowData tempPt = (AliHLTTPCDigitRowData)fDigitRowData;
147	fNClusters = 0;
148	for(Int_t i=fFirstRow; i<=fLastRow; i++)
149	{
150	fCurrentRow = i;
151	if((Int_t)tempPt->fRow!=fCurrentRow){
152	LOG(AliHLTTPCLog::kWarning,"AliHLTTPCClustFinderNew::ProcessDigits","Digits")
153	<<"Row number should match! "<<tempPt->fRow<<" "<<fCurrentRow<<ENDLOG;
154	continue;
155	}
156	ProcessRow(tempPt);
157	Byte_t tmp = (Byte_t)tempPt;
158	Int_t size = sizeof(AliHLTTPCDigitRowData) + tempPt->fNDigit*sizeof(AliHLTTPCDigitData);
159	tmp += size;
160	tempPt = (AliHLTTPCDigitRowData*)tmp;
161	}
162	LOG(AliHLTTPCLog::kInformational,"AliHLTTPCClustFinderNew::ProcessDigits","Space points")
163	<<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
164	<<" patch "<<fCurrentPatch<<ENDLOG;
165	}
166
167	void AliHLTTPCClustFinderNew::ProcessRow(AliHLTTPCDigitRowData *tempPt)
168	{
169	//process row
170	UInt_t lastpad = 123456789;
171
172	AliClusterData *pad1[5000]; //2 lists for internal memory=2pads
173	AliClusterData *pad2[5000]; //2 lists for internal memory=2pads
174	AliClusterData clusterlist[10000]; //Clusterlist
175
176	AliClusterData **currentPt; //List of pointers to the current pad
177	AliClusterData **previousPt; //List of pointers to the previous pad
178	currentPt = pad2;
179	previousPt = pad1;
180	UInt_t nprevious=0,ncurrent=0,ntotal=0;
181
182	//Loop over sequences of this row:
183	for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
184	{
185	AliHLTTPCDigitData *data = tempPt->fDigitData;
186	if(data[bin].fPad != lastpad)
187	{
188	//This is a new pad
189
190	//Switch the lists:
191	if(currentPt == pad2)
192	{
193	currentPt = pad1;
194	previousPt = pad2;
195	}
196	else
197	{
198	currentPt = pad2;
199	previousPt = pad1;
200	}
201	nprevious = ncurrent;
202	ncurrent = 0;
203	if(bin[data].fPad != lastpad+1)
204	{
205	//this happens if there is a pad with no signal.
206	nprevious = ncurrent = 0;
207	}
208	lastpad = data[bin].fPad;
209	}
210
211	Bool_t newcluster = kTRUE;
212	UInt_t seqcharge=0,seqaverage=0,seqerror=0;
213	UInt_t lastcharge=0,lastwas_falling=0;
214	Int_t newbin=-1;
215
216	if(fDeconvTime)
217	{
218	redo: //This is a goto.
219	if(newbin > -1)
220	{
221	bin = newbin;
222	newbin = -1;
223	}
224
225	lastcharge=0;
226	lastwas_falling = 0;
227	}
228
229	while(1) //Loop over current sequence
230	{
231	if(data[bin].fTime >= AliHLTTPCTransform::GetNTimeBins())
232	{
233	LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Digits")
234	<<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
235	break;
236	}
237
238	//Get the current ADC-value
239	UInt_t charge = data[bin].fCharge;
240
241	if(fDeconvTime)
242	{
243	//Check if the last pixel in the sequence is smaller than this
244	if(charge > lastcharge)
245	{
246	if(lastwas_falling)
247	{
248	newbin = bin;
249	break;
250	}
251	}
252	else lastwas_falling = 1; //last pixel was larger than this
253	lastcharge = charge;
254	}
255
256	//Sum the total charge of this sequence
257	seqcharge += charge;
258	seqaverage += data[bin].fTime*charge;
259	seqerror += data[bin].fTimedata[bin].fTimecharge;
260
261	//Check where to stop:
262	if(bin >= tempPt->fNDigit - 1) //out of range
263	break;
264	if(data[bin+1].fPad != data[bin].fPad) //new pad
265	break;
266	if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
267	break;
268
269	bin++;
270	}//end loop over sequence
271
272	//Calculate mean of sequence:
273	Int_t seqmean=0;
274	if(seqcharge)
275	seqmean = seqaverage/seqcharge;
276	else
277	{
278	LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Data")
279	<<"Error in data given to the cluster finder"<<ENDLOG;
280	seqmean = 1;
281	seqcharge = 1;
282	}
283
284	//Calculate mean in pad direction:
285	Int_t padmean = seqcharge*data[bin].fPad;
286	Int_t paderror = data[bin].fPad*padmean;
287
288	//Compare with results on previous pad:
289	for(UInt_t p=0; p<nprevious; p++)
290	{
291	//dont merge sequences on the same pad twice
292	if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
293
294	Int_t difference = seqmean - previousPt[p]->fMean;
295	if(difference < -fMatch) break;
296
297	if(difference <= fMatch) //There is a match here!!
298	{
299	AliClusterData *local = previousPt[p];
300
301	if(fDeconvPad)
302	{
303	if(seqcharge > local->fLastCharge)
304	{
305	if(local->fChargeFalling) //The previous pad was falling
306	{
307	break; //create a new cluster
308	}
309	}
310	else
311	local->fChargeFalling = 1;
312	local->fLastCharge = seqcharge;
313	}
314
315	//Don't create a new cluster, because we found a match
316	newcluster = kFALSE;
317
318	//Update cluster on current pad with the matching one:
319	local->fTotalCharge += seqcharge;
320	local->fPad += padmean;
321	local->fPad2 += paderror;
322	local->fTime += seqaverage;
323	local->fTime2 += seqerror;
324	local->fMean = seqmean;
325	local->fFlags++; //means we have more than one pad
326	local->fLastMergedPad = data[bin].fPad;
327
328	currentPt[ncurrent] = local;
329	ncurrent++;
330
331	break;
332	} //Checking for match at previous pad
333	} //Loop over results on previous pad.
334
335	if(newcluster)
336	{
337	//Start a new cluster. Add it to the clusterlist, and update
338	//the list of pointers to clusters in current pad.
339	//current pad will be previous pad on next pad.
340
341	//Add to the clusterlist:
342	AliClusterData *tmp = &clusterlist[ntotal];
343	tmp->fTotalCharge = seqcharge;
344	tmp->fPad = padmean;
345	tmp->fPad2 = paderror;
346	tmp->fTime = seqaverage;
347	tmp->fTime2 = seqerror;
348	tmp->fMean = seqmean;
349	tmp->fFlags = 0; //flags for single pad clusters
350	tmp->fLastMergedPad = data[bin].fPad;
351
352	if(fDeconvPad)
353	{
354	tmp->fChargeFalling = 0;
355	tmp->fLastCharge = seqcharge;
356	}
357
358	//Update list of pointers to previous pad:
359	currentPt[ncurrent] = &clusterlist[ntotal];
360	ntotal++;
361	ncurrent++;
362	}
363
364	if(fDeconvTime)
365	if(newbin >= 0) goto redo;
366	}//Loop over digits on this padrow
367
368	WriteClusters(ntotal,clusterlist);
369	}
370
371	void AliHLTTPCClustFinderNew::WriteClusters(Int_t nclusters,AliClusterData *list)
372	{
373	//write cluster to output pointer
374	Int_t thisrow,thissector;
375	UInt_t counter = fNClusters;
376
377	for(int j=0; j<nclusters; j++)
378	{
379	if(!list[j].fFlags) continue; //discard single pad clusters
380	if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
381
382	Float_t xyz[3];
383	Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
384	Float_t fpad2=fXYErr*fXYErr; //fixed given error
385	Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
386	Float_t ftime2=fZErr*fZErr; //fixed given error
387
388	if(fCalcerr) { //calc the errors, otherwice take the fixed error
389	Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
390	UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
391	Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
392	sy2/=q2;
393	if(sy2 < 0) {
394	LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
395	<<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
396	continue;
397	} else {
398	if(!fRawSP){
399	fpad2 = (sy2 + 1./12)AliHLTTPCTransform::GetPadPitchWidth(patch)AliHLTTPCTransform::GetPadPitchWidth(patch);
400	if(sy2 != 0){
401	fpad2*=0.108; //constants are from offline studies
402	if(patch<2)
403	fpad2*=2.07;
404	}
405	} else fpad2=sy2; //take the width not the error
406	}
407	Float_t sz2=list[j].fTime2list[j].fTotalCharge - list[j].fTimelist[j].fTime;
408	sz2/=q2;
409	if(sz2 < 0){
410	LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
411	<<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
412	continue;
413	} else {
414	if(!fRawSP){
415	ftime2 = (sz2 + 1./12)AliHLTTPCTransform::GetZWidth()AliHLTTPCTransform::GetZWidth();
416	if(sz2 != 0) {
417	ftime2 *= 0.169; //constants are from offline studies
418	if(patch<2)
419	ftime2 *= 1.77;
420	}
421	} else ftime2=sz2; //take the width, not the error
422	}
423	}
424	if(fStdout==kTRUE)
425	cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
426
427	if(!fRawSP){
428	AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
429	AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
430
431	if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
432	<<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
433	if(fNClusters >= fMaxNClusters)
434	{
435	LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
436	<<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
437	return;
438	}
439
440	fSpacePointData[counter].fX = xyz[0];
441	fSpacePointData[counter].fY = xyz[1];
442	fSpacePointData[counter].fZ = xyz[2];
443
444	} else {
445	fSpacePointData[counter].fX = fCurrentRow;
446	fSpacePointData[counter].fY = fpad;
447	fSpacePointData[counter].fZ = ftime;
448	}
449
450	fSpacePointData[counter].fCharge = list[j].fTotalCharge;
451	fSpacePointData[counter].fPadRow = fCurrentRow;
452	fSpacePointData[counter].fSigmaY2 = fpad2;
453	fSpacePointData[counter].fSigmaZ2 = ftime2;
454
455	Int_t patch=fCurrentPatch;
456	if(patch==-1) patch=0; //never store negative patch number
457	fSpacePointData[counter].fID = counter
458	+((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
459	#ifdef do_mc
460	Int_t trackID[3];
461	GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
462
463	fSpacePointData[counter].fTrackID[0] = trackID[0];
464	fSpacePointData[counter].fTrackID[1] = trackID[1];
465	fSpacePointData[counter].fTrackID[2] = trackID[2];
466
467	//cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
468	#endif
469
470	fNClusters++;
471	counter++;
472	}
473	}
474
475	#ifdef do_mc
476	void AliHLTTPCClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
477	{
478	//get mc id
479	AliHLTTPCDigitRowData rowPt = (AliHLTTPCDigitRowData)fDigitRowData;
480
481	trackID[0]=trackID[1]=trackID[2]=-2;
482	//cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
483	for(Int_t i=fFirstRow; i<=fLastRow; i++){
484	if(rowPt->fRow < (UInt_t)fCurrentRow){
485	AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
486	continue;
487	}
488	AliHLTTPCDigitData digPt = (AliHLTTPCDigitData)rowPt->fDigitData;
489	for(UInt_t j=0; j<rowPt->fNDigit; j++){
490	Int_t cpad = digPt[j].fPad;
491	Int_t ctime = digPt[j].fTime;
492	if(cpad != pad) continue;
493	if(ctime != time) continue;
494
495	trackID[0] = digPt[j].fTrackID[0];
496	trackID[1] = digPt[j].fTrackID[1];
497	trackID[2] = digPt[j].fTrackID[2];
498
499	//cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
500	break;
501	}
502	break;
503	}
504	}
505	#endif