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

// @(#) $Id$
// Original: AliL3ClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp 

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, 
//         Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
//         Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>

//*-- Copyright &copy ALICE HLT Group


#include "AliHLTTPCDigitReader.h"
#include "AliHLTTPCRootTypes.h"
#include "AliHLTTPCLogging.h"
#include "AliHLTTPCClusterFinder.h"
#include "AliHLTTPCDigitData.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCMemHandler.h"
#include "AliHLTTPCPad.h"

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

/** \class AliHLTTPCClusterFinder
<pre>
//_____________________________________________________________
// AliHLTTPCClusterFinder
//
// 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 provided by the AliHLTTPCDigitReader class,
// using the init() funktion, and the next() funktion in order 
// to get the next bin. Either packed or unpacked data can be
// processed, dependent if one uses AliHLTTPCDigitReaderPacked 
// class or AliHLTTPCDigitReaderUnpacked class in the 
// Clusterfinder Component.
// 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);
//     AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
//     cf->SetMatchWidth(2);
//     cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
//     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(iter->fPtr, iter->fSize ); //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(AliHLTTPCClusterFinder)

AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
  :
  fMatch(1),
  fThreshold(10),
  fXYErr(0.2),
  fZErr(0.3),
  fDeconvPad(kTRUE),
  fDeconvTime(kTRUE),
  fStdout(kFALSE),
  fCalcerr(kTRUE),
  fRawSP(kFALSE),
  fFirstRow(0),
  fLastRow(0),
  fDigitReader(NULL)
{
  //constructor
}

AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
  :
  fMatch(src.fMatch),
  fThreshold(src.fThreshold),
  fXYErr(src.fXYErr),
  fZErr(src.fZErr),
  fDeconvPad(src.fDeconvPad),
  fDeconvTime(src.fDeconvTime),
  fStdout(src.fStdout),
  fCalcerr(src.fCalcerr),
  fRawSP(src.fRawSP),
  fFirstRow(src.fFirstRow),
  fLastRow(src.fLastRow),
  fDigitReader(src.fDigitReader)
{
}

AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
{
  fMatch=src.fMatch;
  fThreshold=src.fThreshold;
  fXYErr=src.fXYErr;
  fZErr=src.fZErr;
  fDeconvPad=src.fDeconvPad;
  fDeconvTime=src.fDeconvTime;
  fStdout=src.fStdout;
  fCalcerr=src.fCalcerr;
  fRawSP=src.fRawSP;
  fFirstRow=src.fFirstRow;
  fLastRow=src.fLastRow;
  fDigitReader=src.fDigitReader;
  return (*this);
}

AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
{
  //destructor
}
 
void AliHLTTPCClusterFinder::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 AliHLTTPCClusterFinder::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 AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
{
  //set pointer to output
  fSpacePointData = pt;
}

void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
  //set input pointer
  fPtr = (UChar_t*)ptr;
  fSize = size;
}

void AliHLTTPCClusterFinder::ProcessDigits()
{
  bool readValue = true;
  Int_t newRow = 0;    
  Int_t rowOffset = 0;
  UShort_t time=0,newTime=0;
  UInt_t pad=0,newPad=0;
  AliHLTTPCSignal_t charge=0;

  fNClusters = 0;

  // initialize block for reading packed data
  fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow);
  readValue = fDigitReader->Next();

  if (!readValue)return;
 
  pad = fDigitReader->GetPad();
  time = fDigitReader->GetTime();
  fCurrentRow = fDigitReader->GetRow();

  if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
    rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );

  fCurrentRow += rowOffset;

  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;

  /* quick implementation of baseline calculation and zero suppression
     open a pad object for each pad and delete it after processing.
     later a list of pad objects with base line history can be used
     The whole thing only works if we really get unprocessed raw data, if
     the data is already zero suppressed, there might be gaps in the time
     bins.
   */
  Int_t gatingGridOffset=50;
  AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
  // just to make later conversion to a list of objects easier
  AliHLTTPCPad* pCurrentPad=&baseline;

  while ( readValue ){   // Reads through all digits in block

    if(pad != 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(pad != lastpad+1){
	//this happens if there is a pad with no signal.
	nprevious = ncurrent = 0;
      }
      lastpad = pad;
    }

    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 time bins of current pad
      // read all the values for one pad at once to calculate the base line
      if (pCurrentPad) {
	if (!pCurrentPad->IsStarted()) {
	  //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
	  pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
	  if ((pCurrentPad->StartEvent())>=0) {
	    do {
	      if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
	      if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
	      pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
	      //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
	    } while ((readValue = fDigitReader->Next())!=0);
	  }
	  pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
	  if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
	    HLTDebug("no data available after zero suppression");
	    pCurrentPad->StopEvent();
	    pCurrentPad->ResetHistory();
	    break;
	  }
	  time=pCurrentPad->GetCurrentPosition();
	  if (time>pCurrentPad->GetSize()) {
	    HLTError("invalid time bin for pad");
	    break;
	  }
	}
      }

      if (pCurrentPad) {
	charge = pCurrentPad->GetCorrectedData();
      } else {
      charge = fDigitReader->GetSignal();
      }
      //HLTDebug("get next charge value: position %d charge %d", time, charge);


      // CHARGE DEBUG
      if (fDigitReader->GetRow() == 90){
/////	  LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime() 
	  //					   << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;

      }

      if(time >= AliHLTTPCTransform::GetNTimeBins()){
	HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
	break;
      }


      //Get the current ADC-value
      if(fDeconvTime){

	//Check if the last pixel in the sequence is smaller than this
	if(charge > lastcharge){
	  if(lastwas_falling){
	    newbin = 1;
	    break;
	  }
	}
	else lastwas_falling = 1; //last pixel was larger than this
	lastcharge = charge;
      }
	  
      //Sum the total charge of this sequence
      seqcharge += charge;
      seqaverage += time*charge;
      seqerror += time*time*charge;
      
      if (pCurrentPad) {
	
	if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
	  pCurrentPad->StopEvent();
	  pCurrentPad->ResetHistory();
	  if(readValue) {
	    newPad = fDigitReader->GetPad();
	    newTime = fDigitReader->GetTime();
	    newRow = fDigitReader->GetRow() + rowOffset;
	  }
	  break;
	}

	newPad=pCurrentPad->GetPadNumber();
	newTime=pCurrentPad->GetCurrentPosition();
	newRow=pCurrentPad->GetRowNumber();
      } else {
      readValue = fDigitReader->Next();
      //Check where to stop:
      if(!readValue) break; //No more value

      newPad = fDigitReader->GetPad();
      newTime = fDigitReader->GetTime();
      newRow = fDigitReader->GetRow() + rowOffset;
      }

      if(newPad != pad)break; //new pad
      if(newTime != time+1) break; //end of sequence

      // pad = newpad;    is equal
      time = newTime;

    }//end loop over sequence

    //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
    //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
    if (seqcharge<=0) {
      // with active zero suppression zero values are possible
      continue;
    }

    //Calculate mean of sequence:
    Int_t seqmean=0;
    if(seqcharge)
      seqmean = seqaverage/seqcharge;
    else{
      LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
	<<"Error in data given to the cluster finder"<<ENDLOG;
      seqmean = 1;
      seqcharge = 1;
    }

    //Calculate mean in pad direction:
    Int_t padmean = seqcharge*pad;
    Int_t paderror = pad*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==pad) 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 = pad;

	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 = pad;

      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;
  
    // to prevent endless loop  
    if(time >= AliHLTTPCTransform::GetNTimeBins()){
      HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
      break;
    }


    if(!readValue) break; //No more value

    if(fCurrentRow != newRow){
      WriteClusters(ntotal,clusterlist);

      lastpad = 123456789;

      currentPt = pad2;
      previousPt = pad1;
      nprevious=0;
      ncurrent=0;
      ntotal=0;
      
      fCurrentRow = newRow;
    }

    pad = newPad;
    time = newTime;
  
  } // END while(readValue)

  if (pCurrentPad) pCurrentPad->StopEvent();

  WriteClusters(ntotal,clusterlist);

  HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);

} // ENDEND

