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