]> 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
Added functionality to read the magnetic field value, and setting the recoparam objec...
[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
297174de 4//**************************************************************************
5//* This file is property of and copyright by the ALICE HLT Project *
6//* ALICE Experiment at CERN, All rights reserved. *
7//* *
8//* Primary Authors: Anders Vestbo, Constantin Loizides *
9//* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
10//* Kalliopi Kanaki *
11//* for The ALICE HLT Project. *
12//* *
13//* Permission to use, copy, modify and distribute this software and its *
14//* documentation strictly for non-commercial purposes is hereby granted *
15//* without fee, provided that the above copyright notice appears in all *
16//* copies and that both the copyright notice and this permission notice *
17//* appear in the supporting documentation. The authors make no claims *
18//* about the suitability of this software for any purpose. It is *
19//* provided "as is" without express or implied warranty. *
20//**************************************************************************
a38a7850 21
600e6a1b 22/** @file AliHLTTPCClusterFinder.cxx
98034d0d 23 @author Kenneth Aamodt, Kalliopi Kanaki
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"
a38a7850 32#include "AliHLTTPCSpacePointData.h"
33#include "AliHLTTPCMemHandler.h"
46b33a24 34#include "AliHLTTPCPad.h"
8252a538 35#include <sys/time.h>
deba5d85 36#include <algorithm>
6a5788b9 37#include <cmath>
5ee8e602 38#include "AliTPCcalibDB.h"
39#include "AliTPCTransform.h"
a38a7850 40
41#if __GNUC__ >= 3
42using namespace std;
43#endif
44
a38a7850 45ClassImp(AliHLTTPCClusterFinder)
46
47AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
46b33a24 48 :
2fdb1ae7 49 fClustersHWAddressVector(),
50 fRowPadVector(),
2a083ac4 51 fSpacePointData(NULL),
52 fDigitReader(NULL),
53 fPtr(NULL),
54 fSize(0),
6b15c309 55 fDeconvTime(kFALSE),
56 fDeconvPad(kFALSE),
46b33a24 57 fStdout(kFALSE),
58 fCalcerr(kTRUE),
59 fRawSP(kFALSE),
60 fFirstRow(0),
61 fLastRow(0),
2a083ac4 62 fCurrentRow(0),
63 fCurrentSlice(0),
64 fCurrentPatch(0),
65 fMatch(1),
66 fThreshold(10),
2a083ac4 67 fNClusters(0),
68 fMaxNClusters(0),
69 fXYErr(0.2),
70 fZErr(0.3),
01f43166 71 fOccupancyLimit(1.0),
64defa03 72 fUnsorted(0),
a74855c2 73 fVectorInitialized(kFALSE),
8252a538 74 fClusters(),
deba5d85 75 fClustersMCInfo(),
76 fMCDigits(),
a74855c2 77 fNumberOfPadsInRow(NULL),
78 fNumberOfRows(0),
2fdb1ae7 79 fRowOfFirstCandidate(0),
6b15c309 80 fDoPadSelection(kFALSE),
81 fFirstTimeBin(0),
82 fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
83 fTotalChargeOfPreviousClusterCandidate(0),
deba5d85 84 fChargeOfCandidatesFalling(kFALSE),
85 f32BitFormat(kFALSE),
86 fDoMC(kFALSE),
5ee8e602 87 fClusterMCVector(),
6a5788b9 88 fOfflineTransform(NULL),
f14729cb 89 fOfflineTPCRecoParam(),
6a5788b9 90 fTimeMeanDiff(2)
a38a7850 91{
8252a538 92 //constructor
5ee8e602 93 fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
94 if(!fOfflineTransform){
95 HLTError("AliHLTTPCClusterFinder(): Offline transform not in AliTPCcalibDB.");
96 }
f14729cb 97 else{
98 fOfflineTPCRecoParam.SetUseExBCorrection(1);
99 fOfflineTPCRecoParam.SetUseTOFCorrection(1);
100 fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
101 }
46b33a24 102}
a38a7850 103
a912b63b 104AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder(){
105 // see header file for class documentation
106
a38a7850 107 //destructor
8252a538 108 if(fVectorInitialized){
109 DeInitializePadArray();
110 }
111 if(fNumberOfPadsInRow){
112 delete [] fNumberOfPadsInRow;
113 fNumberOfPadsInRow=NULL;
114 }
a38a7850 115}
a912b63b 116
117void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints){
118 // see header file for class documentation
119
a38a7850 120 //init slice
121 fNClusters = 0;
122 fMaxNClusters = nmaxpoints;
123 fCurrentSlice = slice;
124 fCurrentPatch = patch;
a912b63b 125 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
126 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
37c9ee0f 127
128 fClusters.clear();
129 fClustersMCInfo.clear();
130 fMCDigits.clear();
131 fClusterMCVector.clear();
a38a7850 132}
133
8252a538 134void AliHLTTPCClusterFinder::InitializePadArray(){
135 // see header file for class documentation
136
137 if(fCurrentPatch>5||fCurrentPatch<0){
138 HLTFatal("Patch is not set");
139 return;
140 }
141
142 HLTDebug("Patch number=%d",fCurrentPatch);
143
144 fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
145 fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
146
147 fNumberOfRows=fLastRow-fFirstRow+1;
148 fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
149
150 memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
151
a74855c2 152 for(UInt_t i=0;i<fNumberOfRows;i++){
8252a538 153 fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
154 AliHLTTPCPadVector tmpRow;
288f1d48 155 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
8252a538 156 AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
157 tmpPad->SetID(i,j);
158 tmpRow.push_back(tmpPad);
159 }
160 fRowPadVector.push_back(tmpRow);
161 }
162 fVectorInitialized=kTRUE;
163}
164
a912b63b 165Int_t AliHLTTPCClusterFinder::DeInitializePadArray(){
8252a538 166 // see header file for class documentation
a912b63b 167
a74855c2 168 for(UInt_t i=0;i<fNumberOfRows;i++){
288f1d48 169 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
8252a538 170 delete fRowPadVector[i][j];
171 fRowPadVector[i][j]=NULL;
172 }
173 fRowPadVector[i].clear();
174 }
175 fRowPadVector.clear();
176 return 1;
177}
178
179
a38a7850 180
a912b63b 181void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt){
182 // see header file for class documentation
a38a7850 183 //set pointer to output
184 fSpacePointData = pt;
185}
186
a912b63b 187
188void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
189 // see header file for class documentation
a38a7850 190 //set input pointer
191 fPtr = (UChar_t*)ptr;
192 fSize = size;
a38a7850 193
a912b63b 194 if(!fVectorInitialized){
195 InitializePadArray();
6b15c309 196 }
a912b63b 197
198 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
199 HLTError("failed setting up digit reader (InitBlock)");
200 return;
84645eb0 201 }
a912b63b 202
203 while(fDigitReader->NextChannel()){
204 UInt_t row=fDigitReader->GetRow();
205 UInt_t pad=fDigitReader->GetPad();
a38a7850 206
a912b63b 207 if(row>=fRowPadVector.size()){
208 HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
209 continue;
a38a7850 210 }
a912b63b 211 if(pad>=fRowPadVector[row].size()){
212 HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
213 continue;
a38a7850 214 }
215
a912b63b 216 while(fDigitReader->NextBunch()){
217 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
218 UInt_t time = fDigitReader->GetTime();
219 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
deba5d85 220 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
221 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
222 // The same is true for the function ReadDataUnsortedDeconvoluteTime() below.
223 // In addition the signals are organized in the opposite direction
224 if(f32BitFormat){
225 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
226 AliHLTTPCClusters candidate;
227 for(Int_t i=fDigitReader->GetBunchSize()-1;i>=0;i--){
228 candidate.fTotalCharge+=bunchData[i];
229 candidate.fTime += time*bunchData[i];
230 candidate.fTime2 += time*time*bunchData[i];
231 if(bunchData[i]>candidate.fQMax){
232 candidate.fQMax=bunchData[i];
233 }
234 time++;
235 }
236 if(candidate.fTotalCharge>0){
237 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
238 candidate.fPad=candidate.fTotalCharge*pad;
239 candidate.fPad2=candidate.fPad*pad;
240 candidate.fLastMergedPad=pad;
241 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
242 }
243 if(fRowPadVector[row][pad] != NULL){
244 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
a912b63b 245 }
46b33a24 246 }
deba5d85 247 else{
248 const UInt_t *bunchData= fDigitReader->GetSignals();
249 AliHLTTPCClusters candidate;
025ba4dc 250 const AliHLTTPCDigitData* digits = NULL;
251 if(fDoMC && (digits = fDigitReader->GetBunchDigits())!=NULL){
deba5d85 252 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
253 candidate.fTotalCharge+=bunchData[i];
254 candidate.fTime += time*bunchData[i];
255 candidate.fTime2 += time*time*bunchData[i];
256 if(bunchData[i]>candidate.fQMax){
257 candidate.fQMax=bunchData[i];
258 }
259 fMCDigits.push_back(digits[i]);
260 time++;
261 }
262 }
263 else{
264 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
265 candidate.fTotalCharge+=bunchData[i];
266 candidate.fTime += time*bunchData[i];
267 candidate.fTime2 += time*time*bunchData[i];
268 if(bunchData[i]>candidate.fQMax){
269 candidate.fQMax=bunchData[i];
270 }
271 time++;
272 }
273 }
274 if(candidate.fTotalCharge>0){
275 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
276 candidate.fPad=candidate.fTotalCharge*pad;
277 candidate.fPad2=candidate.fPad*pad;
278 candidate.fLastMergedPad=pad;
279 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
280 }
281 if(fRowPadVector[row][pad] != NULL){
282 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
283 if(fDoMC){
284 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
285 fMCDigits.clear();
286 }
287 }
46b33a24 288 }
289 }
290 }
a912b63b 291 }
292 }
293}
a38a7850 294
a912b63b 295void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size){
296 // see header file for class documentation
738c049f 297
a912b63b 298 //set input pointer
299 fPtr = (UChar_t*)ptr;
300 fSize = size;
738c049f 301
a912b63b 302 if(!fVectorInitialized){
303 InitializePadArray();
304 }
738c049f 305
a912b63b 306 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
307 HLTError("failed setting up digit reader (InitBlock)");
308 return;
309 }
310
311 while(fDigitReader->NextChannel()){
312 UInt_t row=fDigitReader->GetRow();
313 UInt_t pad=fDigitReader->GetPad();
46b33a24 314
a912b63b 315 while(fDigitReader->NextBunch()){
316 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
317 UInt_t time = fDigitReader->GetTime();
318 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
319 Int_t indexInBunchData=0;
320 Bool_t moreDataInBunch=kFALSE;
321 UInt_t prevSignal=0;
322 Bool_t signalFalling=kFALSE;
deba5d85 323
324 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
325 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
326 // The same is true for the function ReadDataUnsorted() above.
327 // In addition the signals are organized in the opposite direction
328 if(f32BitFormat){
329 indexInBunchData = fDigitReader->GetBunchSize();
330 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
5ee8e602 331
deba5d85 332 do{
333 AliHLTTPCClusters candidate;
334 //for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
335 for(Int_t i=indexInBunchData;i>=0;i--){
336 // Checks if one need to deconvolute the signals
337 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
338 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
339 moreDataInBunch=kTRUE;
340 prevSignal=0;
341 }
342 break;
343 }
344
345 // Checks if the signal is 0, then quit processing the data.
346 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
a912b63b 347 moreDataInBunch=kTRUE;
348 prevSignal=0;
deba5d85 349 break;
350 }
351
352 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
353 signalFalling=kTRUE;
354 }
355 candidate.fTotalCharge+=bunchData[i];
356 candidate.fTime += time*bunchData[i];
357 candidate.fTime2 += time*time*bunchData[i];
358 if(bunchData[i]>candidate.fQMax){
359 candidate.fQMax=bunchData[i];
a912b63b 360 }
deba5d85 361
362 prevSignal=bunchData[i];
363 time++;
364 indexInBunchData--;
a912b63b 365 }
deba5d85 366 if(candidate.fTotalCharge>0){
367 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
368 candidate.fPad=candidate.fTotalCharge*pad;
369 candidate.fPad2=candidate.fPad*pad;
370 candidate.fLastMergedPad=pad;
371 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
a912b63b 372 }
deba5d85 373 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
374 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
375 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
376 moreDataInBunch=kFALSE;
a912b63b 377 }
deba5d85 378 }while(moreDataInBunch);
379 }
380 else{
381 const UInt_t *bunchData= fDigitReader->GetSignals();
382 do{
383 AliHLTTPCClusters candidate;
384 for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
385 // Checks if one need to deconvolute the signals
386 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
387 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
388 moreDataInBunch=kTRUE;
389 prevSignal=0;
390 }
391 break;
392 }
393
394 // Checks if the signal is 0, then quit processing the data.
395 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
396 moreDataInBunch=kTRUE;
397 prevSignal=0;
398 break;
399 }
400
401 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
402 signalFalling=kTRUE;
403 }
5ee8e602 404
deba5d85 405 candidate.fTotalCharge+=bunchData[i];
406 candidate.fTime += time*bunchData[i];
407 candidate.fTime2 += time*time*bunchData[i];
408 if(bunchData[i]>candidate.fQMax){
409 candidate.fQMax=bunchData[i];
410 }
411
412 prevSignal=bunchData[i];
413 time++;
414 indexInBunchData++;
a912b63b 415 }
deba5d85 416 if(candidate.fTotalCharge>0){
417 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
418 candidate.fPad=candidate.fTotalCharge*pad;
419 candidate.fPad2=candidate.fPad*pad;
420 candidate.fLastMergedPad=pad;
421 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
422 }
423 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
424 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
425 moreDataInBunch=kFALSE;
426 }
427 }while(moreDataInBunch);
428 }
a912b63b 429 }
430 }
431 }
432 }
433}
434
435Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
436 // see header file for class documentation
437
438 //Checking if we have a match on the next pad
439 for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
2d6b17ff 440 if(nextPad->fUsedClusterCandidates[candidateNumber] == 1){
441 continue;
442 }
a912b63b 443 AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
6a5788b9 444 // if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
445
446 if( abs((Int_t)(cluster->fMean - candidate->fMean)) <= fTimeMeanDiff ){
a912b63b 447 if(fDeconvPad){
448 if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
449 fChargeOfCandidatesFalling=kTRUE;
450 }
451 if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
452 return kFALSE;
453 }
454 }
455 cluster->fMean=candidate->fMean;
456 cluster->fTotalCharge+=candidate->fTotalCharge;
457 cluster->fTime += candidate->fTime;
458 cluster->fTime2 += candidate->fTime2;
459 cluster->fPad+=candidate->fPad;
460 cluster->fPad2+=candidate->fPad2;
461 cluster->fLastMergedPad=candidate->fPad;
462 if(candidate->fQMax>cluster->fQMax){
463 cluster->fQMax=candidate->fQMax;
464 }
deba5d85 465 if(fDoMC){
466 FillMCClusterVector(nextPad->GetCandidateDigits(candidateNumber));
467 }
468
a912b63b 469 if(fDoPadSelection){
470 UInt_t rowNo = nextPad->GetRowNumber();
471 UInt_t padNo = nextPad->GetPadNumber();
472 if(padNo-1>0){
473 fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
474 fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
475 }
476 fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
477 fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
478 fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
479 fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
480 }
481
482 //setting the matched pad to used
483 nextPad->fUsedClusterCandidates[candidateNumber]=1;
484 nextPadToRead++;
485 if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
486 nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
487 ComparePads(nextPad,cluster,nextPadToRead);
488 }
489 else{
490 return kFALSE;
491 }
492 }
a912b63b 493 }
494 return kFALSE;
495}
496
497Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
498 // see header file for class documentation
499
500 Int_t counter=0;
501 for(UInt_t row=0;row<fNumberOfRows;row++){
502 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
503 if(fRowPadVector[row][pad]->fSelectedPad){
504 if(counter<maxHWadd){
505 hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
506 counter++;
507 }
508 else{
509 HLTWarning("To many hardwareaddresses, skip adding");
510 }
511
512 }
513 }
514 }
515 return counter;
516}
deba5d85 517
518
519Int_t AliHLTTPCClusterFinder::FillOutputMCInfo(AliHLTTPCClusterFinder::ClusterMCInfo * outputMCInfo, Int_t maxNumberOfClusterMCInfo){
520 // see header file for class documentation
521
522 Int_t counter=0;
523
524 for(UInt_t mc=0;mc<fClustersMCInfo.size();mc++){
525 if(counter<maxNumberOfClusterMCInfo){
526 outputMCInfo[counter] = fClustersMCInfo[mc];
527 counter++;
528 }
529 else{
530 HLTWarning("To much MCInfo has been added (no more space), skip adding");
531 }
532 }
533 return counter;
534}
a912b63b 535
536void AliHLTTPCClusterFinder::FindClusters(){
537 // see header file for function documentation
538
539 AliHLTTPCClusters* tmpCandidate=NULL;
540 for(UInt_t row=0;row<fNumberOfRows;row++){
541 fRowOfFirstCandidate=row;
542 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
543 AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
544 for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
545 if(tmpPad->fUsedClusterCandidates[candidate]){
546 continue;
547 }
548 tmpCandidate=&tmpPad->fClusterCandidates[candidate];
549 UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
550
deba5d85 551 if(fDoMC){
552 fClusterMCVector.clear();
553 FillMCClusterVector(tmpPad->GetCandidateDigits(candidate));
554 }
555
a912b63b 556 ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
557 if(tmpCandidate->fTotalCharge>tmpTotalCharge){
558 //we have a cluster
559 fClusters.push_back(*tmpCandidate);
deba5d85 560 if(fDoMC){
561 //sort the vector (large->small) according to weight and remove elements above 2 (keep 0 1 and 2)
37c9ee0f 562 sort(fClusterMCVector.begin(),fClusterMCVector.end(), MCWeight::CompareWeights );
deba5d85 563 ClusterMCInfo tmpClusterMCInfo;
564
565 MCWeight zeroMC;
566 zeroMC.fMCID=-1;
567 zeroMC.fWeight=0;
568
569 if(fClusterMCVector.size()>0){
570 tmpClusterMCInfo.fClusterID[0]=fClusterMCVector.at(0);
571 }
572 else{
573 tmpClusterMCInfo.fClusterID[0]=zeroMC;
574 }
575
576 if(fClusterMCVector.size()>1){
577 tmpClusterMCInfo.fClusterID[1]=fClusterMCVector.at(1);
578 }
579 else{
580 tmpClusterMCInfo.fClusterID[1]=zeroMC;
581 }
582
583 if(fClusterMCVector.size()>2){
584 tmpClusterMCInfo.fClusterID[2]=fClusterMCVector.at(2);
585 }
586 else{
587 tmpClusterMCInfo.fClusterID[2]=zeroMC;
588 }
589
590 fClustersMCInfo.push_back(tmpClusterMCInfo);
591 }
592
a912b63b 593 }
594 }
595 tmpPad->ClearCandidates();
596 }
597 fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
598 }
599
600 HLTInfo("Found %d clusters.",fClusters.size());
601
602 //TODO: Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
603
604 AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
605 for(unsigned int i=0;i<fClusters.size();i++){
606 clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
607 clusterlist[i].fPad = fClusters[i].fPad;
608 clusterlist[i].fPad2 = fClusters[i].fPad2;
609 clusterlist[i].fTime = fClusters[i].fTime;
610 clusterlist[i].fTime2 = fClusters[i].fTime2;
611 clusterlist[i].fMean = fClusters[i].fMean;
612 clusterlist[i].fFlags = fClusters[i].fFlags;
613 clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
614 clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
615 clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
616 clusterlist[i].fRow = fClusters[i].fRowNumber;
617 clusterlist[i].fQMax = fClusters[i].fQMax;
618 }
619
620 WriteClusters(fClusters.size(),clusterlist);
621 delete [] clusterlist;
622 fClusters.clear();
623}
624
625
4f43db26 626Bool_t AliHLTTPCClusterFinder::UpdateCalibDB(){
627
628 //update the db
629 AliTPCcalibDB::Instance()->Update();
630
631 //uptate the transform class
632 AliTPCTransform * tmp = AliTPCcalibDB::Instance()->GetTransform();
633 if(!tmp){
634 HLTError("AliHLTTPCClusterFinder::UpdateCAlibDB: Offline transform not in AliTPCcalibDB.");
635 return 0;
636 }
637 fOfflineTransform = tmp;
638 return 1;
639}
640
a912b63b 641//---------------------------------- Under this line the old sorted clusterfinder functions can be found --------------------------------
642
643
644void AliHLTTPCClusterFinder::PrintClusters(){
645 // see header file for class documentation
646
647 for(size_t i=0;i<fClusters.size();i++){
648 HLTInfo("Cluster number: %d",i);
649 HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
650 HLTInfo("Total Charge: %d",fClusters[i].fTotalCharge);
651 HLTInfo("fPad: %d",fClusters[i].fPad);
652 HLTInfo("PadError: %d",fClusters[i].fPad2);
653 HLTInfo("TimeMean: %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
654 HLTInfo("TimeError: %d",fClusters[i].fTime2);
655 HLTInfo("EndOfCluster:");
656 }
657}
658
deba5d85 659void AliHLTTPCClusterFinder::FillMCClusterVector(vector<AliHLTTPCDigitData> digitData){
660
661 for(UInt_t d=0;d<digitData.size();d++){
662 Int_t nIDsInDigit = (digitData.at(d).fTrackID[0]>=0) + (digitData.at(d).fTrackID[1]>=0) + (digitData.at(d).fTrackID[2]>=0);
663 for(Int_t id=0; id<3; id++){
664 if(digitData.at(d).fTrackID[id]>=0){
665 Bool_t matchFound = kFALSE;
666 MCWeight mc;
667 mc.fMCID = digitData.at(d).fTrackID[id];
668 mc.fWeight = ((Float_t)digitData.at(d).fCharge)/nIDsInDigit;
669 for(UInt_t i=0;i<fClusterMCVector.size();i++){
670 if(mc.fMCID == fClusterMCVector.at(i).fMCID){
671 fClusterMCVector.at(i).fWeight += mc.fWeight;
672 matchFound = kTRUE;
673 }
674 }
675 if(matchFound == kFALSE){
676 fClusterMCVector.push_back(mc);
677 }
678 }
679 }
680 }
681}
682
683
a912b63b 684void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
685 //set input pointer
686 fPtr = (UChar_t*)ptr;
687 fSize = size;
688}
689
690void AliHLTTPCClusterFinder::ProcessDigits(){
691 // see header file for class documentation
692
693 int iResult=0;
694 bool readValue = true;
695 Int_t newRow = 0;
696 Int_t rowOffset = 0;
697 UShort_t time=0,newTime=0;
698 UInt_t pad=0,newPad=0;
699 AliHLTTPCSignal_t charge=0;
700
701 fNClusters = 0;
702
703 // initialize block for reading packed data
704 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
705 if (iResult<0) return;
706
707 readValue = fDigitReader->Next();
708
709 // Matthias 08.11.2006 the following return would cause termination without writing the
710 // ClusterData and thus would block the component. I just want to have the commented line
711 // here for information
712 //if (!readValue)return;
713
714 pad = fDigitReader->GetPad();
715 time = fDigitReader->GetTime();
716 fCurrentRow = fDigitReader->GetRow();
717
718 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
719 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
720
721 fCurrentRow += rowOffset;
722
723 UInt_t lastpad = 123456789;
724 const UInt_t kPadArraySize=5000;
725 const UInt_t kClusterListSize=10000;
726 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
727 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
728 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
729
730 AliClusterData **currentPt; //List of pointers to the current pad
731 AliClusterData **previousPt; //List of pointers to the previous pad
732 currentPt = pad2;
733 previousPt = pad1;
734 UInt_t nprevious=0,ncurrent=0,ntotal=0;
735
736 /* quick implementation of baseline calculation and zero suppression
737 open a pad object for each pad and delete it after processing.
738 later a list of pad objects with base line history can be used
739 The whole thing only works if we really get unprocessed raw data, if
740 the data is already zero suppressed, there might be gaps in the time
741 bins.
742 */
743 Int_t gatingGridOffset=50;
744 if(fFirstTimeBin>0){
745 gatingGridOffset=fFirstTimeBin;
746 }
747 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
748 // just to make later conversion to a list of objects easier
749 AliHLTTPCPad* pCurrentPad=NULL;
750 /*
751 if (fSignalThreshold>=0) {
752 pCurrentPad=&baseline;
753 baseline.SetThreshold(fSignalThreshold);
754 }
755 */
756 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
757 iResult=0;
758
759 if(pad != lastpad){
760 //This is a new pad
761
762 //Switch the lists:
763 if(currentPt == pad2){
764 currentPt = pad1;
765 previousPt = pad2;
766 }
767 else {
768 currentPt = pad2;
769 previousPt = pad1;
770 }
771 nprevious = ncurrent;
772 ncurrent = 0;
773 if(pad != lastpad+1){
774 //this happens if there is a pad with no signal.
775 nprevious = ncurrent = 0;
776 }
777 lastpad = pad;
778 }
779
780 Bool_t newcluster = kTRUE;
781 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
782 AliHLTTPCSignal_t lastcharge=0;
783 UInt_t bLastWasFalling=0;
784 Int_t newbin=-1;
785
786
787 if(fDeconvTime){
788 redo: //This is a goto.
789
790 if(newbin > -1){
791 //bin = newbin;
792 newbin = -1;
793 }
794
795 lastcharge=0;
796 bLastWasFalling = 0;
797 }
798
799 while(iResult>=0){ //Loop over time bins of current pad
800 // read all the values for one pad at once to calculate the base line
801 if (pCurrentPad) {
802 if (!pCurrentPad->IsStarted()) {
803 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
804 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
805 if ((pCurrentPad->StartEvent())>=0) {
806 do {
807 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
808 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
809 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
810 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
811 } while ((readValue = fDigitReader->Next())!=0);
812 }
813 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
814 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
815 //HLTDebug("no data available after zero suppression");
816 pCurrentPad->StopEvent();
817 pCurrentPad->ResetHistory();
818 break;
819 }
820 time=pCurrentPad->GetCurrentPosition();
821 if (time>pCurrentPad->GetSize()) {
822 HLTError("invalid time bin for pad");
823 break;
824 }
825 }
826 }
827 if (pCurrentPad) {
828 Float_t occupancy=pCurrentPad->GetOccupancy();
829 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
830 if ( occupancy < fOccupancyLimit ) {
831 charge = pCurrentPad->GetCorrectedData();
832 } else {
833 charge = 0;
834 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
835 }
836 } else {
837 charge = fDigitReader->GetSignal();
838 }
839 //HLTDebug("get next charge value: position %d charge %d", time, charge);
840
841
842 // CHARGE DEBUG
843 if (fDigitReader->GetRow() == 90){
844///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
845 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
846
847 }
848
849 if(time >= AliHLTTPCTransform::GetNTimeBins()){
850 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
851 iResult=-ERANGE;
852 }
853
854
855 //Get the current ADC-value
a38a7850 856 if(fDeconvTime){
857
858 //Check if the last pixel in the sequence is smaller than this
859 if(charge > lastcharge){
2a083ac4 860 if(bLastWasFalling){
a38a7850 861 newbin = 1;
862 break;
863 }
864 }
2a083ac4 865 else bLastWasFalling = 1; //last pixel was larger than this
a38a7850 866 lastcharge = charge;
867 }
868
869 //Sum the total charge of this sequence
870 seqcharge += charge;
871 seqaverage += time*charge;
872 seqerror += time*time*charge;
873
46b33a24 874 if (pCurrentPad) {
875
876 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
877 pCurrentPad->StopEvent();
878 pCurrentPad->ResetHistory();
879 if(readValue) {
880 newPad = fDigitReader->GetPad();
881 newTime = fDigitReader->GetTime();
882 newRow = fDigitReader->GetRow() + rowOffset;
883 }
884 break;
885 }
886
887 newPad=pCurrentPad->GetPadNumber();
888 newTime=pCurrentPad->GetCurrentPosition();
889 newRow=pCurrentPad->GetRowNumber();
890 } else {
a38a7850 891 readValue = fDigitReader->Next();
a38a7850 892 //Check where to stop:
893 if(!readValue) break; //No more value
894
895 newPad = fDigitReader->GetPad();
896 newTime = fDigitReader->GetTime();
897 newRow = fDigitReader->GetRow() + rowOffset;
46b33a24 898 }
a38a7850 899
900 if(newPad != pad)break; //new pad
901 if(newTime != time+1) break; //end of sequence
6b0fdc0e 902 if(iResult<0) break;
a38a7850 903
904 // pad = newpad; is equal
905 time = newTime;
906
907 }//end loop over sequence
908
46b33a24 909 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
910 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
911 if (seqcharge<=0) {
912 // with active zero suppression zero values are possible
913 continue;
914 }
a38a7850 915
916 //Calculate mean of sequence:
917 Int_t seqmean=0;
918 if(seqcharge)
919 seqmean = seqaverage/seqcharge;
920 else{
921 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
922 <<"Error in data given to the cluster finder"<<ENDLOG;
923 seqmean = 1;
924 seqcharge = 1;
925 }
926
927 //Calculate mean in pad direction:
928 Int_t padmean = seqcharge*pad;
929 Int_t paderror = pad*padmean;
930
a38a7850 931 //Compare with results on previous pad:
84645eb0 932 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
a38a7850 933
934 //dont merge sequences on the same pad twice
935 if(previousPt[p]->fLastMergedPad==pad) continue;
936
937 Int_t difference = seqmean - previousPt[p]->fMean;
938 if(difference < -fMatch) break;
939
940 if(difference <= fMatch){ //There is a match here!!
941 AliClusterData *local = previousPt[p];
942
943 if(fDeconvPad){
944 if(seqcharge > local->fLastCharge){
945 if(local->fChargeFalling){ //The previous pad was falling
946 break; //create a new cluster
947 }
948 }
949 else local->fChargeFalling = 1;
950 local->fLastCharge = seqcharge;
951 }
952
953 //Don't create a new cluster, because we found a match
954 newcluster = kFALSE;
955
956 //Update cluster on current pad with the matching one:
957 local->fTotalCharge += seqcharge;
958 local->fPad += padmean;
959 local->fPad2 += paderror;
960 local->fTime += seqaverage;
961 local->fTime2 += seqerror;
962 local->fMean = seqmean;
963 local->fFlags++; //means we have more than one pad
964 local->fLastMergedPad = pad;
965
966 currentPt[ncurrent] = local;
967 ncurrent++;
968
969 break;
970 } //Checking for match at previous pad
971 } //Loop over results on previous pad.
972
84645eb0 973 if(newcluster && ncurrent<kPadArraySize){
a38a7850 974 //Start a new cluster. Add it to the clusterlist, and update
975 //the list of pointers to clusters in current pad.
976 //current pad will be previous pad on next pad.
977
978 //Add to the clusterlist:
979 AliClusterData *tmp = &clusterlist[ntotal];
980 tmp->fTotalCharge = seqcharge;
981 tmp->fPad = padmean;
982 tmp->fPad2 = paderror;
983 tmp->fTime = seqaverage;
984 tmp->fTime2 = seqerror;
985 tmp->fMean = seqmean;
986 tmp->fFlags = 0; //flags for single pad clusters
987 tmp->fLastMergedPad = pad;
988
989 if(fDeconvPad){
990 tmp->fChargeFalling = 0;
991 tmp->fLastCharge = seqcharge;
992 }
993
994 //Update list of pointers to previous pad:
995 currentPt[ncurrent] = &clusterlist[ntotal];
996 ntotal++;
997 ncurrent++;
998 }
999
1000 if(fDeconvTime)
1001 if(newbin >= 0) goto redo;
1002
1003 // to prevent endless loop
1004 if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24 1005 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850 1006 break;
1007 }
1008
1009
1010 if(!readValue) break; //No more value
84645eb0 1011
1012 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
1013 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
1014 break;
1015 }
a38a7850 1016
1017 if(fCurrentRow != newRow){
1018 WriteClusters(ntotal,clusterlist);
1019
1020 lastpad = 123456789;
1021
1022 currentPt = pad2;
1023 previousPt = pad1;
1024 nprevious=0;
1025 ncurrent=0;
1026 ntotal=0;
1027
1028 fCurrentRow = newRow;
1029 }
1030
1031 pad = newPad;
1032 time = newTime;
db16520a 1033
a38a7850 1034 } // END while(readValue)
1035
1036 WriteClusters(ntotal,clusterlist);
1037
46b33a24 1038 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850 1039
1040} // ENDEND
1041
a912b63b 1042void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list){
1043 // see header file for class documentation
1044
a38a7850 1045 //write cluster to output pointer
8f8bf0af 1046 Int_t thisrow=-1,thissector=-1;
a38a7850 1047 UInt_t counter = fNClusters;
1048
1049 for(int j=0; j<nclusters; j++)
1050 {
d8b2f74a 1051
1052
1053
deba5d85 1054 if(!list[j].fFlags){
1055 if(fDoMC){
5ee8e602 1056 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
03efd623 1057 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1058 }
deba5d85 1059 }
1060 continue; //discard single pad clusters
1061 }
1062 if(list[j].fTotalCharge < fThreshold){
1063 if(fDoMC){
5ee8e602 1064 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
03efd623 1065 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1066 }
deba5d85 1067 }
1068 continue; //noise cluster
1069 }
a38a7850 1070 Float_t xyz[3];
1071 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1072 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1073 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1074 Float_t ftime2=fZErr*fZErr; //fixed given error
1075
738c049f 1076
aff6e981 1077
a912b63b 1078 if(fUnsorted){
1079 fCurrentRow=list[j].fRow;
8252a538 1080 }
8252a538 1081
a912b63b 1082
1083 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1084 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1085 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1086 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1087 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1088 sy2/=q2;
1089 if(sy2 < 0) {
1090 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1091 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1092 continue;
1093 } else {
1094 if(!fRawSP){
1095 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1096 if(sy2 != 0){
1097 fpad2*=0.108; //constants are from offline studies
1098 if(patch<2)
1099 fpad2*=2.07;
1100 }
1101 } else fpad2=sy2; //take the width not the error
6b15c309 1102 }
a912b63b 1103 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1104 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1105 sz2/=q2;
1106 if(sz2 < 0){
1107 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1108 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1109 continue;
1110 } else {
1111 if(!fRawSP){
1112 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1113 if(sz2 != 0) {
1114 ftime2 *= 0.169; //constants are from offline studies
1115 if(patch<2)
1116 ftime2 *= 1.77;
1117 }
1118 } else ftime2=sz2; //take the width, not the error
6b15c309 1119 }
1120 }
a912b63b 1121 if(fStdout==kTRUE)
1122 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
6b15c309 1123
a912b63b 1124 if(!fRawSP){
1125 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
5ee8e602 1126
1127 if(fOfflineTransform == NULL){
1128 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1129
1130 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1131 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1132 if(fNClusters >= fMaxNClusters)
1133 {
1134 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1135 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1136 return;
1137 }
a912b63b 1138
5ee8e602 1139 fSpacePointData[counter].fX = xyz[0];
1140 // fSpacePointData[counter].fY = xyz[1];
1141 if(fCurrentSlice<18){
1142 fSpacePointData[counter].fY = xyz[1];
1143 }
1144 else{
1145 fSpacePointData[counter].fY = -1*xyz[1];
1146 }
1147 fSpacePointData[counter].fZ = xyz[2];
6b15c309 1148 }
a912b63b 1149 else{
5ee8e602 1150 Double_t x[3]={thisrow,fpad+.5,ftime};
1151 Int_t iSector[1]={thissector};
1152 fOfflineTransform->Transform(x,iSector,0,1);
1153 fSpacePointData[counter].fX = x[0];
1154 fSpacePointData[counter].fY = x[1];
1155 fSpacePointData[counter].fZ = x[2];
a912b63b 1156 }
8252a538 1157
5ee8e602 1158 }
1159 else {
a912b63b 1160 fSpacePointData[counter].fX = fCurrentRow;
1161 fSpacePointData[counter].fY = fpad;
1162 fSpacePointData[counter].fZ = ftime;
8252a538 1163 }
a912b63b 1164
1165 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1166 fSpacePointData[counter].fPadRow = fCurrentRow;
1167 fSpacePointData[counter].fSigmaY2 = fpad2;
1168 fSpacePointData[counter].fSigmaZ2 = ftime2;
64defa03 1169
a912b63b 1170 fSpacePointData[counter].fQMax = list[j].fQMax;
2fdb1ae7 1171
a912b63b 1172 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1173 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
8252a538 1174
a912b63b 1175 Int_t patch=fCurrentPatch;
1176 if(patch==-1) patch=0; //never store negative patch number
1177 fSpacePointData[counter].fID = counter
1178 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
6b15c309 1179
a912b63b 1180#ifdef do_mc
1181 Int_t trackID[3];
1182 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1183
1184 fSpacePointData[counter].fTrackID[0] = trackID[0];
1185 fSpacePointData[counter].fTrackID[1] = trackID[1];
1186 fSpacePointData[counter].fTrackID[2] = trackID[2];
1187
1188#endif
1189
1190 fNClusters++;
1191 counter++;
8252a538 1192 }
a912b63b 1193}
8252a538 1194
a912b63b 1195// STILL TO FIX ----------------------------------------------------------------------------
8252a538 1196
a912b63b 1197#ifdef do_mc
1198void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID){
1199 // see header file for class documentation
1200
1201 //get mc id
1202 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
8252a538 1203
a912b63b 1204 trackID[0]=trackID[1]=trackID[2]=-2;
1205 for(Int_t i=fFirstRow; i<=fLastRow; i++){
1206 if(rowPt->fRow < (UInt_t)fCurrentRow){
1207 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1208 continue;
1209 }
1210 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1211 for(UInt_t j=0; j<rowPt->fNDigit; j++){
1212 Int_t cpad = digPt[j].fPad;
1213 Int_t ctime = digPt[j].fTime;
1214 if(cpad != pad) continue;
1215 if(ctime != time) continue;
fec96a66 1216
a912b63b 1217 trackID[0] = digPt[j].fTrackID[0];
1218 trackID[1] = digPt[j].fTrackID[1];
1219 trackID[2] = digPt[j].fTrackID[2];
1220
1221 break;
1222 }
1223 break;
8252a538 1224 }
01f43166 1225}
a912b63b 1226#endif
64defa03 1227
a912b63b 1228
1229
1230void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
8252a538 1231 // see header file for class documentation
8252a538 1232
01f43166 1233 //write cluster to output pointer
1234 Int_t thisrow,thissector;
1235 UInt_t counter = fNClusters;
1236
1237 for(int j=0; j<nclusters; j++)
1238 {
1239 if(!list[j].fFlags) continue; //discard single pad clusters
1240 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1241
1242 Float_t xyz[3];
1243 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1244 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1245 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1246 Float_t ftime2=fZErr*fZErr; //fixed given error
1247
d8b2f74a 1248
01f43166 1249 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1250 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1251 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1252 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1253 sy2/=q2;
1254 if(sy2 < 0) {
1255 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1256 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1257 continue;
1258 } else {
1259 if(!fRawSP){
1260 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1261 if(sy2 != 0){
1262 fpad2*=0.108; //constants are from offline studies
1263 if(patch<2)
1264 fpad2*=2.07;
1265 }
1266 } else fpad2=sy2; //take the width not the error
1267 }
1268 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1269 sz2/=q2;
1270 if(sz2 < 0){
1271 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1272 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1273 continue;
1274 } else {
1275 if(!fRawSP){
1276 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1277 if(sz2 != 0) {
1278 ftime2 *= 0.169; //constants are from offline studies
1279 if(patch<2)
1280 ftime2 *= 1.77;
1281 }
1282 } else ftime2=sz2; //take the width, not the error
1283 }
1284 }
1285 if(fStdout==kTRUE)
8252a538 1286 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1287
01f43166 1288 if(!fRawSP){
1289 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1290 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1291
1292 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1293 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1294 if(fNClusters >= fMaxNClusters)
1295 {
1296 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1297 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1298 return;
1299 }
1300
1301 fSpacePointData[counter].fX = xyz[0];
6b15c309 1302 // fSpacePointData[counter].fY = xyz[1];
1303 if(fCurrentSlice<18){
1304 fSpacePointData[counter].fY = xyz[1];
1305 }
1306 else{
1307 fSpacePointData[counter].fY = -1*xyz[1];
1308 }
01f43166 1309 fSpacePointData[counter].fZ = xyz[2];
1310
1311 } else {
1312 fSpacePointData[counter].fX = fCurrentRow;
1313 fSpacePointData[counter].fY = fpad;
1314 fSpacePointData[counter].fZ = ftime;
1315 }
1316
1317 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1318 fSpacePointData[counter].fPadRow = fCurrentRow;
1319 fSpacePointData[counter].fSigmaY2 = fpad2;
1320 fSpacePointData[counter].fSigmaZ2 = ftime2;
1321
a912b63b 1322 fSpacePointData[counter].fQMax = list[j].fQMax;
6b15c309 1323
01f43166 1324 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1325 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1326
1327 Int_t patch=fCurrentPatch;
1328 if(patch==-1) patch=0; //never store negative patch number
1329 fSpacePointData[counter].fID = counter
1330 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1331
1332#ifdef do_mc
1333 Int_t trackID[3];
1334 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1335
1336 fSpacePointData[counter].fTrackID[0] = trackID[0];
1337 fSpacePointData[counter].fTrackID[1] = trackID[1];
1338 fSpacePointData[counter].fTrackID[2] = trackID[2];
1339
01f43166 1340#endif
1341
1342 fNClusters++;
1343 counter++;
1344 }
1345}