void AliHLTTPCClusterFinder::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,"AliHLTTPCClusterFinder::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,"AliHLTTPCClusterFinder::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;

      fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
      fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay

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

// STILL TO FIX  ----------------------------------------------------------------------------

#ifdef do_mc
void AliHLTTPCClusterFinder::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
a38a7850	1	// @(#) $Id$
c2f37813	2	// Original: AliL3ClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
a38a7850	3
	4	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>,
	5	// Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
	6	// Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
	7
	8	//*-- Copyright &copy ALICE HLT Group
	9
	10
	11	#include "AliHLTTPCDigitReader.h"
a38a7850	12	#include "AliHLTTPCRootTypes.h"
	13	#include "AliHLTTPCLogging.h"
	14	#include "AliHLTTPCClusterFinder.h"
	15	#include "AliHLTTPCDigitData.h"
	16	#include "AliHLTTPCTransform.h"
	17	#include "AliHLTTPCSpacePointData.h"
	18	#include "AliHLTTPCMemHandler.h"
46b33a24	19	#include "AliHLTTPCPad.h"
a38a7850	20
	21	#if __GNUC__ >= 3
	22	using namespace std;
	23	#endif
	24
	25	/** \class AliHLTTPCClusterFinder
	26	<pre>
	27	//_____________________________________________________________
	28	// AliHLTTPCClusterFinder
	29	//
	30	// The current cluster finder for HLT
	31	// (Based on STAR L3)
	32	//
	33	// The cluster finder is initialized with the Init function,
	34	// providing the slice and patch information to work on.
	35	//
	36	// The input is a provided by the AliHLTTPCDigitReader class,
	37	// using the init() funktion, and the next() funktion in order
	38	// to get the next bin. Either packed or unpacked data can be
	39	// processed, dependent if one uses AliHLTTPCDigitReaderPacked
	40	// class or AliHLTTPCDigitReaderUnpacked class in the
	41	// Clusterfinder Component.
	42	// The resulting space points will be in the
	43	// array given by the SetOutputArray function.
	44	//
	45	// There are several setters which control the behaviour:
	46	//
	47	// - SetXYError(Float_t): set fixed error in XY direction
	48	// - SetZError(Float_t): set fixed error in Z direction
	49	// (used if errors are not calculated)
	50	// - SetDeconv(Bool_t): switch on/off deconvolution
	51	// - SetThreshold(UInt_t): set charge threshold for cluster
	52	// - SetMatchWidth(UInt_t): set the match distance in
	53	// time for sequences to be merged
	54	// - SetSTDOutput(Bool_t): switch on/off output about found clusters
	55	// - SetCalcErr(Bool_t): switch on/off calculation of
	56	// space point errors (or widths in raw system)
	57	// - SetRawSP(Bool_t): switch on/off convertion to raw system
	58	//
	59	//
	60	// Example Usage:
	61	//
	62	// AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
	63	// file->SetAliInput(digitfile); //give some input file
	64	// for(int slice=0; slice<=35; slice++){
	65	// for(int patch=0; pat<6; pat++){
	66	// file->Init(slice,patch);
	67	// UInt_t ndigits=0;
	68	// UInt_t maxclusters=100000;
	69	// UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
	70	// AliHLTTPCSpacePointData points = (AliHLTTPCSpacePointData)memory->Allocate(pointsize);
	71	// AliHLTTPCDigitRowData digits = (AliHLTTPCDigitRowData)file->AliAltroDigits2Memory(ndigits,event);
	72	// AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
	73	// cf->SetMatchWidth(2);
	74	// cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
	75	// cf->SetSTDOutput(kTRUE); //Some output to standard IO
	76	// cf->SetRawSP(kFALSE); //Convert space points to local system
	77	// cf->SetThreshold(5); //Threshold of cluster charge
	78	// cf->SetDeconv(kTRUE); //Deconv in pad and time direction
	79	// cf->SetCalcErr(kTRUE); //Calculate the errors of the spacepoints
	80	// cf->SetOutputArray(points); //Move the spacepoints to the array
	81	// cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
	82	// cf->ProcessDigits(); //process the rows given by init
	83	// Int_t npoints = cf->GetNumberOfClusters();
84	// AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
85	// out->SetBinaryOutput(fname);
86	// out->Memory2Binary(npoints,points); //store the spacepoints
87	// out->CloseBinaryOutput();
88	// delete out;
89	// file->free();
90	// delete cf;
91	// }
92	// }
93	</pre>
94	*/
95
96	ClassImp(AliHLTTPCClusterFinder)
97
98	AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
46b33a24	99	:
	100	fMatch(1),
	101	fThreshold(10),
	102	fXYErr(0.2),
	103	fZErr(0.3),
	104	fDeconvPad(kTRUE),
	105	fDeconvTime(kTRUE),
	106	fStdout(kFALSE),
	107	fCalcerr(kTRUE),
	108	fRawSP(kFALSE),
	109	fFirstRow(0),
	110	fLastRow(0),
	111	fDigitReader(NULL)
