2 // Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
4 //**************************************************************************
5 //* This file is property of and copyright by the ALICE HLT Project *
6 //* ALICE Experiment at CERN, All rights reserved. *
8 //* Primary Authors: Anders Vestbo, Constantin Loizides *
9 //* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
11 //* for The ALICE HLT Project. *
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 //**************************************************************************
22 // @file AliHLTTPCClusterFinder.cxx
23 // @author Kenneth Aamodt, Kalliopi Kanaki
25 // @brief Cluster Finder for the TPC
28 #include "AliHLTTPCDigitReader.h"
29 #include "AliHLTTPCRootTypes.h"
30 #include "AliHLTTPCLogging.h"
31 #include "AliHLTTPCClusterFinder.h"
32 #include "AliHLTTPCSpacePointData.h"
33 #include "AliHLTTPCMemHandler.h"
34 #include "AliHLTTPCPad.h"
38 #include "AliTPCcalibDB.h"
39 #include "AliTPCTransform.h"
40 #include "AliTPCParam.h"
46 ClassImp(AliHLTTPCClusterFinder)
48 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
50 fClustersHWAddressVector(),
52 fSpacePointData(NULL),
74 fVectorInitialized(kFALSE),
78 fNumberOfPadsInRow(NULL),
80 fRowOfFirstCandidate(0),
81 fDoPadSelection(kFALSE),
83 fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
84 fTotalChargeOfPreviousClusterCandidate(0),
85 fChargeOfCandidatesFalling(kFALSE),
89 fOfflineTransform(NULL),
90 fOfflineTPCParam( NULL ),
91 fOfflineTPCRecoParam(*AliTPCRecoParam::GetHLTParam()),
97 //uptate the transform class
99 fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
100 if(!fOfflineTransform){
101 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline transform not in AliTPCcalibDB.");
104 fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
107 fOfflineTPCParam = AliTPCcalibDB::Instance()->GetParameters();
108 if( !fOfflineTPCParam ){
109 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline TPC parameters not in AliTPCcalibDB.");
111 fOfflineTPCParam->Update();
112 fOfflineTPCParam->ReadGeoMatrices();
117 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder(){
118 // see header file for class documentation
121 if(fVectorInitialized){
122 DeInitializePadArray();
124 if(fNumberOfPadsInRow){
125 delete [] fNumberOfPadsInRow;
126 fNumberOfPadsInRow=NULL;
130 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints){
131 // see header file for class documentation
135 fMaxNClusters = nmaxpoints;
136 fCurrentSlice = slice;
137 fCurrentPatch = patch;
138 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
139 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
142 fClustersMCInfo.clear();
144 fClusterMCVector.clear();
147 void AliHLTTPCClusterFinder::InitializePadArray(){
148 // see header file for class documentation
150 if(fCurrentPatch>5||fCurrentPatch<0){
151 HLTFatal("Patch is not set");
155 HLTDebug("Patch number=%d",fCurrentPatch);
157 fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
158 fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
160 fNumberOfRows=fLastRow-fFirstRow+1;
161 fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
163 memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
165 fRowPadVector.clear();
167 for(UInt_t i=0;i<fNumberOfRows;i++){
168 fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
169 AliHLTTPCPadVector tmpRow;
170 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
171 AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
173 tmpRow.push_back(tmpPad);
175 fRowPadVector.push_back(tmpRow);
177 fVectorInitialized=kTRUE;
180 Int_t AliHLTTPCClusterFinder::DeInitializePadArray(){
181 // see header file for class documentation
183 if( fVectorInitialized ){
184 for(UInt_t i=0;i<fNumberOfRows;i++){
185 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
186 delete fRowPadVector[i][j];
187 fRowPadVector[i][j]=NULL;
189 fRowPadVector[i].clear();
191 fRowPadVector.clear();
192 delete[] fNumberOfPadsInRow;
193 fNumberOfPadsInRow = 0;
195 fVectorInitialized=kFALSE;
200 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt){
201 // see header file for class documentation
202 //set pointer to output
203 fSpacePointData = pt;
207 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
208 // see header file for class documentation
210 fPtr = (UChar_t*)ptr;
213 if(!fVectorInitialized){
214 InitializePadArray();
217 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
218 HLTError("failed setting up digit reader (InitBlock)");
222 while(fDigitReader->NextChannel()){
223 UInt_t row=fDigitReader->GetRow();
224 UInt_t pad=fDigitReader->GetPad();
226 if(row>=fRowPadVector.size()){
227 HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
230 if(pad>=fRowPadVector[row].size()){
231 HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
235 while(fDigitReader->NextBunch()){
236 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
237 UInt_t time = fDigitReader->GetTime();
238 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
239 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
240 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
241 // The same is true for the function ReadDataUnsortedDeconvoluteTime() below.
242 // In addition the signals are organized in the opposite direction
244 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
245 AliHLTTPCClusters candidate;
246 for(Int_t i=fDigitReader->GetBunchSize()-1;i>=0;i--){
247 candidate.fTotalCharge+=bunchData[i];
248 candidate.fTime += time*bunchData[i];
249 candidate.fTime2 += time*time*bunchData[i];
250 if(bunchData[i]>candidate.fQMax){
251 candidate.fQMax=bunchData[i];
255 if(candidate.fTotalCharge>0){
256 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
257 candidate.fPad=candidate.fTotalCharge*pad;
258 candidate.fPad2=candidate.fPad*pad;
259 candidate.fLastMergedPad=pad;
260 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
262 if(fRowPadVector[row][pad] != NULL){
263 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
267 const UInt_t *bunchData= fDigitReader->GetSignals();
268 AliHLTTPCClusters candidate;
269 const AliHLTTPCDigitData* digits = NULL;
270 if(fDoMC && (digits = fDigitReader->GetBunchDigits())!=NULL){
271 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
272 candidate.fTotalCharge+=bunchData[i];
273 candidate.fTime += time*bunchData[i];
274 candidate.fTime2 += time*time*bunchData[i];
275 if(bunchData[i]>candidate.fQMax){
276 candidate.fQMax=bunchData[i];
278 fMCDigits.push_back(digits[i]);
283 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
284 candidate.fTotalCharge+=bunchData[i];
285 candidate.fTime += time*bunchData[i];
286 candidate.fTime2 += time*time*bunchData[i];
287 if(bunchData[i]>candidate.fQMax){
288 candidate.fQMax=bunchData[i];
293 if(candidate.fTotalCharge>0){
294 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
295 candidate.fPad=candidate.fTotalCharge*pad;
296 candidate.fPad2=candidate.fPad*pad;
297 candidate.fLastMergedPad=pad;
298 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
300 if(fRowPadVector[row][pad] != NULL){
301 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
303 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
314 void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size){
315 // see header file for class documentation
318 fPtr = (UChar_t*)ptr;
321 if(!fVectorInitialized){
322 InitializePadArray();
325 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
326 HLTError("failed setting up digit reader (InitBlock)");
330 while(fDigitReader->NextChannel()){
331 UInt_t row=fDigitReader->GetRow();
332 UInt_t pad=fDigitReader->GetPad();
334 while(fDigitReader->NextBunch()){
335 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
336 UInt_t time = fDigitReader->GetTime();
337 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
338 Int_t indexInBunchData=0;
339 Bool_t moreDataInBunch=kFALSE;
341 Bool_t signalFalling=kFALSE;
343 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
344 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
345 // The same is true for the function ReadDataUnsorted() above.
346 // In addition the signals are organized in the opposite direction
348 indexInBunchData = fDigitReader->GetBunchSize()-1;
349 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
352 AliHLTTPCClusters candidate;
353 //for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
354 for(Int_t i=indexInBunchData;i>=0;i--){
356 // Checks if one need to deconvolute the signals
357 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
358 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
359 moreDataInBunch=kTRUE;
365 // Checks if the signal is 0, then quit processing the data.
366 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
367 moreDataInBunch=kTRUE;
372 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
375 candidate.fTotalCharge+=bunchData[i];
376 candidate.fTime += time*bunchData[i];
377 candidate.fTime2 += time*time*bunchData[i];
378 if(bunchData[i]>candidate.fQMax){
379 candidate.fQMax=bunchData[i];
381 prevSignal=bunchData[i];
385 if(candidate.fTotalCharge>0){
386 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
387 candidate.fPad=candidate.fTotalCharge*pad;
388 candidate.fPad2=candidate.fPad*pad;
389 candidate.fLastMergedPad=pad;
390 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
392 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
393 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
394 moreDataInBunch=kFALSE;
396 }while(moreDataInBunch);
399 const UInt_t *bunchData= fDigitReader->GetSignals();
401 AliHLTTPCClusters candidate;
402 const AliHLTTPCDigitData* digits = fDigitReader->GetBunchDigits();
403 if(fDoMC) fMCDigits.clear();
405 for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
406 // Checks if one need to deconvolute the signals
407 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
408 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
409 moreDataInBunch=kTRUE;
415 // Checks if the signal is 0, then quit processing the data.
416 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
417 moreDataInBunch=kTRUE;
422 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
425 candidate.fTotalCharge+=bunchData[i];
426 candidate.fTime += time*bunchData[i];
427 candidate.fTime2 += time*time*bunchData[i];
428 if(bunchData[i]>candidate.fQMax){
429 candidate.fQMax=bunchData[i];
431 if( fDoMC ) fMCDigits.push_back(digits[i]);
433 prevSignal=bunchData[i];
437 if(candidate.fTotalCharge>0){
438 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
439 candidate.fPad=candidate.fTotalCharge*pad;
440 candidate.fPad2=candidate.fPad*pad;
441 candidate.fLastMergedPad=pad;
442 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
444 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
446 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
449 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
450 moreDataInBunch=kFALSE;
452 }while(moreDataInBunch);
460 Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
461 // see header file for class documentation
463 //Checking if we have a match on the next pad
464 for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
465 if(nextPad->fUsedClusterCandidates[candidateNumber] == 1){
468 AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
469 // if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
471 if( abs((Int_t)(cluster->fMean - candidate->fMean)) <= fTimeMeanDiff ){
473 if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
474 fChargeOfCandidatesFalling=kTRUE;
476 if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
480 cluster->fMean=candidate->fMean;
481 cluster->fTotalCharge+=candidate->fTotalCharge;
482 cluster->fTime += candidate->fTime;
483 cluster->fTime2 += candidate->fTime2;
484 cluster->fPad+=candidate->fPad;
485 cluster->fPad2+=candidate->fPad2;
486 cluster->fLastMergedPad=candidate->fPad;
487 if(candidate->fQMax>cluster->fQMax){
488 cluster->fQMax=candidate->fQMax;
491 FillMCClusterVector(nextPad->GetCandidateDigits(candidateNumber));
495 UInt_t rowNo = nextPad->GetRowNumber();
496 UInt_t padNo = nextPad->GetPadNumber();
498 fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
499 fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
501 fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
502 fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
503 fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
504 fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
507 //setting the matched pad to used
508 nextPad->fUsedClusterCandidates[candidateNumber]=1;
510 if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
511 nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
512 ComparePads(nextPad,cluster,nextPadToRead);
522 Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
523 // see header file for class documentation
526 for(UInt_t row=0;row<fNumberOfRows;row++){
527 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
528 if(fRowPadVector[row][pad]->fSelectedPad){
529 if(counter<maxHWadd){
530 hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
534 HLTWarning("To many hardwareaddresses, skip adding");
544 Int_t AliHLTTPCClusterFinder::FillOutputMCInfo(AliHLTTPCClusterFinder::ClusterMCInfo * outputMCInfo, Int_t maxNumberOfClusterMCInfo){
545 // see header file for class documentation
549 for(UInt_t mc=0;mc<fClustersMCInfo.size();mc++){
550 if(counter<maxNumberOfClusterMCInfo){
551 outputMCInfo[counter] = fClustersMCInfo[mc];
555 HLTWarning("To much MCInfo has been added (no more space), skip adding");
561 void AliHLTTPCClusterFinder::FindClusters(){
562 // see header file for function documentation
564 AliHLTTPCClusters* tmpCandidate=NULL;
565 for(UInt_t row=0;row<fNumberOfRows;row++){
566 fRowOfFirstCandidate=row;
567 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
568 AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
569 for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
570 if(tmpPad->fUsedClusterCandidates[candidate]){
573 tmpCandidate=&tmpPad->fClusterCandidates[candidate];
574 UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
577 fClusterMCVector.clear();
578 FillMCClusterVector(tmpPad->GetCandidateDigits(candidate));
581 ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
582 if(tmpCandidate->fTotalCharge>tmpTotalCharge){
584 fClusters.push_back(*tmpCandidate);
586 //sort the vector (large->small) according to weight and remove elements above 2 (keep 0 1 and 2)
587 sort(fClusterMCVector.begin(),fClusterMCVector.end(), MCWeight::CompareWeights );
588 ClusterMCInfo tmpClusterMCInfo;
594 if(fClusterMCVector.size()>0){
595 tmpClusterMCInfo.fClusterID[0]=fClusterMCVector.at(0);
598 tmpClusterMCInfo.fClusterID[0]=zeroMC;
601 if(fClusterMCVector.size()>1){
602 tmpClusterMCInfo.fClusterID[1]=fClusterMCVector.at(1);
605 tmpClusterMCInfo.fClusterID[1]=zeroMC;
608 if(fClusterMCVector.size()>2){
609 tmpClusterMCInfo.fClusterID[2]=fClusterMCVector.at(2);
612 tmpClusterMCInfo.fClusterID[2]=zeroMC;
615 fClustersMCInfo.push_back(tmpClusterMCInfo);
620 tmpPad->ClearCandidates();
622 fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
625 HLTInfo("Found %d clusters.",fClusters.size());
627 //TODO: Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
629 AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
630 for(unsigned int i=0;i<fClusters.size();i++){
631 clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
632 clusterlist[i].fPad = fClusters[i].fPad;
633 clusterlist[i].fPad2 = fClusters[i].fPad2;
634 clusterlist[i].fTime = fClusters[i].fTime;
635 clusterlist[i].fTime2 = fClusters[i].fTime2;
636 clusterlist[i].fMean = fClusters[i].fMean;
637 clusterlist[i].fFlags = fClusters[i].fFlags;
638 clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
639 clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
640 clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
641 clusterlist[i].fRow = fClusters[i].fRowNumber;
642 clusterlist[i].fQMax = fClusters[i].fQMax;
645 WriteClusters(fClusters.size(),clusterlist);
646 delete [] clusterlist;
648 if( fReleaseMemory ) DeInitializePadArray();// call this when the -releaseMemory flag is set
652 Bool_t AliHLTTPCClusterFinder::UpdateCalibDB(){
655 AliTPCcalibDB::Instance()->Update();
659 //uptate the transform class
661 fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
662 if(!fOfflineTransform){
663 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline transform not in AliTPCcalibDB.");
667 fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
670 fOfflineTPCParam = AliTPCcalibDB::Instance()->GetParameters();
671 if( !fOfflineTPCParam ){
672 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline TPC parameters not in AliTPCcalibDB.");
675 fOfflineTPCParam->Update();
676 fOfflineTPCParam->ReadGeoMatrices();
682 //---------------------------------- Under this line the old sorted clusterfinder functions can be found --------------------------------
685 void AliHLTTPCClusterFinder::PrintClusters(){
686 // see header file for class documentation
688 for(size_t i=0;i<fClusters.size();i++){
689 HLTInfo("Cluster number: %d",i);
690 HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
691 HLTInfo("Total Charge: %d",fClusters[i].fTotalCharge);
692 HLTInfo("fPad: %d",fClusters[i].fPad);
693 HLTInfo("PadError: %d",fClusters[i].fPad2);
694 HLTInfo("TimeMean: %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
695 HLTInfo("TimeError: %d",fClusters[i].fTime2);
696 HLTInfo("EndOfCluster:");
700 void AliHLTTPCClusterFinder::FillMCClusterVector(vector<AliHLTTPCDigitData> *digitData){
701 // see header file for class documentation
702 if( !digitData ) return;
703 for(UInt_t d=0;d<digitData->size();d++){
704 Int_t nIDsInDigit = (digitData->at(d).fTrackID[0]>=0) + (digitData->at(d).fTrackID[1]>=0) + (digitData->at(d).fTrackID[2]>=0);
705 for(Int_t id=0; id<3; id++){
706 if(digitData->at(d).fTrackID[id]>=0){
707 Bool_t matchFound = kFALSE;
709 mc.fMCID = digitData->at(d).fTrackID[id];
710 mc.fWeight = ((Float_t)digitData->at(d).fCharge)/nIDsInDigit;
711 for(UInt_t i=0;i<fClusterMCVector.size();i++){
712 if(mc.fMCID == fClusterMCVector.at(i).fMCID){
713 fClusterMCVector.at(i).fWeight += mc.fWeight;
717 if(matchFound == kFALSE){
718 fClusterMCVector.push_back(mc);
726 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
728 fPtr = (UChar_t*)ptr;
732 void AliHLTTPCClusterFinder::ProcessDigits(){
733 // see header file for class documentation
736 bool readValue = true;
739 UShort_t time=0,newTime=0;
740 UInt_t pad=0,newPad=0;
741 AliHLTTPCSignal_t charge=0;
745 // initialize block for reading packed data
746 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
747 if (iResult<0) return;
749 readValue = fDigitReader->Next();
751 // Matthias 08.11.2006 the following return would cause termination without writing the
752 // ClusterData and thus would block the component. I just want to have the commented line
753 // here for information
754 //if (!readValue)return;
756 pad = fDigitReader->GetPad();
757 time = fDigitReader->GetTime();
758 fCurrentRow = fDigitReader->GetRow();
760 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
761 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
763 fCurrentRow += rowOffset;
765 UInt_t lastpad = 123456789;
766 const UInt_t kPadArraySize=5000;
767 const UInt_t kClusterListSize=10000;
768 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
769 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
770 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
772 AliClusterData **currentPt; //List of pointers to the current pad
773 AliClusterData **previousPt; //List of pointers to the previous pad
776 UInt_t nprevious=0,ncurrent=0,ntotal=0;
778 /* quick implementation of baseline calculation and zero suppression
779 open a pad object for each pad and delete it after processing.
780 later a list of pad objects with base line history can be used
781 The whole thing only works if we really get unprocessed raw data, if
782 the data is already zero suppressed, there might be gaps in the time
785 Int_t gatingGridOffset=50;
787 gatingGridOffset=fFirstTimeBin;
789 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
790 // just to make later conversion to a list of objects easier
791 AliHLTTPCPad* pCurrentPad=NULL;
793 if (fSignalThreshold>=0) {
794 pCurrentPad=&baseline;
795 baseline.SetThreshold(fSignalThreshold);
798 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
805 if(currentPt == pad2){
813 nprevious = ncurrent;
815 if(pad != lastpad+1){
816 //this happens if there is a pad with no signal.
817 nprevious = ncurrent = 0;
822 Bool_t newcluster = kTRUE;
823 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
824 AliHLTTPCSignal_t lastcharge=0;
825 UInt_t bLastWasFalling=0;
830 redo: //This is a goto.
841 while(iResult>=0){ //Loop over time bins of current pad
842 // read all the values for one pad at once to calculate the base line
844 if (!pCurrentPad->IsStarted()) {
845 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
846 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
847 if ((pCurrentPad->StartEvent())>=0) {
849 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
850 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
851 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
852 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
853 } while ((readValue = fDigitReader->Next())!=0);
855 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
856 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
857 //HLTDebug("no data available after zero suppression");
858 pCurrentPad->StopEvent();
859 pCurrentPad->ResetHistory();
862 time=pCurrentPad->GetCurrentPosition();
863 if (time>pCurrentPad->GetSize()) {
864 HLTError("invalid time bin for pad");
870 Float_t occupancy=pCurrentPad->GetOccupancy();
871 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
872 if ( occupancy < fOccupancyLimit ) {
873 charge = pCurrentPad->GetCorrectedData();
876 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
879 charge = fDigitReader->GetSignal();
881 //HLTDebug("get next charge value: position %d charge %d", time, charge);
885 if (fDigitReader->GetRow() == 90){
886 ///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
887 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
891 if(time >= AliHLTTPCTransform::GetNTimeBins()){
892 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
897 //Get the current ADC-value
900 //Check if the last pixel in the sequence is smaller than this
901 if(charge > lastcharge){
907 else bLastWasFalling = 1; //last pixel was larger than this
911 //Sum the total charge of this sequence
913 seqaverage += time*charge;
914 seqerror += time*time*charge;
918 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
919 pCurrentPad->StopEvent();
920 pCurrentPad->ResetHistory();
922 newPad = fDigitReader->GetPad();
923 newTime = fDigitReader->GetTime();
924 newRow = fDigitReader->GetRow() + rowOffset;
929 newPad=pCurrentPad->GetPadNumber();
930 newTime=pCurrentPad->GetCurrentPosition();
931 newRow=pCurrentPad->GetRowNumber();
933 readValue = fDigitReader->Next();
934 //Check where to stop:
935 if(!readValue) break; //No more value
937 newPad = fDigitReader->GetPad();
938 newTime = fDigitReader->GetTime();
939 newRow = fDigitReader->GetRow() + rowOffset;
942 if(newPad != pad)break; //new pad
943 if(newTime != time+1) break; //end of sequence
946 // pad = newpad; is equal
949 }//end loop over sequence
951 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
952 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
954 // with active zero suppression zero values are possible
958 //Calculate mean of sequence:
961 seqmean = seqaverage/seqcharge;
963 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
964 <<"Error in data given to the cluster finder"<<ENDLOG;
969 //Calculate mean in pad direction:
970 Int_t padmean = seqcharge*pad;
971 Int_t paderror = pad*padmean;
973 //Compare with results on previous pad:
974 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
976 //dont merge sequences on the same pad twice
977 if(previousPt[p]->fLastMergedPad==pad) continue;
979 Int_t difference = seqmean - previousPt[p]->fMean;
980 if(difference < -fMatch) break;
982 if(difference <= fMatch){ //There is a match here!!
983 AliClusterData *local = previousPt[p];
986 if(seqcharge > local->fLastCharge){
987 if(local->fChargeFalling){ //The previous pad was falling
988 break; //create a new cluster
991 else local->fChargeFalling = 1;
992 local->fLastCharge = seqcharge;
995 //Don't create a new cluster, because we found a match
998 //Update cluster on current pad with the matching one:
999 local->fTotalCharge += seqcharge;
1000 local->fPad += padmean;
1001 local->fPad2 += paderror;
1002 local->fTime += seqaverage;
1003 local->fTime2 += seqerror;
1004 local->fMean = seqmean;
1005 local->fFlags++; //means we have more than one pad
1006 local->fLastMergedPad = pad;
1008 currentPt[ncurrent] = local;
1012 } //Checking for match at previous pad
1013 } //Loop over results on previous pad.
1015 if(newcluster && ncurrent<kPadArraySize){
1016 //Start a new cluster. Add it to the clusterlist, and update
1017 //the list of pointers to clusters in current pad.
1018 //current pad will be previous pad on next pad.
1020 //Add to the clusterlist:
1021 AliClusterData *tmp = &clusterlist[ntotal];
1022 tmp->fTotalCharge = seqcharge;
1023 tmp->fPad = padmean;
1024 tmp->fPad2 = paderror;
1025 tmp->fTime = seqaverage;
1026 tmp->fTime2 = seqerror;
1027 tmp->fMean = seqmean;
1028 tmp->fFlags = 0; //flags for single pad clusters
1029 tmp->fLastMergedPad = pad;
1032 tmp->fChargeFalling = 0;
1033 tmp->fLastCharge = seqcharge;
1036 //Update list of pointers to previous pad:
1037 currentPt[ncurrent] = &clusterlist[ntotal];
1043 if(newbin >= 0) goto redo;
1045 // to prevent endless loop
1046 if(time >= AliHLTTPCTransform::GetNTimeBins()){
1047 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
1052 if(!readValue) break; //No more value
1054 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
1055 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
1059 if(fCurrentRow != newRow){
1060 WriteClusters(ntotal,clusterlist);
1062 lastpad = 123456789;
1070 fCurrentRow = newRow;
1076 } // END while(readValue)
1078 WriteClusters(ntotal,clusterlist);
1080 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
1084 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list){
1085 // see header file for class documentation
1087 //write cluster to output pointer
1088 Int_t thisrow=-1,thissector=-1;
1089 UInt_t counter = fNClusters;
1091 for(int j=0; j<nclusters; j++)
1096 if(!list[j].fFlags){
1098 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
1099 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1102 continue; //discard single pad clusters
1104 if(list[j].fTotalCharge < fThreshold){
1106 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
1107 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1110 continue; //noise cluster
1113 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1114 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1115 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1116 Float_t ftime2=fZErr*fZErr; //fixed given error
1121 fCurrentRow=list[j].fRow;
1125 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1126 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1127 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1128 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1129 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1132 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1133 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1137 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1139 fpad2*=0.108; //constants are from offline studies
1143 } else fpad2=sy2; //take the width not the error
1145 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1146 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1149 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1150 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1154 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1156 ftime2 *= 0.169; //constants are from offline studies
1160 } else ftime2=sz2; //take the width, not the error
1164 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1167 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1169 if(fOfflineTransform == NULL){
1170 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1172 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1173 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1174 if(fNClusters >= fMaxNClusters)
1176 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1177 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1181 fSpacePointData[counter].fX = xyz[0];
1182 // fSpacePointData[counter].fY = xyz[1];
1183 if(fCurrentSlice<18){
1184 fSpacePointData[counter].fY = xyz[1];
1187 fSpacePointData[counter].fY = -1*xyz[1];
1189 fSpacePointData[counter].fZ = xyz[2];
1192 Double_t x[3]={thisrow,fpad+.5,ftime};
1193 Int_t iSector[1]={thissector};
1194 fOfflineTransform->Transform(x,iSector,0,1);
1195 double y[3] = {x[0], x[1], x[2] };
1197 if( fOfflineTPCParam && thissector<fOfflineTPCParam->GetNSector() ){
1198 TGeoHMatrix *alignment = fOfflineTPCParam->GetClusterMatrix( thissector );
1199 if ( alignment ) alignment->LocalToMaster( x, y);
1202 fSpacePointData[counter].fX = y[0];
1203 fSpacePointData[counter].fY = y[1];
1204 fSpacePointData[counter].fZ = y[2];
1209 fSpacePointData[counter].fX = fCurrentRow;
1210 fSpacePointData[counter].fY = fpad;
1211 fSpacePointData[counter].fZ = ftime;
1214 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1215 fSpacePointData[counter].fPadRow = fCurrentRow;
1216 fSpacePointData[counter].fSigmaY2 = fpad2;
1217 fSpacePointData[counter].fSigmaZ2 = ftime2;
1219 fSpacePointData[counter].fQMax = list[j].fQMax;
1221 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1222 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1224 Int_t patch=fCurrentPatch;
1225 if(patch==-1) patch=0; //never store negative patch number
1226 fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
1230 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1232 fSpacePointData[counter].fTrackID[0] = trackID[0];
1233 fSpacePointData[counter].fTrackID[1] = trackID[1];
1234 fSpacePointData[counter].fTrackID[2] = trackID[2];
1243 // STILL TO FIX ----------------------------------------------------------------------------
1246 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID) const {
1247 // see header file for class documentation
1250 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
1252 trackID[0]=trackID[1]=trackID[2]=-2;
1253 for(Int_t i=fFirstRow; i<=fLastRow; i++){
1254 if(rowPt->fRow < (UInt_t)fCurrentRow){
1255 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1258 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1259 for(UInt_t j=0; j<rowPt->fNDigit; j++){
1260 Int_t cpad = digPt[j].fPad;
1261 Int_t ctime = digPt[j].fTime;
1262 if(cpad != pad) continue;
1263 if(ctime != time) continue;
1265 trackID[0] = digPt[j].fTrackID[0];
1266 trackID[1] = digPt[j].fTrackID[1];
1267 trackID[2] = digPt[j].fTrackID[2];
1278 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
1279 // see header file for class documentation
1281 //write cluster to output pointer
1282 Int_t thisrow,thissector;
1283 UInt_t counter = fNClusters;
1285 for(int j=0; j<nclusters; j++)
1287 if(!list[j].fFlags) continue; //discard single pad clusters
1288 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1291 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1292 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1293 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1294 Float_t ftime2=fZErr*fZErr; //fixed given error
1297 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1298 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1299 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1300 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1303 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1304 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1308 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1310 fpad2*=0.108; //constants are from offline studies
1314 } else fpad2=sy2; //take the width not the error
1316 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1319 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1320 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1324 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1326 ftime2 *= 0.169; //constants are from offline studies
1330 } else ftime2=sz2; //take the width, not the error
1334 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1337 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1338 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1340 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1341 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1342 if(fNClusters >= fMaxNClusters)
1344 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1345 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1349 fSpacePointData[counter].fX = xyz[0];
1350 // fSpacePointData[counter].fY = xyz[1];
1351 if(fCurrentSlice<18){
1352 fSpacePointData[counter].fY = xyz[1];
1355 fSpacePointData[counter].fY = -1*xyz[1];
1357 fSpacePointData[counter].fZ = xyz[2];
1360 fSpacePointData[counter].fX = fCurrentRow;
1361 fSpacePointData[counter].fY = fpad;
1362 fSpacePointData[counter].fZ = ftime;
1365 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1366 fSpacePointData[counter].fPadRow = fCurrentRow;
1367 fSpacePointData[counter].fSigmaY2 = fpad2;
1368 fSpacePointData[counter].fSigmaZ2 = ftime2;
1370 fSpacePointData[counter].fQMax = list[j].fQMax;
1372 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1373 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1375 Int_t patch=fCurrentPatch;
1376 if(patch==-1) patch=0; //never store negative patch number
1377 fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
1381 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1383 fSpacePointData[counter].fTrackID[0] = trackID[0];
1384 fSpacePointData[counter].fTrackID[1] = trackID[1];
1385 fSpacePointData[counter].fTrackID[2] = trackID[2];