]> 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
Bug fix related to reconstruction of simulated data in raw format (E. Fragiacomo)
[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
aff6e981 556
557 if(fUnsorted){
558 fCurrentRow=list[j].fRow;
559 }
560
738c049f 561
a38a7850 562 if(fCalcerr) { //calc the errors, otherwice take the fixed error
563 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
564 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
565 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
566 sy2/=q2;
567 if(sy2 < 0) {
568 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
569 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
570 continue;
571 } else {
572 if(!fRawSP){
573 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
574 if(sy2 != 0){
575 fpad2*=0.108; //constants are from offline studies
576 if(patch<2)
577 fpad2*=2.07;
578 }
579 } else fpad2=sy2; //take the width not the error
580 }
581 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
582 sz2/=q2;
583 if(sz2 < 0){
584 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
585 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
586 continue;
587 } else {
588 if(!fRawSP){
589 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
590 if(sz2 != 0) {
591 ftime2 *= 0.169; //constants are from offline studies
592 if(patch<2)
593 ftime2 *= 1.77;
594 }
595 } else ftime2=sz2; //take the width, not the error
596 }
597 }
598 if(fStdout==kTRUE)
599 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
600
601 if(!fRawSP){
602 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
603 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
604
605 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
606 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
607 if(fNClusters >= fMaxNClusters)
608 {
609 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
610 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
611 return;
612 }
613
614 fSpacePointData[counter].fX = xyz[0];
615 fSpacePointData[counter].fY = xyz[1];
616 fSpacePointData[counter].fZ = xyz[2];
617
618 } else {
619 fSpacePointData[counter].fX = fCurrentRow;
620 fSpacePointData[counter].fY = fpad;
621 fSpacePointData[counter].fZ = ftime;
622 }
623
624 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
625 fSpacePointData[counter].fPadRow = fCurrentRow;
626 fSpacePointData[counter].fSigmaY2 = fpad2;
627 fSpacePointData[counter].fSigmaZ2 = ftime2;
628
44be0fde 629 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
738c049f 630 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
44be0fde 631
a38a7850 632 Int_t patch=fCurrentPatch;
633 if(patch==-1) patch=0; //never store negative patch number
634 fSpacePointData[counter].fID = counter
635 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
636
637#ifdef do_mc
638 Int_t trackID[3];
639 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
640
641 fSpacePointData[counter].fTrackID[0] = trackID[0];
642 fSpacePointData[counter].fTrackID[1] = trackID[1];
643 fSpacePointData[counter].fTrackID[2] = trackID[2];
644
645 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
646#endif
647
648 fNClusters++;
649 counter++;
650 }
651}
652
653// STILL TO FIX ----------------------------------------------------------------------------
654
655#ifdef do_mc
656void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
657{
658 //get mc id
659 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
660
661 trackID[0]=trackID[1]=trackID[2]=-2;
662 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
663 for(Int_t i=fFirstRow; i<=fLastRow; i++){
664 if(rowPt->fRow < (UInt_t)fCurrentRow){
665 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
666 continue;
667 }
668 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
669 for(UInt_t j=0; j<rowPt->fNDigit; j++){
670 Int_t cpad = digPt[j].fPad;
671 Int_t ctime = digPt[j].fTime;
672 if(cpad != pad) continue;
673 if(ctime != time) continue;
674
675 trackID[0] = digPt[j].fTrackID[0];
676 trackID[1] = digPt[j].fTrackID[1];
677 trackID[2] = digPt[j].fTrackID[2];
678
679 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
680 break;
681 }
682 break;
683 }
684}
685#endif
01f43166 686
687//----------------------------------Methods for the new unsorted way of reading the data --------------------------------
688
689void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
b1c46961 690 // see header file for function documentation
01f43166 691 fPadArray=padArray;
692}
693
694void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
695 //set input pointer
696 fPtr = (UChar_t*)ptr;
697 fSize = size;
698
699 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
700
701 fPadArray->SetDigitReader(fDigitReader);
b1c46961 702
703 if(fSignalThreshold>0){
704 fPadArray->SetSignalThreshold(fSignalThreshold);
705 }
706 if(fNSigmaThreshold>0){
707 fPadArray->SetNSigmaThreshold(fNSigmaThreshold);
708 }
01f43166 709 fPadArray->ReadData();
710}
711void AliHLTTPCClusterFinder::FindClusters(){
b1c46961 712 // see header file for function documentation
01f43166 713 fPadArray->FindClusterCandidates();
714 fPadArray->FindClusters(fMatch);
d8b2f74a 715
716 AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
a1dbf058 717 for(unsigned int i=0;i<fPadArray->fClusters.size();i++){
d8b2f74a 718 clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
719 clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
720 clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
721 clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
722 clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
723 clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
724 clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
725 clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
726 clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
727 clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
728 clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
01f43166 729 }
730 WriteClusters(fPadArray->fClusters.size(),clusterlist);
e83e889b 731 delete [] clusterlist;
6be21538 732 fPadArray->DataToDefault();
01f43166 733}
b1c46961 734Int_t AliHLTTPCClusterFinder::GetActivePads(AliHLTTPCPadArray::AliHLTTPCActivePads* activePads,Int_t maxActivePads){
735 // see header file for function documentation
736 return fPadArray->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)activePads , maxActivePads);
737}
01f43166 738void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
739{
740 //write cluster to output pointer
741 Int_t thisrow,thissector;
742 UInt_t counter = fNClusters;
743
744 for(int j=0; j<nclusters; j++)
745 {
746 if(!list[j].fFlags) continue; //discard single pad clusters
747 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
748
749 Float_t xyz[3];
750 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
751 Float_t fpad2=fXYErr*fXYErr; //fixed given error
752 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
753 Float_t ftime2=fZErr*fZErr; //fixed given error
754
d8b2f74a 755
01f43166 756 if(fCalcerr) { //calc the errors, otherwice take the fixed error
757 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
758 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
759 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
760 sy2/=q2;
761 if(sy2 < 0) {
762 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
763 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
764 continue;
765 } else {
766 if(!fRawSP){
767 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
768 if(sy2 != 0){
769 fpad2*=0.108; //constants are from offline studies
770 if(patch<2)
771 fpad2*=2.07;
772 }
773 } else fpad2=sy2; //take the width not the error
774 }
775 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
776 sz2/=q2;
777 if(sz2 < 0){
778 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
779 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
780 continue;
781 } else {
782 if(!fRawSP){
783 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
784 if(sz2 != 0) {
785 ftime2 *= 0.169; //constants are from offline studies
786 if(patch<2)
787 ftime2 *= 1.77;
788 }
789 } else ftime2=sz2; //take the width, not the error
790 }
791 }
792 if(fStdout==kTRUE)
793 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
794
795 if(!fRawSP){
796 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
797 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
798
799 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
800 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
801 if(fNClusters >= fMaxNClusters)
802 {
803 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
804 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
805 return;
806 }
807
808 fSpacePointData[counter].fX = xyz[0];
809 fSpacePointData[counter].fY = xyz[1];
810 fSpacePointData[counter].fZ = xyz[2];
811
812 } else {
813 fSpacePointData[counter].fX = fCurrentRow;
814 fSpacePointData[counter].fY = fpad;
815 fSpacePointData[counter].fZ = ftime;
816 }
817
818 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
819 fSpacePointData[counter].fPadRow = fCurrentRow;
820 fSpacePointData[counter].fSigmaY2 = fpad2;
821 fSpacePointData[counter].fSigmaZ2 = ftime2;
822
823 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
824 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
825
826 Int_t patch=fCurrentPatch;
827 if(patch==-1) patch=0; //never store negative patch number
828 fSpacePointData[counter].fID = counter
829 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
830
831#ifdef do_mc
832 Int_t trackID[3];
833 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
834
835 fSpacePointData[counter].fTrackID[0] = trackID[0];
836 fSpacePointData[counter].fTrackID[1] = trackID[1];
837 fSpacePointData[counter].fTrackID[2] = trackID[2];
838
839 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
840#endif
841
842 fNClusters++;
843 counter++;
844 }
845}