coding conventions and compilation warnings
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClusterFinder.cxx
CommitLineData
a38a7850 1// @(#) $Id$
4aa41877 2// Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
a38a7850 3
600e6a1b 4/**************************************************************************
9be2600f 5 * This file is property of and copyright by the ALICE HLT Project *
6 * ALICE Experiment at CERN, All rights reserved. *
600e6a1b 7 * *
01f43166 8 * Primary Authors: Anders Vestbo, Constantin Loizides, Jochen Thaeder *
9 * Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
10 * for The ALICE HLT Project. *
600e6a1b 11 * *
12 * Permission to use, copy, modify and distribute this software and its *
13 * documentation strictly for non-commercial purposes is hereby granted *
14 * without fee, provided that the above copyright notice appears in all *
15 * copies and that both the copyright notice and this permission notice *
16 * appear in the supporting documentation. The authors make no claims *
17 * about the suitability of this software for any purpose. It is *
18 * provided "as is" without express or implied warranty. *
19 **************************************************************************/
a38a7850 20
600e6a1b 21/** @file AliHLTTPCClusterFinder.cxx
01f43166 22 @author Anders Vestbo, Constantin Loizides, Jochen Thaeder
23 Kenneth Aamodt kenneth.aamodt@student.uib.no
600e6a1b 24 @date
25 @brief Cluster Finder for the TPC
26*/
a38a7850 27
28#include "AliHLTTPCDigitReader.h"
a38a7850 29#include "AliHLTTPCRootTypes.h"
30#include "AliHLTTPCLogging.h"
31#include "AliHLTTPCClusterFinder.h"
32#include "AliHLTTPCDigitData.h"
33#include "AliHLTTPCTransform.h"
34#include "AliHLTTPCSpacePointData.h"
35#include "AliHLTTPCMemHandler.h"
46b33a24 36#include "AliHLTTPCPad.h"
01f43166 37#include "AliHLTTPCPadArray.h"
a38a7850 38
39#if __GNUC__ >= 3
40using namespace std;
41#endif
42
43/** \class AliHLTTPCClusterFinder
a38a7850 44//
45// The current cluster finder for HLT
46// (Based on STAR L3)
01f43166 47//
48//Basically we have two versions for the cluster finder now.
49//The default version, reads the data pad by pad, and find the
50//clusters as it reads the data. The other version has now been
51//developed to cope with unsorted data. New methods for the unsorted
52//version can be found at the end of the default one i the source file.
53//Currently the new version is only build to manage zero-suppressed data.
54//More functionality will be added later.
55//
a38a7850 56// The cluster finder is initialized with the Init function,
57// providing the slice and patch information to work on.
58//
59// The input is a provided by the AliHLTTPCDigitReader class,
60// using the init() funktion, and the next() funktion in order
61// to get the next bin. Either packed or unpacked data can be
62// processed, dependent if one uses AliHLTTPCDigitReaderPacked
63// class or AliHLTTPCDigitReaderUnpacked class in the
64// Clusterfinder Component.
65// The resulting space points will be in the
66// array given by the SetOutputArray function.
67//
68// There are several setters which control the behaviour:
69//
70// - SetXYError(Float_t): set fixed error in XY direction
71// - SetZError(Float_t): set fixed error in Z direction
72// (used if errors are not calculated)
73// - SetDeconv(Bool_t): switch on/off deconvolution
74// - SetThreshold(UInt_t): set charge threshold for cluster
75// - SetMatchWidth(UInt_t): set the match distance in
76// time for sequences to be merged
77// - SetSTDOutput(Bool_t): switch on/off output about found clusters
78// - SetCalcErr(Bool_t): switch on/off calculation of
79// space point errors (or widths in raw system)
80// - SetRawSP(Bool_t): switch on/off convertion to raw system
81//
82//
83// Example Usage:
84//
85// AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
86// file->SetAliInput(digitfile); //give some input file
87// for(int slice=0; slice<=35; slice++){
88// for(int patch=0; pat<6; pat++){
89// file->Init(slice,patch);
90// UInt_t ndigits=0;
91// UInt_t maxclusters=100000;
92// UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
93// AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
94// AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
95// AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
96// cf->SetMatchWidth(2);
97// cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
98// cf->SetSTDOutput(kTRUE); //Some output to standard IO
99// cf->SetRawSP(kFALSE); //Convert space points to local system
100// cf->SetThreshold(5); //Threshold of cluster charge
101// cf->SetDeconv(kTRUE); //Deconv in pad and time direction
102// cf->SetCalcErr(kTRUE); //Calculate the errors of the spacepoints
103// cf->SetOutputArray(points); //Move the spacepoints to the array
104// cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
105// cf->ProcessDigits(); //process the rows given by init
106// Int_t npoints = cf->GetNumberOfClusters();
107// AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
108// out->SetBinaryOutput(fname);
109// out->Memory2Binary(npoints,points); //store the spacepoints
110// out->CloseBinaryOutput();
111// delete out;
112// file->free();
113// delete cf;
114// }
115// }
a38a7850 116*/
117
118ClassImp(AliHLTTPCClusterFinder)
119
120AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
46b33a24 121 :
2a083ac4 122 fSpacePointData(NULL),
123 fDigitReader(NULL),
124 fPtr(NULL),
125 fSize(0),
46b33a24 126 fDeconvTime(kTRUE),
2a083ac4 127 fDeconvPad(kTRUE),
46b33a24 128 fStdout(kFALSE),
129 fCalcerr(kTRUE),
130 fRawSP(kFALSE),
131 fFirstRow(0),
132 fLastRow(0),
2a083ac4 133 fCurrentRow(0),
134 fCurrentSlice(0),
135 fCurrentPatch(0),
136 fMatch(1),
137 fThreshold(10),
138 fSignalThreshold(-1),
139 fNClusters(0),
140 fMaxNClusters(0),
141 fXYErr(0.2),
142 fZErr(0.3),
01f43166 143 fOccupancyLimit(1.0),
144 fPadArray(NULL)
a38a7850 145{
146 //constructor
46b33a24 147}
a38a7850 148
46b33a24 149AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
150 :
2a083ac4 151 fSpacePointData(NULL),
152 fDigitReader(src.fDigitReader),
153 fPtr(NULL),
154 fSize(src.fSize),
46b33a24 155 fDeconvTime(src.fDeconvTime),
2a083ac4 156 fDeconvPad(src.fDeconvPad),
46b33a24 157 fStdout(src.fStdout),
158 fCalcerr(src.fCalcerr),
159 fRawSP(src.fRawSP),
160 fFirstRow(src.fFirstRow),
161 fLastRow(src.fLastRow),
2a083ac4 162 fCurrentRow(src.fCurrentRow),
163 fCurrentSlice(src.fCurrentSlice),
164 fCurrentPatch(src.fCurrentPatch),
165 fMatch(src.fMatch),
166 fThreshold(src.fThreshold),
167 fSignalThreshold(src.fSignalThreshold),
168 fNClusters(src.fNClusters),
169 fMaxNClusters(src.fMaxNClusters),
170 fXYErr(src.fXYErr),
171 fZErr(src.fZErr),
01f43166 172 fOccupancyLimit(src.fOccupancyLimit),
173 fPadArray(NULL)
46b33a24 174{
175}
176
177AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
178{
179 fMatch=src.fMatch;
180 fThreshold=src.fThreshold;
84645eb0 181 fSignalThreshold=src.fSignalThreshold;
5235c3e9 182 fOccupancyLimit=src.fOccupancyLimit;
46b33a24 183 fXYErr=src.fXYErr;
184 fZErr=src.fZErr;
185 fDeconvPad=src.fDeconvPad;
186 fDeconvTime=src.fDeconvTime;
187 fStdout=src.fStdout;
188 fCalcerr=src.fCalcerr;
189 fRawSP=src.fRawSP;
190 fFirstRow=src.fFirstRow;
191 fLastRow=src.fLastRow;
192 fDigitReader=src.fDigitReader;
193 return (*this);
a38a7850 194}
195
196AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
197{
198 //destructor
199}
200
201void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
202{
203 //init slice
204 fNClusters = 0;
205 fMaxNClusters = nmaxpoints;
206 fCurrentSlice = slice;
207 fCurrentPatch = patch;
208 fFirstRow = firstrow;
209 fLastRow = lastrow;
210}
211
212void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
213{
214 //init slice
215 fNClusters = 0;
216 fMaxNClusters = nmaxpoints;
217 fCurrentSlice = slice;
218 fCurrentPatch = patch;
219 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
220 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
221}
222
223void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
224{
225 //set pointer to output
226 fSpacePointData = pt;
227}
228
229void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
230 //set input pointer
231 fPtr = (UChar_t*)ptr;
232 fSize = size;
233}
234
235void AliHLTTPCClusterFinder::ProcessDigits()
236{
6b0fdc0e 237 int iResult=0;
a38a7850 238 bool readValue = true;
239 Int_t newRow = 0;
240 Int_t rowOffset = 0;
46b33a24 241 UShort_t time=0,newTime=0;
242 UInt_t pad=0,newPad=0;
243 AliHLTTPCSignal_t charge=0;
a38a7850 244
245 fNClusters = 0;
246
247 // initialize block for reading packed data
db16520a 248 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
a38a7850 249 readValue = fDigitReader->Next();
250
5235c3e9 251 // Matthias 08.11.2006 the following return would cause termination without writing the
6c1a9d9e 252 // ClusterData and thus would block the component. I just want to have the commented line
5235c3e9 253 // here for information
254 //if (!readValue)return;
db16520a 255
a38a7850 256 pad = fDigitReader->GetPad();
257 time = fDigitReader->GetTime();
258 fCurrentRow = fDigitReader->GetRow();
259
260 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
261 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
262
263 fCurrentRow += rowOffset;
264
265 UInt_t lastpad = 123456789;
2a083ac4 266 const UInt_t kPadArraySize=5000;
267 const UInt_t kClusterListSize=10000;
84645eb0 268 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
269 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
270 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
a38a7850 271
272 AliClusterData **currentPt; //List of pointers to the current pad
273 AliClusterData **previousPt; //List of pointers to the previous pad
274 currentPt = pad2;
275 previousPt = pad1;
276 UInt_t nprevious=0,ncurrent=0,ntotal=0;
277
46b33a24 278 /* quick implementation of baseline calculation and zero suppression
279 open a pad object for each pad and delete it after processing.
280 later a list of pad objects with base line history can be used
281 The whole thing only works if we really get unprocessed raw data, if
282 the data is already zero suppressed, there might be gaps in the time
283 bins.
284 */
285 Int_t gatingGridOffset=50;
286 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
287 // just to make later conversion to a list of objects easier
84645eb0 288 AliHLTTPCPad* pCurrentPad=NULL;
289 if (fSignalThreshold>=0) {
290 pCurrentPad=&baseline;
291 baseline.SetThreshold(fSignalThreshold);
292 }
46b33a24 293
6b0fdc0e 294 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
295 iResult=0;
a38a7850 296
297 if(pad != lastpad){
298 //This is a new pad
299
300 //Switch the lists:
301 if(currentPt == pad2){
302 currentPt = pad1;
303 previousPt = pad2;
304 }
305 else {
306 currentPt = pad2;
307 previousPt = pad1;
308 }
309 nprevious = ncurrent;
310 ncurrent = 0;
311 if(pad != lastpad+1){
312 //this happens if there is a pad with no signal.
313 nprevious = ncurrent = 0;
314 }
315 lastpad = pad;
316 }
317
318 Bool_t newcluster = kTRUE;
319 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
2a083ac4 320 AliHLTTPCSignal_t lastcharge=0;
321 UInt_t bLastWasFalling=0;
a38a7850 322 Int_t newbin=-1;
323
324
325 if(fDeconvTime){
326 redo: //This is a goto.
327
328 if(newbin > -1){
329 //bin = newbin;
330 newbin = -1;
331 }
332
333 lastcharge=0;
2a083ac4 334 bLastWasFalling = 0;
a38a7850 335 }
336
6b0fdc0e 337 while(iResult>=0){ //Loop over time bins of current pad
46b33a24 338 // read all the values for one pad at once to calculate the base line
339 if (pCurrentPad) {
340 if (!pCurrentPad->IsStarted()) {
341 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
342 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
343 if ((pCurrentPad->StartEvent())>=0) {
344 do {
345 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
346 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
347 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
348 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
349 } while ((readValue = fDigitReader->Next())!=0);
350 }
351 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
352 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
84645eb0 353 //HLTDebug("no data available after zero suppression");
46b33a24 354 pCurrentPad->StopEvent();
355 pCurrentPad->ResetHistory();
356 break;
357 }
358 time=pCurrentPad->GetCurrentPosition();
359 if (time>pCurrentPad->GetSize()) {
360 HLTError("invalid time bin for pad");
361 break;
362 }
363 }
364 }
a38a7850 365
46b33a24 366 if (pCurrentPad) {
5235c3e9 367 Float_t occupancy=pCurrentPad->GetOccupancy();
368 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
369 if ( occupancy < fOccupancyLimit ) {
370 charge = pCurrentPad->GetCorrectedData();
371 } else {
372 charge = 0;
373 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
374 }
46b33a24 375 } else {
5235c3e9 376 charge = fDigitReader->GetSignal();
46b33a24 377 }
378 //HLTDebug("get next charge value: position %d charge %d", time, charge);
a38a7850 379
738c049f 380
381 // CHARGE DEBUG
382 if (fDigitReader->GetRow() == 90){
383///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
384 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
385
386 }
387
a38a7850 388 if(time >= AliHLTTPCTransform::GetNTimeBins()){
6b0fdc0e 389 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
390 iResult=-ERANGE;
a38a7850 391 }
46b33a24 392
393
a38a7850 394 //Get the current ADC-value
395 if(fDeconvTime){
396
397 //Check if the last pixel in the sequence is smaller than this
398 if(charge > lastcharge){
2a083ac4 399 if(bLastWasFalling){
a38a7850 400 newbin = 1;
401 break;
402 }
403 }
2a083ac4 404 else bLastWasFalling = 1; //last pixel was larger than this
a38a7850 405 lastcharge = charge;
406 }
407
408 //Sum the total charge of this sequence
409 seqcharge += charge;
410 seqaverage += time*charge;
411 seqerror += time*time*charge;
412
46b33a24 413 if (pCurrentPad) {
414
415 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
416 pCurrentPad->StopEvent();
417 pCurrentPad->ResetHistory();
418 if(readValue) {
419 newPad = fDigitReader->GetPad();
420 newTime = fDigitReader->GetTime();
421 newRow = fDigitReader->GetRow() + rowOffset;
422 }
423 break;
424 }
425
426 newPad=pCurrentPad->GetPadNumber();
427 newTime=pCurrentPad->GetCurrentPosition();
428 newRow=pCurrentPad->GetRowNumber();
429 } else {
a38a7850 430 readValue = fDigitReader->Next();
a38a7850 431 //Check where to stop:
432 if(!readValue) break; //No more value
433
434 newPad = fDigitReader->GetPad();
435 newTime = fDigitReader->GetTime();
436 newRow = fDigitReader->GetRow() + rowOffset;
46b33a24 437 }
a38a7850 438
439 if(newPad != pad)break; //new pad
440 if(newTime != time+1) break; //end of sequence
6b0fdc0e 441 if(iResult<0) break;
a38a7850 442
443 // pad = newpad; is equal
444 time = newTime;
445
446 }//end loop over sequence
447
46b33a24 448 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
449 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
450 if (seqcharge<=0) {
451 // with active zero suppression zero values are possible
452 continue;
453 }
a38a7850 454
455 //Calculate mean of sequence:
456 Int_t seqmean=0;
457 if(seqcharge)
458 seqmean = seqaverage/seqcharge;
459 else{
460 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
461 <<"Error in data given to the cluster finder"<<ENDLOG;
462 seqmean = 1;
463 seqcharge = 1;
464 }
465
466 //Calculate mean in pad direction:
467 Int_t padmean = seqcharge*pad;
468 Int_t paderror = pad*padmean;
469
a38a7850 470 //Compare with results on previous pad:
84645eb0 471 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
a38a7850 472
473 //dont merge sequences on the same pad twice
474 if(previousPt[p]->fLastMergedPad==pad) continue;
475
476 Int_t difference = seqmean - previousPt[p]->fMean;
477 if(difference < -fMatch) break;
478
479 if(difference <= fMatch){ //There is a match here!!
480 AliClusterData *local = previousPt[p];
481
482 if(fDeconvPad){
483 if(seqcharge > local->fLastCharge){
484 if(local->fChargeFalling){ //The previous pad was falling
485 break; //create a new cluster
486 }
487 }
488 else local->fChargeFalling = 1;
489 local->fLastCharge = seqcharge;
490 }
491
492 //Don't create a new cluster, because we found a match
493 newcluster = kFALSE;
494
495 //Update cluster on current pad with the matching one:
496 local->fTotalCharge += seqcharge;
497 local->fPad += padmean;
498 local->fPad2 += paderror;
499 local->fTime += seqaverage;
500 local->fTime2 += seqerror;
501 local->fMean = seqmean;
502 local->fFlags++; //means we have more than one pad
503 local->fLastMergedPad = pad;
504
505 currentPt[ncurrent] = local;
506 ncurrent++;
507
508 break;
509 } //Checking for match at previous pad
510 } //Loop over results on previous pad.
511
84645eb0 512 if(newcluster && ncurrent<kPadArraySize){
a38a7850 513 //Start a new cluster. Add it to the clusterlist, and update
514 //the list of pointers to clusters in current pad.
515 //current pad will be previous pad on next pad.
516
517 //Add to the clusterlist:
518 AliClusterData *tmp = &clusterlist[ntotal];
519 tmp->fTotalCharge = seqcharge;
520 tmp->fPad = padmean;
521 tmp->fPad2 = paderror;
522 tmp->fTime = seqaverage;
523 tmp->fTime2 = seqerror;
524 tmp->fMean = seqmean;
525 tmp->fFlags = 0; //flags for single pad clusters
526 tmp->fLastMergedPad = pad;
527
528 if(fDeconvPad){
529 tmp->fChargeFalling = 0;
530 tmp->fLastCharge = seqcharge;
531 }
532
533 //Update list of pointers to previous pad:
534 currentPt[ncurrent] = &clusterlist[ntotal];
535 ntotal++;
536 ncurrent++;
537 }
538
539 if(fDeconvTime)
540 if(newbin >= 0) goto redo;
541
542 // to prevent endless loop
543 if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24 544 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850 545 break;
546 }
547
548
549 if(!readValue) break; //No more value
84645eb0 550
551 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
552 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
553 break;
554 }
a38a7850 555
556 if(fCurrentRow != newRow){
557 WriteClusters(ntotal,clusterlist);
558
559 lastpad = 123456789;
560
561 currentPt = pad2;
562 previousPt = pad1;
563 nprevious=0;
564 ncurrent=0;
565 ntotal=0;
566
567 fCurrentRow = newRow;
568 }
569
570 pad = newPad;
571 time = newTime;
db16520a 572
a38a7850 573 } // END while(readValue)
574
575 WriteClusters(ntotal,clusterlist);
576
46b33a24 577 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850 578
579} // ENDEND
580
581void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
582{
583 //write cluster to output pointer
584 Int_t thisrow,thissector;
585 UInt_t counter = fNClusters;
586
587 for(int j=0; j<nclusters; j++)
588 {
d8b2f74a 589
590
591
a38a7850 592 if(!list[j].fFlags) continue; //discard single pad clusters
593 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
594
595 Float_t xyz[3];
596 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
597 Float_t fpad2=fXYErr*fXYErr; //fixed given error
598 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
599 Float_t ftime2=fZErr*fZErr; //fixed given error
600
738c049f 601
d8b2f74a 602#if UNSORTED
603 fCurrentRow=list[j].fRow;
604#endif
738c049f 605
a38a7850 606 if(fCalcerr) { //calc the errors, otherwice take the fixed error
607 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
608 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
609 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
610 sy2/=q2;
611 if(sy2 < 0) {
612 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
613 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
614 continue;
615 } else {
616 if(!fRawSP){
617 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
618 if(sy2 != 0){
619 fpad2*=0.108; //constants are from offline studies
620 if(patch<2)
621 fpad2*=2.07;
622 }
623 } else fpad2=sy2; //take the width not the error
624 }
625 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
626 sz2/=q2;
627 if(sz2 < 0){
628 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
629 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
630 continue;
631 } else {
632 if(!fRawSP){
633 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
634 if(sz2 != 0) {
635 ftime2 *= 0.169; //constants are from offline studies
636 if(patch<2)
637 ftime2 *= 1.77;
638 }
639 } else ftime2=sz2; //take the width, not the error
640 }
641 }
642 if(fStdout==kTRUE)
643 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
644
645 if(!fRawSP){
646 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
647 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
648
649 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
650 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
651 if(fNClusters >= fMaxNClusters)
652 {
653 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
654 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
655 return;
656 }
657
658 fSpacePointData[counter].fX = xyz[0];
659 fSpacePointData[counter].fY = xyz[1];
660 fSpacePointData[counter].fZ = xyz[2];
661
662 } else {
663 fSpacePointData[counter].fX = fCurrentRow;
664 fSpacePointData[counter].fY = fpad;
665 fSpacePointData[counter].fZ = ftime;
666 }
667
668 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
669 fSpacePointData[counter].fPadRow = fCurrentRow;
670 fSpacePointData[counter].fSigmaY2 = fpad2;
671 fSpacePointData[counter].fSigmaZ2 = ftime2;
672
44be0fde 673 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
738c049f 674 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
44be0fde 675
a38a7850 676 Int_t patch=fCurrentPatch;
677 if(patch==-1) patch=0; //never store negative patch number
678 fSpacePointData[counter].fID = counter
679 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
680
681#ifdef do_mc
682 Int_t trackID[3];
683 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
684
685 fSpacePointData[counter].fTrackID[0] = trackID[0];
686 fSpacePointData[counter].fTrackID[1] = trackID[1];
687 fSpacePointData[counter].fTrackID[2] = trackID[2];
688
689 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
690#endif
691
692 fNClusters++;
693 counter++;
694 }
695}
696
697// STILL TO FIX ----------------------------------------------------------------------------
698
699#ifdef do_mc
700void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
701{
702 //get mc id
703 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
704
705 trackID[0]=trackID[1]=trackID[2]=-2;
706 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
707 for(Int_t i=fFirstRow; i<=fLastRow; i++){
708 if(rowPt->fRow < (UInt_t)fCurrentRow){
709 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
710 continue;
711 }
712 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
713 for(UInt_t j=0; j<rowPt->fNDigit; j++){
714 Int_t cpad = digPt[j].fPad;
715 Int_t ctime = digPt[j].fTime;
716 if(cpad != pad) continue;
717 if(ctime != time) continue;
718
719 trackID[0] = digPt[j].fTrackID[0];
720 trackID[1] = digPt[j].fTrackID[1];
721 trackID[2] = digPt[j].fTrackID[2];
722
723 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
724 break;
725 }
726 break;
727 }
728}
729#endif
01f43166 730
731//----------------------------------Methods for the new unsorted way of reading the data --------------------------------
732
733void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
734 fPadArray=padArray;
735}
736
737void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
738 //set input pointer
739 fPtr = (UChar_t*)ptr;
740 fSize = size;
741
742 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
743
744 fPadArray->SetDigitReader(fDigitReader);
745 fPadArray->ReadData();
746}
747void AliHLTTPCClusterFinder::FindClusters(){
748 fPadArray->FindClusterCandidates();
749 fPadArray->FindClusters(fMatch);
d8b2f74a 750
751 AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
a1dbf058 752 for(unsigned int i=0;i<fPadArray->fClusters.size();i++){
d8b2f74a 753 clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
754 clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
755 clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
756 clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
757 clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
758 clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
759 clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
760 clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
761 clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
762 clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
763 clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
01f43166 764 }
765 WriteClusters(fPadArray->fClusters.size(),clusterlist);
e83e889b 766 delete [] clusterlist;
6be21538 767 fPadArray->DataToDefault();
01f43166 768}
769void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
770{
771 //write cluster to output pointer
772 Int_t thisrow,thissector;
773 UInt_t counter = fNClusters;
774
775 for(int j=0; j<nclusters; j++)
776 {
777 if(!list[j].fFlags) continue; //discard single pad clusters
778 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
779
780 Float_t xyz[3];
781 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
782 Float_t fpad2=fXYErr*fXYErr; //fixed given error
783 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
784 Float_t ftime2=fZErr*fZErr; //fixed given error
785
d8b2f74a 786
01f43166 787 if(fCalcerr) { //calc the errors, otherwice take the fixed error
788 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
789 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
790 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
791 sy2/=q2;
792 if(sy2 < 0) {
793 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
794 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
795 continue;
796 } else {
797 if(!fRawSP){
798 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
799 if(sy2 != 0){
800 fpad2*=0.108; //constants are from offline studies
801 if(patch<2)
802 fpad2*=2.07;
803 }
804 } else fpad2=sy2; //take the width not the error
805 }
806 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
807 sz2/=q2;
808 if(sz2 < 0){
809 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
810 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
811 continue;
812 } else {
813 if(!fRawSP){
814 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
815 if(sz2 != 0) {
816 ftime2 *= 0.169; //constants are from offline studies
817 if(patch<2)
818 ftime2 *= 1.77;
819 }
820 } else ftime2=sz2; //take the width, not the error
821 }
822 }
823 if(fStdout==kTRUE)
824 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
825
826 if(!fRawSP){
827 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
828 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
829
830 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
831 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
832 if(fNClusters >= fMaxNClusters)
833 {
834 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
835 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
836 return;
837 }
838
839 fSpacePointData[counter].fX = xyz[0];
840 fSpacePointData[counter].fY = xyz[1];
841 fSpacePointData[counter].fZ = xyz[2];
842
843 } else {
844 fSpacePointData[counter].fX = fCurrentRow;
845 fSpacePointData[counter].fY = fpad;
846 fSpacePointData[counter].fZ = ftime;
847 }
848
849 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
850 fSpacePointData[counter].fPadRow = fCurrentRow;
851 fSpacePointData[counter].fSigmaY2 = fpad2;
852 fSpacePointData[counter].fSigmaZ2 = ftime2;
853
854 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
855 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
856
857 Int_t patch=fCurrentPatch;
858 if(patch==-1) patch=0; //never store negative patch number
859 fSpacePointData[counter].fID = counter
860 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
861
862#ifdef do_mc
863 Int_t trackID[3];
864 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
865
866 fSpacePointData[counter].fTrackID[0] = trackID[0];
867 fSpacePointData[counter].fTrackID[1] = trackID[1];
868 fSpacePointData[counter].fTrackID[2] = trackID[2];
869
870 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
871#endif
872
873 fNClusters++;
874 counter++;
875 }
876}