a38a7850	112	{
a38a7850	113	//constructor
46b33a24	114	}
a38a7850	115
46b33a24	116	AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
	117	:
	118	fMatch(src.fMatch),
	119	fThreshold(src.fThreshold),
	120	fXYErr(src.fXYErr),
	121	fZErr(src.fZErr),
	122	fDeconvPad(src.fDeconvPad),
	123	fDeconvTime(src.fDeconvTime),
	124	fStdout(src.fStdout),
	125	fCalcerr(src.fCalcerr),
	126	fRawSP(src.fRawSP),
	127	fFirstRow(src.fFirstRow),
	128	fLastRow(src.fLastRow),
	129	fDigitReader(src.fDigitReader)
	130	{
	131	}
	132
	133	AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
	134	{
	135	fMatch=src.fMatch;
	136	fThreshold=src.fThreshold;
	137	fXYErr=src.fXYErr;
	138	fZErr=src.fZErr;
	139	fDeconvPad=src.fDeconvPad;
	140	fDeconvTime=src.fDeconvTime;
	141	fStdout=src.fStdout;
	142	fCalcerr=src.fCalcerr;
	143	fRawSP=src.fRawSP;
	144	fFirstRow=src.fFirstRow;
	145	fLastRow=src.fLastRow;
	146	fDigitReader=src.fDigitReader;
	147	return (*this);
a38a7850	148	}
	149
	150	AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
	151	{
	152	//destructor
	153	}
	154
	155	void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
	156	{
	157	//init slice
	158	fNClusters = 0;
	159	fMaxNClusters = nmaxpoints;
	160	fCurrentSlice = slice;
	161	fCurrentPatch = patch;
	162	fFirstRow = firstrow;
	163	fLastRow = lastrow;
	164	}
	165
	166	void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
	167	{
	168	//init slice
	169	fNClusters = 0;
	170	fMaxNClusters = nmaxpoints;
	171	fCurrentSlice = slice;
	172	fCurrentPatch = patch;
	173	fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
	174	fLastRow=AliHLTTPCTransform::GetLastRow(patch);
	175	}
	176
	177	void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
	178	{
	179	//set pointer to output
	180	fSpacePointData = pt;
	181	}
	182
	183	void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
	184	//set input pointer
	185	fPtr = (UChar_t*)ptr;
	186	fSize = size;
	187	}
	188
	189	void AliHLTTPCClusterFinder::ProcessDigits()
	190	{
	191	bool readValue = true;
	192	Int_t newRow = 0;
	193	Int_t rowOffset = 0;
46b33a24	194	UShort_t time=0,newTime=0;
	195	UInt_t pad=0,newPad=0;
	196	AliHLTTPCSignal_t charge=0;
a38a7850	197
	198	fNClusters = 0;
	199
	200	// initialize block for reading packed data
	201	fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow);
	202	readValue = fDigitReader->Next();
	203
	204	if (!readValue)return;
	205
	206	pad = fDigitReader->GetPad();
	207	time = fDigitReader->GetTime();
	208	fCurrentRow = fDigitReader->GetRow();
	209
	210	if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
	211	rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
	212
	213	fCurrentRow += rowOffset;
	214
	215	UInt_t lastpad = 123456789;
	216	AliClusterData *pad1[5000]; //2 lists for internal memory=2pads
	217	AliClusterData *pad2[5000]; //2 lists for internal memory=2pads
	218	AliClusterData clusterlist[10000]; //Clusterlist
	219
	220	AliClusterData **currentPt; //List of pointers to the current pad
	221	AliClusterData **previousPt; //List of pointers to the previous pad
	222	currentPt = pad2;
	223	previousPt = pad1;
	224	UInt_t nprevious=0,ncurrent=0,ntotal=0;
	225
