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