L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / src / AliHLTClustFinderNew.cxx
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 "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"
14
15 #if __GNUC__ >= 3
16 using namespace std;
17 #endif
18
19 /** \class AliHLTClustFinderNew
20 <pre>
21 //_____________________________________________________________
22 // AliHLTClustFinderNew
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 AliHLTDigitRowData 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 // 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);
55 //     UInt_t ndigits=0;
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();
76 //     delete out;
77 //     file->free();
78 //     delete cf;
79 //   }
80 // }
81 </pre> 
82 */
83
84 ClassImp(AliHLTClustFinderNew)
85
86 AliHLTClustFinderNew::AliHLTClustFinderNew()
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 AliHLTClustFinderNew::~AliHLTClustFinderNew()
103 {
104   //destructor
105   ;
106 }
107
108 void AliHLTClustFinderNew::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 AliHLTClustFinderNew::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=AliHLTTransform::GetFirstRow(patch);
127   fLastRow=AliHLTTransform::GetLastRow(patch);
128 }
129
130 void AliHLTClustFinderNew::SetOutputArray(AliHLTSpacePointData *pt)
131 {
132   //set pointer to output
133   fSpacePointData = pt;
134 }
135
136 void AliHLTClustFinderNew::Read(UInt_t ndigits,AliHLTDigitRowData *ptr)
137 {
138   //set input pointer
139   fNDigitRowData = ndigits;
140   fDigitRowData = ptr;
141 }
142
143 void AliHLTClustFinderNew::ProcessDigits()
144 {
145   //Loop over rows, and call processrow
146   AliHLTDigitRowData *tempPt = (AliHLTDigitRowData*)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(AliHLTLog::kWarning,"AliHLTClustFinderNew::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(AliHLTDigitRowData) + tempPt->fNDigit*sizeof(AliHLTDigitData);
159       tmp += size;
160       tempPt = (AliHLTDigitRowData*)tmp;
161     }
162   LOG(AliHLTLog::kInformational,"AliHLTClustFinderNew::ProcessDigits","Space points")
163     <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
164     <<" patch "<<fCurrentPatch<<ENDLOG;
165 }
166
167 void AliHLTClustFinderNew::ProcessRow(AliHLTDigitRowData *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       AliHLTDigitData *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 >= AliHLTTransform::GetNTimeBins())
232             {
233               LOG(AliHLTLog::kFatal,"AliHLTClustFinderNew::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].fTime*data[bin].fTime*charge;
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(AliHLTLog::kFatal,"AliHLTClustFinderNew::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 AliHLTClustFinderNew::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 = AliHLTTransform::GetPatch(fCurrentRow);
390         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
391         if (!q2) {
392             LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Total charge")
393               <<"Zero total charge "<<q2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
394             continue;
395
396         }
397         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
398         sy2/=q2;
399         if(sy2 < 0) {
400             LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Cluster width")
401               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
402             continue;
403         } else {
404           if(!fRawSP){
405             fpad2 = (sy2 + 1./12)*AliHLTTransform::GetPadPitchWidth(patch)*AliHLTTransform::GetPadPitchWidth(patch);
406             if(sy2 != 0){
407               fpad2*=0.108; //constants are from offline studies
408               if(patch<2)
409                 fpad2*=2.07;
410             }
411           } else fpad2=sy2; //take the width not the error
412         }
413         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
414         sz2/=q2;
415         if(sz2 < 0){
416           LOG(AliHLTLog::kError,"AliHLTClustFinderNew::WriteClusters","Cluster width")
417             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
418           continue;
419         } else {
420           if(!fRawSP){
421             ftime2 = (sz2 + 1./12)*AliHLTTransform::GetZWidth()*AliHLTTransform::GetZWidth();
422             if(sz2 != 0) {
423               ftime2 *= 0.169; //constants are from offline studies
424               if(patch<2)
425                 ftime2 *= 1.77;
426             }
427           } else ftime2=sz2; //take the width, not the error
428         }
429       }
430       if(fStdout==kTRUE)
431         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
432       
433       if(!fRawSP){
434         AliHLTTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
435         AliHLTTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
436         
437         if(xyz[0]==0) LOG(AliHLTLog::kError,"AliHLTClustFinder","Cluster Finder")
438           <<AliHLTLog::kDec<<"Zero cluster"<<ENDLOG;
439         if(fNClusters >= fMaxNClusters)
440           {
441             LOG(AliHLTLog::kError,"AliHLTClustFinder::WriteClusters","Cluster Finder")
442               <<AliHLTLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
443             return;
444           }  
445         
446         fSpacePointData[counter].fX = xyz[0];
447         fSpacePointData[counter].fY = xyz[1];
448         fSpacePointData[counter].fZ = xyz[2];
449         
450       } else {
451         fSpacePointData[counter].fX = fCurrentRow;
452         fSpacePointData[counter].fY = fpad;
453         fSpacePointData[counter].fZ = ftime;
454       }
455       
456       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
457       fSpacePointData[counter].fPadRow = fCurrentRow;
458       fSpacePointData[counter].fSigmaY2 = fpad2;
459       fSpacePointData[counter].fSigmaZ2  = ftime2;
460
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
465 #ifdef do_mc
466       Int_t trackID[3];
467       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
468
469       fSpacePointData[counter].fTrackID[0] = trackID[0];
470       fSpacePointData[counter].fTrackID[1] = trackID[1];
471       fSpacePointData[counter].fTrackID[2] = trackID[2];
472
473       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
474 #endif
475       
476       fNClusters++;
477       counter++;
478     }
479 }
480
481 #ifdef do_mc
482 void AliHLTClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
483 {
484   //get mc id
485   AliHLTDigitRowData *rowPt = (AliHLTDigitRowData*)fDigitRowData;
486   
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);
492       continue;
493     }
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;
500
501       trackID[0] = digPt[j].fTrackID[0];
502       trackID[1] = digPt[j].fTrackID[1];
503       trackID[2] = digPt[j].fTrackID[2];
504       
505       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
506       break;
507     }
508     break;
509   }
510 }
511 #endif