46b33a24	226	/* quick implementation of baseline calculation and zero suppression
	227	open a pad object for each pad and delete it after processing.
	228	later a list of pad objects with base line history can be used
	229	The whole thing only works if we really get unprocessed raw data, if
	230	the data is already zero suppressed, there might be gaps in the time
	231	bins.
	232	*/
	233	Int_t gatingGridOffset=50;
	234	AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
	235	// just to make later conversion to a list of objects easier
	236	AliHLTTPCPad* pCurrentPad=&baseline;
	237
a38a7850	238	while ( readValue ){ // Reads through all digits in block
	239
	240	if(pad != lastpad){
	241	//This is a new pad
	242
	243	//Switch the lists:
	244	if(currentPt == pad2){
	245	currentPt = pad1;
	246	previousPt = pad2;
	247	}
	248	else {
	249	currentPt = pad2;
	250	previousPt = pad1;
	251	}
	252	nprevious = ncurrent;
	253	ncurrent = 0;
	254	if(pad != lastpad+1){
	255	//this happens if there is a pad with no signal.
	256	nprevious = ncurrent = 0;
	257	}
	258	lastpad = pad;
	259	}
	260
	261	Bool_t newcluster = kTRUE;
	262	UInt_t seqcharge=0,seqaverage=0,seqerror=0;
	263	UInt_t lastcharge=0,lastwas_falling=0;
	264	Int_t newbin=-1;
	265
	266
	267	if(fDeconvTime){
	268	redo: //This is a goto.
	269
	270	if(newbin > -1){
	271	//bin = newbin;
	272	newbin = -1;
	273	}
	274
	275	lastcharge=0;
	276	lastwas_falling = 0;
	277	}
	278
46b33a24	279	while(1){ //Loop over time bins of current pad
	280	// read all the values for one pad at once to calculate the base line
	281	if (pCurrentPad) {
	282	if (!pCurrentPad->IsStarted()) {
	283	//HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
	284	pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
	285	if ((pCurrentPad->StartEvent())>=0) {
	286	do {
	287	if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
	288	if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
	289	pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
	290	//HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
	291	} while ((readValue = fDigitReader->Next())!=0);
	292	}
	293	pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
	294	if (pCurrentPad->Next(kTRUE/do zero suppression/)==0) {
	295	HLTDebug("no data available after zero suppression");
	296	pCurrentPad->StopEvent();
	297	pCurrentPad->ResetHistory();
	298	break;
	299	}
	300	time=pCurrentPad->GetCurrentPosition();
	301	if (time>pCurrentPad->GetSize()) {
	302	HLTError("invalid time bin for pad");
	303	break;
	304	}
	305	}
	306	}
a38a7850	307
46b33a24	308	if (pCurrentPad) {
	309	charge = pCurrentPad->GetCorrectedData();
	310	} else {
a38a7850	311	charge = fDigitReader->GetSignal();
46b33a24	312	}
46b33a24	313	//HLTDebug("get next charge value: position %d charge %d", time, charge);
a38a7850	314
738c049f	315
	316	// CHARGE DEBUG
	317	if (fDigitReader->GetRow() == 90){
	318	///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
	319	// << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
	320
	321	}
	322
a38a7850	323	if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24	324	HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850	325	break;
a38a7850	326	}
46b33a24	327
46b33a24	328
a38a7850	329	//Get the current ADC-value
	330	if(fDeconvTime){
	331
	332	//Check if the last pixel in the sequence is smaller than this
	333	if(charge > lastcharge){
	334	if(lastwas_falling){
	335	newbin = 1;
	336	break;
	337	}
	338	}
	339	else lastwas_falling = 1; //last pixel was larger than this
	340	lastcharge = charge;
	341	}
	342
	343	//Sum the total charge of this sequence
	344	seqcharge += charge;
	345	seqaverage += time*charge;
	346	seqerror += timetimecharge;
	347
46b33a24	348	if (pCurrentPad) {
	349
	350	if((pCurrentPad->Next(kTRUE/do zero suppression/))==0) {
	351	pCurrentPad->StopEvent();
	352	pCurrentPad->ResetHistory();
	353	if(readValue) {
	354	newPad = fDigitReader->GetPad();
	355	newTime = fDigitReader->GetTime();
	356	newRow = fDigitReader->GetRow() + rowOffset;
	357	}
	358	break;
	359	}
	360
	361	newPad=pCurrentPad->GetPadNumber();
	362	newTime=pCurrentPad->GetCurrentPosition();
	363	newRow=pCurrentPad->GetRowNumber();
	364	} else {
a38a7850	365	readValue = fDigitReader->Next();
a38a7850	366	//Check where to stop:
	367	if(!readValue) break; //No more value
	368
	369	newPad = fDigitReader->GetPad();
	370	newTime = fDigitReader->GetTime();
	371	newRow = fDigitReader->GetRow() + rowOffset;
46b33a24	372	}
a38a7850	373
	374	if(newPad != pad)break; //new pad
	375	if(newTime != time+1) break; //end of sequence
	376
	377	// pad = newpad; is equal
	378	time = newTime;
	379
	380	}//end loop over sequence
	381
