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