2 // Original: AliL3ClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
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>
8 //*-- Copyright © ALICE HLT Group
11 #include "AliHLTTPCDigitReader.h"
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"
19 #include "AliHLTTPCPad.h"
25 /** \class AliHLTTPCClusterFinder
27 //_____________________________________________________________
28 // AliHLTTPCClusterFinder
30 // The current cluster finder for HLT
33 // The cluster finder is initialized with the Init function,
34 // providing the slice and patch information to work on.
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.
45 // There are several setters which control the behaviour:
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
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);
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();
96 ClassImp(AliHLTTPCClusterFinder)
98 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
116 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
119 fThreshold(src.fThreshold),
122 fDeconvPad(src.fDeconvPad),
123 fDeconvTime(src.fDeconvTime),
124 fStdout(src.fStdout),
125 fCalcerr(src.fCalcerr),
127 fFirstRow(src.fFirstRow),
128 fLastRow(src.fLastRow),
129 fDigitReader(src.fDigitReader)
133 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
136 fThreshold=src.fThreshold;
139 fDeconvPad=src.fDeconvPad;
140 fDeconvTime=src.fDeconvTime;
142 fCalcerr=src.fCalcerr;
144 fFirstRow=src.fFirstRow;
145 fLastRow=src.fLastRow;
146 fDigitReader=src.fDigitReader;
150 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
155 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
159 fMaxNClusters = nmaxpoints;
160 fCurrentSlice = slice;
161 fCurrentPatch = patch;
162 fFirstRow = firstrow;
166 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
170 fMaxNClusters = nmaxpoints;
171 fCurrentSlice = slice;
172 fCurrentPatch = patch;
173 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
174 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
177 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
179 //set pointer to output
180 fSpacePointData = pt;
183 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
185 fPtr = (UChar_t*)ptr;
189 void AliHLTTPCClusterFinder::ProcessDigits()
191 bool readValue = true;
194 UShort_t time=0,newTime=0;
195 UInt_t pad=0,newPad=0;
196 AliHLTTPCSignal_t charge=0;
200 // initialize block for reading packed data
201 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow);
202 readValue = fDigitReader->Next();
204 if (!readValue)return;
206 pad = fDigitReader->GetPad();
207 time = fDigitReader->GetTime();
208 fCurrentRow = fDigitReader->GetRow();
210 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
211 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
213 fCurrentRow += rowOffset;
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
220 AliClusterData **currentPt; //List of pointers to the current pad
221 AliClusterData **previousPt; //List of pointers to the previous pad
224 UInt_t nprevious=0,ncurrent=0,ntotal=0;
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
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;
238 while ( readValue ){ // Reads through all digits in block
244 if(currentPt == pad2){
252 nprevious = ncurrent;
254 if(pad != lastpad+1){
255 //this happens if there is a pad with no signal.
256 nprevious = ncurrent = 0;
261 Bool_t newcluster = kTRUE;
262 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
263 UInt_t lastcharge=0,lastwas_falling=0;
268 redo: //This is a goto.
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
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) {
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);
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();
300 time=pCurrentPad->GetCurrentPosition();
301 if (time>pCurrentPad->GetSize()) {
302 HLTError("invalid time bin for pad");
309 charge = pCurrentPad->GetCorrectedData();
311 charge = fDigitReader->GetSignal();
313 //HLTDebug("get next charge value: position %d charge %d", time, charge);
317 if (fDigitReader->GetRow() == 90){
318 ///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
319 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
323 if(time >= AliHLTTPCTransform::GetNTimeBins()){
324 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
329 //Get the current ADC-value
332 //Check if the last pixel in the sequence is smaller than this
333 if(charge > lastcharge){
339 else lastwas_falling = 1; //last pixel was larger than this
343 //Sum the total charge of this sequence
345 seqaverage += time*charge;
346 seqerror += time*time*charge;
350 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
351 pCurrentPad->StopEvent();
352 pCurrentPad->ResetHistory();
354 newPad = fDigitReader->GetPad();
355 newTime = fDigitReader->GetTime();
356 newRow = fDigitReader->GetRow() + rowOffset;
361 newPad=pCurrentPad->GetPadNumber();
362 newTime=pCurrentPad->GetCurrentPosition();
363 newRow=pCurrentPad->GetRowNumber();
365 readValue = fDigitReader->Next();
366 //Check where to stop:
367 if(!readValue) break; //No more value
369 newPad = fDigitReader->GetPad();
370 newTime = fDigitReader->GetTime();
371 newRow = fDigitReader->GetRow() + rowOffset;
374 if(newPad != pad)break; //new pad
375 if(newTime != time+1) break; //end of sequence
377 // pad = newpad; is equal
380 }//end loop over sequence
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);
385 // with active zero suppression zero values are possible
389 //Calculate mean of sequence:
392 seqmean = seqaverage/seqcharge;
394 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
395 <<"Error in data given to the cluster finder"<<ENDLOG;
400 //Calculate mean in pad direction:
401 Int_t padmean = seqcharge*pad;
402 Int_t paderror = pad*padmean;
405 //Compare with results on previous pad:
406 for(UInt_t p=0; p<nprevious; p++){
408 //dont merge sequences on the same pad twice
409 if(previousPt[p]->fLastMergedPad==pad) continue;
411 Int_t difference = seqmean - previousPt[p]->fMean;
412 if(difference < -fMatch) break;
414 if(difference <= fMatch){ //There is a match here!!
415 AliClusterData *local = previousPt[p];
418 if(seqcharge > local->fLastCharge){
419 if(local->fChargeFalling){ //The previous pad was falling
420 break; //create a new cluster
423 else local->fChargeFalling = 1;
424 local->fLastCharge = seqcharge;
427 //Don't create a new cluster, because we found a match
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;
440 currentPt[ncurrent] = local;
444 } //Checking for match at previous pad
445 } //Loop over results on previous pad.
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.
453 //Add to the clusterlist:
454 AliClusterData *tmp = &clusterlist[ntotal];
455 tmp->fTotalCharge = seqcharge;
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;
465 tmp->fChargeFalling = 0;
466 tmp->fLastCharge = seqcharge;
469 //Update list of pointers to previous pad:
470 currentPt[ncurrent] = &clusterlist[ntotal];
476 if(newbin >= 0) goto redo;
478 // to prevent endless loop
479 if(time >= AliHLTTPCTransform::GetNTimeBins()){
480 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
485 if(!readValue) break; //No more value
487 if(fCurrentRow != newRow){
488 WriteClusters(ntotal,clusterlist);
498 fCurrentRow = newRow;
504 } // END while(readValue)
506 if (pCurrentPad) pCurrentPad->StopEvent();
508 WriteClusters(ntotal,clusterlist);
510 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
514 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
516 //write cluster to output pointer
517 Int_t thisrow,thissector;
518 UInt_t counter = fNClusters;
520 for(int j=0; j<nclusters; j++)
522 if(!list[j].fFlags) continue; //discard single pad clusters
523 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
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
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;
542 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
543 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
547 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
549 fpad2*=0.108; //constants are from offline studies
553 } else fpad2=sy2; //take the width not the error
555 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
558 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
559 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
563 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
565 ftime2 *= 0.169; //constants are from offline studies
569 } else ftime2=sz2; //take the width, not the error
573 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
576 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
577 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
579 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
580 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
581 if(fNClusters >= fMaxNClusters)
583 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
584 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
588 fSpacePointData[counter].fX = xyz[0];
589 fSpacePointData[counter].fY = xyz[1];
590 fSpacePointData[counter].fZ = xyz[2];
593 fSpacePointData[counter].fX = fCurrentRow;
594 fSpacePointData[counter].fY = fpad;
595 fSpacePointData[counter].fZ = ftime;
598 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
599 fSpacePointData[counter].fPadRow = fCurrentRow;
600 fSpacePointData[counter].fSigmaY2 = fpad2;
601 fSpacePointData[counter].fSigmaZ2 = ftime2;
603 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
604 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
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
613 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
615 fSpacePointData[counter].fTrackID[0] = trackID[0];
616 fSpacePointData[counter].fTrackID[1] = trackID[1];
617 fSpacePointData[counter].fTrackID[2] = trackID[2];
619 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
627 // STILL TO FIX ----------------------------------------------------------------------------
630 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
633 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
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);
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;
649 trackID[0] = digPt[j].fTrackID[0];
650 trackID[1] = digPt[j].fTrackID[1];
651 trackID[2] = digPt[j].fTrackID[2];
653 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;