46b33a24	382	//HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
	383	//HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
	384	if (seqcharge<=0) {
	385	// with active zero suppression zero values are possible
	386	continue;
	387	}
a38a7850	388
	389	//Calculate mean of sequence:
	390	Int_t seqmean=0;
	391	if(seqcharge)
	392	seqmean = seqaverage/seqcharge;
	393	else{
	394	LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
	395	<<"Error in data given to the cluster finder"<<ENDLOG;
	396	seqmean = 1;
	397	seqcharge = 1;
	398	}
	399
	400	//Calculate mean in pad direction:
	401	Int_t padmean = seqcharge*pad;
	402	Int_t paderror = pad*padmean;
	403
	404
	405	//Compare with results on previous pad:
	406	for(UInt_t p=0; p<nprevious; p++){
	407
	408	//dont merge sequences on the same pad twice
	409	if(previousPt[p]->fLastMergedPad==pad) continue;
	410
	411	Int_t difference = seqmean - previousPt[p]->fMean;
	412	if(difference < -fMatch) break;
	413
	414	if(difference <= fMatch){ //There is a match here!!
	415	AliClusterData *local = previousPt[p];
	416
	417	if(fDeconvPad){
	418	if(seqcharge > local->fLastCharge){
	419	if(local->fChargeFalling){ //The previous pad was falling
	420	break; //create a new cluster
	421	}
	422	}
	423	else local->fChargeFalling = 1;
	424	local->fLastCharge = seqcharge;
	425	}
	426
	427	//Don't create a new cluster, because we found a match
	428	newcluster = kFALSE;
	429
	430	//Update cluster on current pad with the matching one:
	431	local->fTotalCharge += seqcharge;
	432	local->fPad += padmean;
	433	local->fPad2 += paderror;
	434	local->fTime += seqaverage;
	435	local->fTime2 += seqerror;
	436	local->fMean = seqmean;
	437	local->fFlags++; //means we have more than one pad
	438	local->fLastMergedPad = pad;
	439
	440	currentPt[ncurrent] = local;
	441	ncurrent++;
	442
	443	break;
	444	} //Checking for match at previous pad
	445	} //Loop over results on previous pad.
	446
	447
	448	if(newcluster){
	449	//Start a new cluster. Add it to the clusterlist, and update
	450	//the list of pointers to clusters in current pad.
	451	//current pad will be previous pad on next pad.
452
453	//Add to the clusterlist:
454	AliClusterData *tmp = &clusterlist[ntotal];
455	tmp->fTotalCharge = seqcharge;
456	tmp->fPad = padmean;
457	tmp->fPad2 = paderror;
458	tmp->fTime = seqaverage;
459	tmp->fTime2 = seqerror;
460	tmp->fMean = seqmean;
461	tmp->fFlags = 0; //flags for single pad clusters
462	tmp->fLastMergedPad = pad;
463
464	if(fDeconvPad){
465	tmp->fChargeFalling = 0;
466	tmp->fLastCharge = seqcharge;
467	}
468
469	//Update list of pointers to previous pad:
470	currentPt[ncurrent] = &clusterlist[ntotal];
471	ntotal++;
472	ncurrent++;
473	}
474
475	if(fDeconvTime)
476	if(newbin >= 0) goto redo;
477
478	// to prevent endless loop
479	if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24	480	HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850	481	break;
	482	}
	483
	484
	485	if(!readValue) break; //No more value
	486
	487	if(fCurrentRow != newRow){
	488	WriteClusters(ntotal,clusterlist);
	489
	490	lastpad = 123456789;
	491
	492	currentPt = pad2;
	493	previousPt = pad1;
	494	nprevious=0;
	495	ncurrent=0;
	496	ntotal=0;
	497
	498	fCurrentRow = newRow;
	499	}
	500
	501	pad = newPad;
	502	time = newTime;
	503
	504	} // END while(readValue)
	505
