3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliHLTStandardIncludes.h"
7 #include "AliHLTRootTypes.h"
8 #include "AliHLTLogging.h"
9 #include "AliHLTClustFinderNew.h"
10 #include "AliHLTDigitData.h"
11 #include "AliHLTTransform.h"
12 #include "AliHLTSpacePointData.h"
13 #include "AliHLTMemHandler.h"
19 /** \class AliHLTClustFinderNew
21 //_____________________________________________________________
22 // AliHLTClustFinderNew
24 // The current cluster finder for HLT
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 AliHLTDigitRowData structure using the
30 // Read function. The resulting space points will be in the
31 // array given by the SetOutputArray function.
33 // There are several setters which control the behaviour:
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
50 // AliHLTFileHandler *file = new AliHLTFileHandler();
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);
56 // UInt_t maxclusters=100000;
57 // UInt_t pointsize = maxclusters*sizeof(AliHLTSpacePointData);
58 // AliHLTSpacePointData *points = (AliHLTSpacePointData*)memory->Allocate(pointsize);
59 // AliHLTDigitRowData *digits = (AliHLTDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
60 // AliHLTClustFinderNew *cf = new AliHLTClustFinderNew();
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 // AliHLTMemHandler *out= new AliHLTMemHandler();
73 // out->SetBinaryOutput(fname);
74 // out->Memory2Binary(npoints,points); //store the spacepoints
75 // out->CloseBinaryOutput();
84 ClassImp(AliHLTClustFinderNew)
86 AliHLTClustFinderNew::AliHLTClustFinderNew()
102 AliHLTClustFinderNew::~AliHLTClustFinderNew()
108 void AliHLTClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
112 fMaxNClusters = nmaxpoints;
113 fCurrentSlice = slice;
114 fCurrentPatch = patch;
115 fFirstRow = firstrow;
119 void AliHLTClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
123 fMaxNClusters = nmaxpoints;
124 fCurrentSlice = slice;
125 fCurrentPatch = patch;
126 fFirstRow=AliHLTTransform::GetFirstRow(patch);
127 fLastRow=AliHLTTransform::GetLastRow(patch);
130 void AliHLTClustFinderNew::SetOutputArray(AliHLTSpacePointData *pt)
132 //set pointer to output
133 fSpacePointData = pt;
136 void AliHLTClustFinderNew::Read(UInt_t ndigits,AliHLTDigitRowData *ptr)
139 fNDigitRowData = ndigits;
143 void AliHLTClustFinderNew::ProcessDigits()
145 //Loop over rows, and call processrow
146 AliHLTDigitRowData *tempPt = (AliHLTDigitRowData*)fDigitRowData;
148 for(Int_t i=fFirstRow; i<=fLastRow; i++)
151 if((Int_t)tempPt->fRow!=fCurrentRow){
152 LOG(AliHLTLog::kWarning,"AliHLTClustFinderNew::ProcessDigits","Digits")
153 <<"Row number should match! "<<tempPt->fRow<<" "<<fCurrentRow<<ENDLOG;
157 Byte_t *tmp = (Byte_t*)tempPt;
158 Int_t size = sizeof(AliHLTDigitRowData) + tempPt->fNDigit*sizeof(AliHLTDigitData);
160 tempPt = (AliHLTDigitRowData*)tmp;
162 LOG(AliHLTLog::kInformational,"AliHLTClustFinderNew::ProcessDigits","Space points")
163 <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
164 <<" patch "<<fCurrentPatch<<ENDLOG;
167 void AliHLTClustFinderNew::ProcessRow(AliHLTDigitRowData *tempPt)
170 UInt_t lastpad = 123456789;
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
176 AliClusterData **currentPt; //List of pointers to the current pad
177 AliClusterData **previousPt; //List of pointers to the previous pad
180 UInt_t nprevious=0,ncurrent=0,ntotal=0;
182 //Loop over sequences of this row:
183 for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
185 AliHLTDigitData *data = tempPt->fDigitData;
186 if(data[bin].fPad != lastpad)
191 if(currentPt == pad2)
201 nprevious = ncurrent;
203 if(bin[data].fPad != lastpad+1)
205 //this happens if there is a pad with no signal.
206 nprevious = ncurrent = 0;
208 lastpad = data[bin].fPad;
211 Bool_t newcluster = kTRUE;
212 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
213 UInt_t lastcharge=0,lastwas_falling=0;
218 redo: //This is a goto.
229 while(1) //Loop over current sequence
231 if(data[bin].fTime >= AliHLTTransform::GetNTimeBins())
233 LOG(AliHLTLog::kFatal,"AliHLTClustFinderNew::ProcessRow","Digits")
234 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
238 //Get the current ADC-value
239 UInt_t charge = data[bin].fCharge;
243 //Check if the last pixel in the sequence is smaller than this
244 if(charge > lastcharge)
252 else lastwas_falling = 1; //last pixel was larger than this
256 //Sum the total charge of this sequence
258 seqaverage += data[bin].fTime*charge;
259 seqerror += data[bin].fTime*data[bin].fTime*charge;
261 //Check where to stop:
262 if(bin >= tempPt->fNDigit - 1) //out of range
264 if(data[bin+1].fPad != data[bin].fPad) //new pad
266 if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
270 }//end loop over sequence
272 //Calculate mean of sequence:
275 seqmean = seqaverage/seqcharge;
278 LOG(AliHLTLog::kFatal,"AliHLTClustFinderNew::ProcessRow","Data")
279 <<"Error in data given to the cluster finder"<<ENDLOG;
284 //Calculate mean in pad direction:
285 Int_t padmean = seqcharge*data[bin].fPad;
286 Int_t paderror = data[bin].fPad*padmean;
288 //Compare with results on previous pad:
289 for(UInt_t p=0; p<nprevious; p++)
291 //dont merge sequences on the same pad twice
292 if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
294 Int_t difference = seqmean - previousPt[p]->fMean;
295 if(difference < -fMatch) break;
297 if(difference <= fMatch) //There is a match here!!
299 AliClusterData *local = previousPt[p];
303 if(seqcharge > local->fLastCharge)
305 if(local->fChargeFalling) //The previous pad was falling
307 break; //create a new cluster
311 local->fChargeFalling = 1;
312 local->fLastCharge = seqcharge;
315 //Don't create a new cluster, because we found a match
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;
328 currentPt[ncurrent] = local;
332 } //Checking for match at previous pad
333 } //Loop over results on previous pad.
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.
341 //Add to the clusterlist:
342 AliClusterData *tmp = &clusterlist[ntotal];
343 tmp->fTotalCharge = seqcharge;
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;
354 tmp->fChargeFalling = 0;
355 tmp->fLastCharge = seqcharge;
358 //Update list of pointers to previous pad:
359 currentPt[ncurrent] = &clusterlist[ntotal];
365 if(newbin >= 0) goto redo;
366 }//Loop over digits on this padrow
368 WriteClusters(ntotal,clusterlist);
371 void AliHLTClustFinderNew::WriteClusters(Int_t nclusters,AliClusterData *list)
373 //write cluster to output pointer
374 Int_t thisrow,thissector;
375 UInt_t counter = fNClusters;
377 for(int j=0; j<nclusters; j++)
379 if(!list[j].fFlags) continue; //discard single pad clusters
380 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
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
388 if(fCalcerr) { //calc the errors, otherwice take the fixed error
389 Int_t patch = AliHLTTransform::GetPatch(fCurrentRow);
390 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
392 LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Total charge")
393 <<"Zero total charge "<<q2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
397 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
400 LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Cluster width")
401 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
405 fpad2 = (sy2 + 1./12)*AliHLTTransform::GetPadPitchWidth(patch)*AliHLTTransform::GetPadPitchWidth(patch);
407 fpad2*=0.108; //constants are from offline studies
411 } else fpad2=sy2; //take the width not the error
413 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
416 LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Cluster width")
417 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
421 ftime2 = (sz2 + 1./12)*AliHLTTransform::GetZWidth()*AliHLTTransform::GetZWidth();
423 ftime2 *= 0.169; //constants are from offline studies
427 } else ftime2=sz2; //take the width, not the error
431 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
434 AliHLTTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
435 AliHLTTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
437 if(xyz[0]==0) LOG(AliHLTLog::kError,"AliHLTClustFinder","Cluster Finder")
438 <<AliHLTLog::kDec<<"Zero cluster"<<ENDLOG;
439 if(fNClusters >= fMaxNClusters)
441 LOG(AliHLTLog::kError,"AliHLTClustFinder::WriteClusters","Cluster Finder")
442 <<AliHLTLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
446 fSpacePointData[counter].fX = xyz[0];
447 fSpacePointData[counter].fY = xyz[1];
448 fSpacePointData[counter].fZ = xyz[2];
451 fSpacePointData[counter].fX = fCurrentRow;
452 fSpacePointData[counter].fY = fpad;
453 fSpacePointData[counter].fZ = ftime;
456 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
457 fSpacePointData[counter].fPadRow = fCurrentRow;
458 fSpacePointData[counter].fSigmaY2 = fpad2;
459 fSpacePointData[counter].fSigmaZ2 = ftime2;
461 Int_t patch=fCurrentPatch;
462 if(patch==-1) patch=0; //never store negative patch number
463 fSpacePointData[counter].fID = counter
464 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
467 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
469 fSpacePointData[counter].fTrackID[0] = trackID[0];
470 fSpacePointData[counter].fTrackID[1] = trackID[1];
471 fSpacePointData[counter].fTrackID[2] = trackID[2];
473 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
482 void AliHLTClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
485 AliHLTDigitRowData *rowPt = (AliHLTDigitRowData*)fDigitRowData;
487 trackID[0]=trackID[1]=trackID[2]=-2;
488 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
489 for(Int_t i=fFirstRow; i<=fLastRow; i++){
490 if(rowPt->fRow < (UInt_t)fCurrentRow){
491 AliHLTMemHandler::UpdateRowPointer(rowPt);
494 AliHLTDigitData *digPt = (AliHLTDigitData*)rowPt->fDigitData;
495 for(UInt_t j=0; j<rowPt->fNDigit; j++){
496 Int_t cpad = digPt[j].fPad;
497 Int_t ctime = digPt[j].fTime;
498 if(cpad != pad) continue;
499 if(ctime != time) continue;
501 trackID[0] = digPt[j].fTrackID[0];
502 trackID[1] = digPt[j].fTrackID[1];
503 trackID[2] = digPt[j].fTrackID[2];
505 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;