Adding setting default CDB storage and run number,
[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),
64defa03 139 fNSigmaThreshold(0),
2a083ac4 140 fNClusters(0),
141 fMaxNClusters(0),
142 fXYErr(0.2),
143 fZErr(0.3),
01f43166 144 fOccupancyLimit(1.0),
b1c46961 145 fPadArray(NULL),
64defa03 146 fUnsorted(0),
147 fActivePads()
a38a7850 148{
149 //constructor
46b33a24 150}
a38a7850 151
a38a7850 152AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
153{
154 //destructor
155}
156
157void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
158{
159 //init slice
160 fNClusters = 0;
161 fMaxNClusters = nmaxpoints;
162 fCurrentSlice = slice;
163 fCurrentPatch = patch;
164 fFirstRow = firstrow;
165 fLastRow = lastrow;
166}
167
168void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
169{
170 //init slice
171 fNClusters = 0;
172 fMaxNClusters = nmaxpoints;
173 fCurrentSlice = slice;
174 fCurrentPatch = patch;
175 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
176 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
177}
178
179void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
180{
181 //set pointer to output
182 fSpacePointData = pt;
183}
184
185void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
186 //set input pointer
187 fPtr = (UChar_t*)ptr;
188 fSize = size;
189}
190
191void AliHLTTPCClusterFinder::ProcessDigits()
192{
6b0fdc0e 193 int iResult=0;
a38a7850 194 bool readValue = true;
195 Int_t newRow = 0;
196 Int_t rowOffset = 0;
46b33a24 197 UShort_t time=0,newTime=0;
198 UInt_t pad=0,newPad=0;
199 AliHLTTPCSignal_t charge=0;
a38a7850 200
201 fNClusters = 0;
64defa03 202 fActivePads.clear();
a38a7850 203
204 // initialize block for reading packed data
67fada6b 205 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
206 if (iResult<0) return;
207
a38a7850 208 readValue = fDigitReader->Next();
209
5235c3e9 210 // Matthias 08.11.2006 the following return would cause termination without writing the
6c1a9d9e 211 // ClusterData and thus would block the component. I just want to have the commented line
5235c3e9 212 // here for information
213 //if (!readValue)return;
db16520a 214
a38a7850 215 pad = fDigitReader->GetPad();
216 time = fDigitReader->GetTime();
217 fCurrentRow = fDigitReader->GetRow();
218
219 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
220 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
221
222 fCurrentRow += rowOffset;
223
224 UInt_t lastpad = 123456789;
2a083ac4 225 const UInt_t kPadArraySize=5000;
226 const UInt_t kClusterListSize=10000;
84645eb0 227 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
228 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
229 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
a38a7850 230
231 AliClusterData **currentPt; //List of pointers to the current pad
232 AliClusterData **previousPt; //List of pointers to the previous pad
233 currentPt = pad2;
234 previousPt = pad1;
235 UInt_t nprevious=0,ncurrent=0,ntotal=0;
236
46b33a24 237 /* quick implementation of baseline calculation and zero suppression
238 open a pad object for each pad and delete it after processing.
239 later a list of pad objects with base line history can be used
240 The whole thing only works if we really get unprocessed raw data, if
241 the data is already zero suppressed, there might be gaps in the time
242 bins.
243 */
244 Int_t gatingGridOffset=50;
245 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
246 // just to make later conversion to a list of objects easier
84645eb0 247 AliHLTTPCPad* pCurrentPad=NULL;
248 if (fSignalThreshold>=0) {
249 pCurrentPad=&baseline;
250 baseline.SetThreshold(fSignalThreshold);
251 }
46b33a24 252
6b0fdc0e 253 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
254 iResult=0;
a38a7850 255
256 if(pad != lastpad){
257 //This is a new pad
258
259 //Switch the lists:
260 if(currentPt == pad2){
261 currentPt = pad1;
262 previousPt = pad2;
263 }
264 else {
265 currentPt = pad2;
266 previousPt = pad1;
267 }
268 nprevious = ncurrent;
269 ncurrent = 0;
270 if(pad != lastpad+1){
271 //this happens if there is a pad with no signal.
272 nprevious = ncurrent = 0;
273 }
274 lastpad = pad;
275 }
276
277 Bool_t newcluster = kTRUE;
278 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
2a083ac4 279 AliHLTTPCSignal_t lastcharge=0;
280 UInt_t bLastWasFalling=0;
a38a7850 281 Int_t newbin=-1;
282
283
284 if(fDeconvTime){
285 redo: //This is a goto.
286
287 if(newbin > -1){
288 //bin = newbin;
289 newbin = -1;
290 }
291
292 lastcharge=0;
2a083ac4 293 bLastWasFalling = 0;
a38a7850 294 }
295
6b0fdc0e 296 while(iResult>=0){ //Loop over time bins of current pad
46b33a24 297 // read all the values for one pad at once to calculate the base line
298 if (pCurrentPad) {
299 if (!pCurrentPad->IsStarted()) {
300 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
301 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
302 if ((pCurrentPad->StartEvent())>=0) {
303 do {
304 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
305 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
306 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
307 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
308 } while ((readValue = fDigitReader->Next())!=0);
309 }
310 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
311 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
84645eb0 312 //HLTDebug("no data available after zero suppression");
46b33a24 313 pCurrentPad->StopEvent();
314 pCurrentPad->ResetHistory();
315 break;
316 }
317 time=pCurrentPad->GetCurrentPosition();
318 if (time>pCurrentPad->GetSize()) {
319 HLTError("invalid time bin for pad");
320 break;
321 }
322 }
323 }
a38a7850 324
64defa03 325 if (fActivePads.size()==0 ||
326 fActivePads.back().fRow!=fCurrentRow-rowOffset ||
327 fActivePads.back().fPad!=pad) {
328 AliHLTTPCPadArray::AliHLTTPCActivePads entry;
329 entry.fRow=fCurrentRow-rowOffset;
330 entry.fPad=pad;
331 fActivePads.push_back(entry);
332 }
333
46b33a24 334 if (pCurrentPad) {
5235c3e9 335 Float_t occupancy=pCurrentPad->GetOccupancy();
336 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
337 if ( occupancy < fOccupancyLimit ) {
338 charge = pCurrentPad->GetCorrectedData();
339 } else {
340 charge = 0;
341 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
342 }
46b33a24 343 } else {
5235c3e9 344 charge = fDigitReader->GetSignal();
46b33a24 345 }
346 //HLTDebug("get next charge value: position %d charge %d", time, charge);
a38a7850 347
738c049f 348
349 // CHARGE DEBUG
350 if (fDigitReader->GetRow() == 90){
351///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
352 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
353
354 }
355
a38a7850 356 if(time >= AliHLTTPCTransform::GetNTimeBins()){
6b0fdc0e 357 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
358 iResult=-ERANGE;
a38a7850 359 }
46b33a24 360
361
a38a7850 362 //Get the current ADC-value
363 if(fDeconvTime){
364
365 //Check if the last pixel in the sequence is smaller than this
366 if(charge > lastcharge){
2a083ac4 367 if(bLastWasFalling){
a38a7850 368 newbin = 1;
369 break;
370 }
371 }
2a083ac4 372 else bLastWasFalling = 1; //last pixel was larger than this
a38a7850 373 lastcharge = charge;
374 }
375
376 //Sum the total charge of this sequence
377 seqcharge += charge;
378 seqaverage += time*charge;
379 seqerror += time*time*charge;
380
46b33a24 381 if (pCurrentPad) {
382
383 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
384 pCurrentPad->StopEvent();
385 pCurrentPad->ResetHistory();
386 if(readValue) {
387 newPad = fDigitReader->GetPad();
388 newTime = fDigitReader->GetTime();
389 newRow = fDigitReader->GetRow() + rowOffset;
390 }
391 break;
392 }
393
394 newPad=pCurrentPad->GetPadNumber();
395 newTime=pCurrentPad->GetCurrentPosition();
396 newRow=pCurrentPad->GetRowNumber();
397 } else {
a38a7850 398 readValue = fDigitReader->Next();
a38a7850 399 //Check where to stop:
400 if(!readValue) break; //No more value
401
402 newPad = fDigitReader->GetPad();
403 newTime = fDigitReader->GetTime();
404 newRow = fDigitReader->GetRow() + rowOffset;
46b33a24 405 }
a38a7850 406
407 if(newPad != pad)break; //new pad
408 if(newTime != time+1) break; //end of sequence
6b0fdc0e 409 if(iResult<0) break;
a38a7850 410
411 // pad = newpad; is equal
412 time = newTime;
413
414 }//end loop over sequence
415
46b33a24 416 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
417 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
418 if (seqcharge<=0) {
419 // with active zero suppression zero values are possible
420 continue;
421 }
a38a7850 422
423 //Calculate mean of sequence:
424 Int_t seqmean=0;
425 if(seqcharge)
426 seqmean = seqaverage/seqcharge;
427 else{
428 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
429 <<"Error in data given to the cluster finder"<<ENDLOG;
430 seqmean = 1;
431 seqcharge = 1;
432 }
433
434 //Calculate mean in pad direction:
435 Int_t padmean = seqcharge*pad;
436 Int_t paderror = pad*padmean;
437
a38a7850 438 //Compare with results on previous pad:
84645eb0 439 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
a38a7850 440
441 //dont merge sequences on the same pad twice
442 if(previousPt[p]->fLastMergedPad==pad) continue;
443
444 Int_t difference = seqmean - previousPt[p]->fMean;
445 if(difference < -fMatch) break;
446
447 if(difference <= fMatch){ //There is a match here!!
448 AliClusterData *local = previousPt[p];
449
450 if(fDeconvPad){
451 if(seqcharge > local->fLastCharge){
452 if(local->fChargeFalling){ //The previous pad was falling
453 break; //create a new cluster
454 }
455 }
456 else local->fChargeFalling = 1;
457 local->fLastCharge = seqcharge;
458 }
459
460 //Don't create a new cluster, because we found a match
461 newcluster = kFALSE;
462
463 //Update cluster on current pad with the matching one:
464 local->fTotalCharge += seqcharge;
465 local->fPad += padmean;
466 local->fPad2 += paderror;
467 local->fTime += seqaverage;
468 local->fTime2 += seqerror;
469 local->fMean = seqmean;
470 local->fFlags++; //means we have more than one pad
471 local->fLastMergedPad = pad;
472
473 currentPt[ncurrent] = local;
474 ncurrent++;
475
476 break;
477 } //Checking for match at previous pad
478 } //Loop over results on previous pad.
479
84645eb0 480 if(newcluster && ncurrent<kPadArraySize){
a38a7850 481 //Start a new cluster. Add it to the clusterlist, and update
482 //the list of pointers to clusters in current pad.
483 //current pad will be previous pad on next pad.
484
485 //Add to the clusterlist:
486 AliClusterData *tmp = &clusterlist[ntotal];
487 tmp->fTotalCharge = seqcharge;
488 tmp->fPad = padmean;
489 tmp->fPad2 = paderror;
490 tmp->fTime = seqaverage;
491 tmp->fTime2 = seqerror;
492 tmp->fMean = seqmean;
493 tmp->fFlags = 0; //flags for single pad clusters
494 tmp->fLastMergedPad = pad;
495
496 if(fDeconvPad){
497 tmp->fChargeFalling = 0;
498 tmp->fLastCharge = seqcharge;
499 }
500
501 //Update list of pointers to previous pad:
502 currentPt[ncurrent] = &clusterlist[ntotal];
503 ntotal++;
504 ncurrent++;
505 }
506
507 if(fDeconvTime)
508 if(newbin >= 0) goto redo;
509
510 // to prevent endless loop
511 if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24 512 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850 513 break;
514 }
515
516
517 if(!readValue) break; //No more value
84645eb0 518
519 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
520 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
521 break;
522 }
a38a7850 523
524 if(fCurrentRow != newRow){
525 WriteClusters(ntotal,clusterlist);
526
527 lastpad = 123456789;
528
529 currentPt = pad2;
530 previousPt = pad1;
531 nprevious=0;
532 ncurrent=0;
533 ntotal=0;
534
535 fCurrentRow = newRow;
536 }
537
538 pad = newPad;
539 time = newTime;
db16520a 540
a38a7850 541 } // END while(readValue)
542
543 WriteClusters(ntotal,clusterlist);
544
46b33a24 545 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850 546
547} // ENDEND
548
549void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
550{
551 //write cluster to output pointer
8f8bf0af 552 Int_t thisrow=-1,thissector=-1;
a38a7850 553 UInt_t counter = fNClusters;
554
555 for(int j=0; j<nclusters; j++)
556 {
d8b2f74a 557
558
559
a38a7850 560 if(!list[j].fFlags) continue; //discard single pad clusters
561 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
562
563 Float_t xyz[3];
564 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
565 Float_t fpad2=fXYErr*fXYErr; //fixed given error
566 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
567 Float_t ftime2=fZErr*fZErr; //fixed given error
568
738c049f 569
aff6e981 570
571 if(fUnsorted){
572 fCurrentRow=list[j].fRow;
573 }
574
738c049f 575
a38a7850 576 if(fCalcerr) { //calc the errors, otherwice take the fixed error
577 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
578 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
579 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
580 sy2/=q2;
581 if(sy2 < 0) {
582 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
583 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
584 continue;
585 } else {
586 if(!fRawSP){
587 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
588 if(sy2 != 0){
589 fpad2*=0.108; //constants are from offline studies
590 if(patch<2)
591 fpad2*=2.07;
592 }
593 } else fpad2=sy2; //take the width not the error
594 }
595 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
596 sz2/=q2;
597 if(sz2 < 0){
598 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
599 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
600 continue;
601 } else {
602 if(!fRawSP){
603 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
604 if(sz2 != 0) {
605 ftime2 *= 0.169; //constants are from offline studies
606 if(patch<2)
607 ftime2 *= 1.77;
608 }
609 } else ftime2=sz2; //take the width, not the error
610 }
611 }
612 if(fStdout==kTRUE)
613 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
614
615 if(!fRawSP){
616 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
617 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
618
619 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
620 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
621 if(fNClusters >= fMaxNClusters)
622 {
623 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
624 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
625 return;
626 }
627
628 fSpacePointData[counter].fX = xyz[0];
629 fSpacePointData[counter].fY = xyz[1];
630 fSpacePointData[counter].fZ = xyz[2];
631
632 } else {
633 fSpacePointData[counter].fX = fCurrentRow;
634 fSpacePointData[counter].fY = fpad;
635 fSpacePointData[counter].fZ = ftime;
636 }
637
638 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
639 fSpacePointData[counter].fPadRow = fCurrentRow;
640 fSpacePointData[counter].fSigmaY2 = fpad2;
641 fSpacePointData[counter].fSigmaZ2 = ftime2;
642
44be0fde 643 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
738c049f 644 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
44be0fde 645
a38a7850 646 Int_t patch=fCurrentPatch;
647 if(patch==-1) patch=0; //never store negative patch number
648 fSpacePointData[counter].fID = counter
649 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
650
651#ifdef do_mc
652 Int_t trackID[3];
653 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
654
655 fSpacePointData[counter].fTrackID[0] = trackID[0];
656 fSpacePointData[counter].fTrackID[1] = trackID[1];
657 fSpacePointData[counter].fTrackID[2] = trackID[2];
658
659 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
660#endif
661
662 fNClusters++;
663 counter++;
664 }
665}
666
667// STILL TO FIX ----------------------------------------------------------------------------
668
669#ifdef do_mc
670void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
671{
672 //get mc id
673 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
674
675 trackID[0]=trackID[1]=trackID[2]=-2;
676 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
677 for(Int_t i=fFirstRow; i<=fLastRow; i++){
678 if(rowPt->fRow < (UInt_t)fCurrentRow){
679 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
680 continue;
681 }
682 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
683 for(UInt_t j=0; j<rowPt->fNDigit; j++){
684 Int_t cpad = digPt[j].fPad;
685 Int_t ctime = digPt[j].fTime;
686 if(cpad != pad) continue;
687 if(ctime != time) continue;
688
689 trackID[0] = digPt[j].fTrackID[0];
690 trackID[1] = digPt[j].fTrackID[1];
691 trackID[2] = digPt[j].fTrackID[2];
692
693 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
694 break;
695 }
696 break;
697 }
698}
699#endif
01f43166 700
701//----------------------------------Methods for the new unsorted way of reading the data --------------------------------
702
64defa03 703void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray)
704{
b1c46961 705 // see header file for function documentation
01f43166 706 fPadArray=padArray;
707}
708
64defa03 709void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
710{
01f43166 711 //set input pointer
712 fPtr = (UChar_t*)ptr;
713 fSize = size;
714
715 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
716
717 fPadArray->SetDigitReader(fDigitReader);
b1c46961 718
719 if(fSignalThreshold>0){
720 fPadArray->SetSignalThreshold(fSignalThreshold);
721 }
722 if(fNSigmaThreshold>0){
723 fPadArray->SetNSigmaThreshold(fNSigmaThreshold);
724 }
01f43166 725 fPadArray->ReadData();
726}
64defa03 727
728void AliHLTTPCClusterFinder::FindClusters()
729{
b1c46961 730 // see header file for function documentation
01f43166 731 fPadArray->FindClusterCandidates();
732 fPadArray->FindClusters(fMatch);
d8b2f74a 733
734 AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
a1dbf058 735 for(unsigned int i=0;i<fPadArray->fClusters.size();i++){
d8b2f74a 736 clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
737 clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
738 clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
739 clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
740 clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
741 clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
742 clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
743 clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
744 clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
745 clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
746 clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
01f43166 747 }
748 WriteClusters(fPadArray->fClusters.size(),clusterlist);
e83e889b 749 delete [] clusterlist;
6be21538 750 fPadArray->DataToDefault();
01f43166 751}
64defa03 752
753Int_t AliHLTTPCClusterFinder::GetActivePads(AliHLTTPCPadArray::AliHLTTPCActivePads* activePads,Int_t maxActivePads)
754{
b1c46961 755 // see header file for function documentation
64defa03 756 Int_t iResult=0;
757 if (fPadArray) {
758 iResult=fPadArray->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)activePads , maxActivePads);
759 } else if ((iResult=fActivePads.size())>0) {
760 if (iResult>maxActivePads) {
761 HLTWarning("target array (%d) not big enough to receive %d active pad descriptors", maxActivePads, iResult);
762 iResult=maxActivePads;
763 }
764 memcpy(activePads, &fActivePads[0], iResult*sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads));
765 }
766
767 return iResult;
b1c46961 768}
64defa03 769
01f43166 770void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
771{
772 //write cluster to output pointer
773 Int_t thisrow,thissector;
774 UInt_t counter = fNClusters;
775
776 for(int j=0; j<nclusters; j++)
777 {
778 if(!list[j].fFlags) continue; //discard single pad clusters
779 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
780
781 Float_t xyz[3];
782 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
783 Float_t fpad2=fXYErr*fXYErr; //fixed given error
784 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
785 Float_t ftime2=fZErr*fZErr; //fixed given error
786
d8b2f74a 787
01f43166 788 if(fCalcerr) { //calc the errors, otherwice take the fixed error
789 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
790 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
791 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
792 sy2/=q2;
793 if(sy2 < 0) {
794 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
795 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
796 continue;
797 } else {
798 if(!fRawSP){
799 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
800 if(sy2 != 0){
801 fpad2*=0.108; //constants are from offline studies
802 if(patch<2)
803 fpad2*=2.07;
804 }
805 } else fpad2=sy2; //take the width not the error
806 }
807 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
808 sz2/=q2;
809 if(sz2 < 0){
810 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
811 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
812 continue;
813 } else {
814 if(!fRawSP){
815 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
816 if(sz2 != 0) {
817 ftime2 *= 0.169; //constants are from offline studies
818 if(patch<2)
819 ftime2 *= 1.77;
820 }
821 } else ftime2=sz2; //take the width, not the error
822 }
823 }
824 if(fStdout==kTRUE)
825 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
826
827 if(!fRawSP){
828 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
829 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
830
831 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
832 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
833 if(fNClusters >= fMaxNClusters)
834 {
835 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
836 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
837 return;
838 }
839
840 fSpacePointData[counter].fX = xyz[0];
841 fSpacePointData[counter].fY = xyz[1];
842 fSpacePointData[counter].fZ = xyz[2];
843
844 } else {
845 fSpacePointData[counter].fX = fCurrentRow;
846 fSpacePointData[counter].fY = fpad;
847 fSpacePointData[counter].fZ = ftime;
848 }
849
850 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
851 fSpacePointData[counter].fPadRow = fCurrentRow;
852 fSpacePointData[counter].fSigmaY2 = fpad2;
853 fSpacePointData[counter].fSigmaZ2 = ftime2;
854
855 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
856 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
857
858 Int_t patch=fCurrentPatch;
859 if(patch==-1) patch=0; //never store negative patch number
860 fSpacePointData[counter].fID = counter
861 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
862
863#ifdef do_mc
864 Int_t trackID[3];
865 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
866
867 fSpacePointData[counter].fTrackID[0] = trackID[0];
868 fSpacePointData[counter].fTrackID[1] = trackID[1];
869 fSpacePointData[counter].fTrackID[2] = trackID[2];
870
871 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
872#endif
873
874 fNClusters++;
875 counter++;
876 }
877}