46b33a24	506	if (pCurrentPad) pCurrentPad->StopEvent();
46b33a24	507
a38a7850	508	WriteClusters(ntotal,clusterlist);
a38a7850	509
46b33a24	510	HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850	511
	512	} // ENDEND
	513
	514	void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
	515	{
	516	//write cluster to output pointer
	517	Int_t thisrow,thissector;
	518	UInt_t counter = fNClusters;
	519
	520	for(int j=0; j<nclusters; j++)
	521	{
	522	if(!list[j].fFlags) continue; //discard single pad clusters
	523	if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
	524
	525	Float_t xyz[3];
	526	Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
	527	Float_t fpad2=fXYErr*fXYErr; //fixed given error
	528	Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
	529	Float_t ftime2=fZErr*fZErr; //fixed given error
	530
738c049f	531
	532
	533
	534
	535
a38a7850	536	if(fCalcerr) { //calc the errors, otherwice take the fixed error
	537	Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
	538	UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
	539	Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
	540	sy2/=q2;
	541	if(sy2 < 0) {
	542	LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
	543	<<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	544	continue;
	545	} else {
	546	if(!fRawSP){
	547	fpad2 = (sy2 + 1./12)AliHLTTPCTransform::GetPadPitchWidth(patch)AliHLTTPCTransform::GetPadPitchWidth(patch);
	548	if(sy2 != 0){
	549	fpad2*=0.108; //constants are from offline studies
	550	if(patch<2)
	551	fpad2*=2.07;
	552	}
	553	} else fpad2=sy2; //take the width not the error
	554	}
	555	Float_t sz2=list[j].fTime2list[j].fTotalCharge - list[j].fTimelist[j].fTime;
	556	sz2/=q2;
	557	if(sz2 < 0){
	558	LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
	559	<<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
	560	continue;
	561	} else {
	562	if(!fRawSP){
	563	ftime2 = (sz2 + 1./12)AliHLTTPCTransform::GetZWidth()AliHLTTPCTransform::GetZWidth();
	564	if(sz2 != 0) {
	565	ftime2 *= 0.169; //constants are from offline studies
	566	if(patch<2)
	567	ftime2 *= 1.77;
	568	}
	569	} else ftime2=sz2; //take the width, not the error
	570	}
	571	}
	572	if(fStdout==kTRUE)
	573	cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
	574
	575	if(!fRawSP){
	576	AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
	577	AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
	578
	579	if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
	580	<<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
	581	if(fNClusters >= fMaxNClusters)
	582	{
	583	LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
	584	<<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
	585	return;
	586	}
	587
	588	fSpacePointData[counter].fX = xyz[0];
	589	fSpacePointData[counter].fY = xyz[1];
	590	fSpacePointData[counter].fZ = xyz[2];
	591
	592	} else {
	593	fSpacePointData[counter].fX = fCurrentRow;
	594	fSpacePointData[counter].fY = fpad;
	595	fSpacePointData[counter].fZ = ftime;
	596	}
	597
	598	fSpacePointData[counter].fCharge = list[j].fTotalCharge;
	599	fSpacePointData[counter].fPadRow = fCurrentRow;
600	fSpacePointData[counter].fSigmaY2 = fpad2;
601	fSpacePointData[counter].fSigmaZ2 = ftime2;
602
44be0fde	603	fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
738c049f	604	fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
44be0fde	605
a38a7850	606	Int_t patch=fCurrentPatch;
	607	if(patch==-1) patch=0; //never store negative patch number
	608	fSpacePointData[counter].fID = counter
	609	+((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
	610
	611	#ifdef do_mc
	612	Int_t trackID[3];
	613	GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
	614
	615	fSpacePointData[counter].fTrackID[0] = trackID[0];
	616	fSpacePointData[counter].fTrackID[1] = trackID[1];
	617	fSpacePointData[counter].fTrackID[2] = trackID[2];
	618
	619	//cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
	620	#endif
	621
	622	fNClusters++;
	623	counter++;
	624	}
	625	}
	626
	627	// STILL TO FIX ----------------------------------------------------------------------------
	628
	629	#ifdef do_mc
	630	void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
	631	{
	632	//get mc id
	633	AliHLTTPCDigitRowData rowPt = (AliHLTTPCDigitRowData)fDigitRowData;
	634
	635	trackID[0]=trackID[1]=trackID[2]=-2;
	636	//cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
	637	for(Int_t i=fFirstRow; i<=fLastRow; i++){
	638	if(rowPt->fRow < (UInt_t)fCurrentRow){
	639	AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
	640	continue;
	641	}
	642	AliHLTTPCDigitData digPt = (AliHLTTPCDigitData)rowPt->fDigitData;
	643	for(UInt_t j=0; j<rowPt->fNDigit; j++){
	644	Int_t cpad = digPt[j].fPad;
	645	Int_t ctime = digPt[j].fTime;
	646	if(cpad != pad) continue;
	647	if(ctime != time) continue;
	648
	649	trackID[0] = digPt[j].fTrackID[0];
	650	trackID[1] = digPt[j].fTrackID[1];
	651	trackID[2] = digPt[j].fTrackID[2];
	652
	653	//cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
	654	break;
	655	}
	656	break;
	657	}
	658	}
	659	#endif