HLT/TPCLib/AliHLTTPCKryptonClusterFinder.cxx

   1 // $Id$
   2 //**************************************************************************
   3 //* This file is property of and copyright by the ALICE HLT Project        *
   4 //* ALICE Experiment at CERN, All rights reserved.                         *
   5 //*                                                                        *
   6 //* Primary Authors: Kenneth Aamodt <kenneth.aamodt@student.uib.no>        *
   7 //*                  for The ALICE HLT Project.                            *
   8 //*                                                                        *
   9 //* Permission to use, copy, modify and distribute this software and its   *
  10 //* documentation strictly for non-commercial purposes is hereby granted   *
  11 //* without fee, provided that the above copyright notice appears in all   *
  12 //* copies and that both the copyright notice and this permission notice   *
  13 //* appear in the supporting documentation. The authors make no claims     *
  14 //* about the suitability of this software for any purpose. It is          *
  15 //* provided "as is" without express or implied warranty.                  *
  16 //**************************************************************************
  17
  18 /** @file    AliHLTTPCKryptonClusterFinder.cxx
  19     @author  Kenneth Aamodt, Kalliopi Kanaki
  20     @date
  21     @brief  Krypton Cluster Finder for the TPC
  22 */
  23
  24 #include "AliHLTTPCDigitReader.h"
  25 #include "AliHLTTPCKryptonClusterFinder.h"
  26 #include "AliHLTTPCTransform.h"
  27 #include "AliHLTTPCPad.h"
  28 #include "AliHLTTPCClusters.h"
  29 #include "TFile.h"
  30 #include "AliHLTTPCSpacePointData.h"
  31
  32
  33 #if __GNUC__ >= 3
  34 using namespace std;
  35 #endif
  36
  37 ClassImp(AliHLTTPCKryptonClusterFinder)
  38
  39 AliHLTTPCKryptonClusterFinder::AliHLTTPCKryptonClusterFinder()
  40   : AliHLTTPCClusterFinder(),
  41     fHWAddressVector(),
  42     fTimebinsInBunch(),
  43     fIndexOfBunchStart(),
  44     fMaxQOfCluster(0),
  45     fSelectionMinRowNumber(0),
  46     fSelectionMaxRowNumber(0),
  47     fNKryptonClusters(0),
  48     fMaxOutputSize(0)
  49 {
  50   //constructor
  51 }
  52
  53 AliHLTTPCKryptonClusterFinder::~AliHLTTPCKryptonClusterFinder(){
  54 }
  55
  56 void AliHLTTPCKryptonClusterFinder::ReBunch(const UInt_t *bunchData , Int_t bunchSize){
  57   fTimebinsInBunch.clear();
  58   fIndexOfBunchStart.clear();
  59   Bool_t newBunch=kTRUE;
  60   Int_t currentBunchNumber=-1;
  61   for(Int_t i=0;i<bunchSize;i++){
  62     if(newBunch){
  63       if(bunchData[i]>5){
  64         fIndexOfBunchStart.push_back(i);
  65         fTimebinsInBunch.push_back(1);
  66         currentBunchNumber++;
  67         newBunch=kFALSE;
  68       }
  69     }
  70     else{
  71       if(bunchData[i]>0){
  72         fTimebinsInBunch[currentBunchNumber]++;
  73       }
  74       else{
  75         newBunch=kTRUE;
  76       }
  77     }
  78   }
  79 }
  80
  81 void AliHLTTPCKryptonClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
  82 {
  83   //set input pointer
  84   fPtr = (UChar_t*)ptr;
  85   fSize = size;
  86
  87   if(!fVectorInitialized){
  88     InitializePadArray();
  89    }
  90
  91   fHWAddressVector.clear();
  92   fNKryptonClusters=0;
  93
  94   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
  95
  96   while(fDigitReader->NextChannel()){
  97     UInt_t row=fDigitReader->GetRow();
  98     UInt_t pad=fDigitReader->GetPad();
  99
 100     fRowPadVector[row][pad]->ClearCandidates();
 101
 102     while(fDigitReader->NextBunch()){
 103       if(fDigitReader->GetBunchSize()>1){
 104         const UInt_t *bunchData= fDigitReader->GetSignals();
 105
 106         ReBunch(bunchData,fDigitReader->GetBunchSize());
 107         Int_t rebunchCount=fIndexOfBunchStart.size();
 108
 109         for(Int_t r=0;r<rebunchCount;r++){
 110           UInt_t time = fDigitReader->GetTime()+fIndexOfBunchStart[r];
 111           AliHLTTPCClusters candidate;
 112           candidate.fTotalCharge=0;
 113           candidate.fQMax=0;
 114           if(fTimebinsInBunch[r]>2){
 115             for(Int_t i=0;i<fTimebinsInBunch[r];i++){
 116               candidate.fTotalCharge+=bunchData[i + fIndexOfBunchStart[r]];
 117               candidate.fTime += time*bunchData[i + fIndexOfBunchStart[r]];
 118               candidate.fTime2 += time*time*bunchData[i + fIndexOfBunchStart[r]];
 119               if(bunchData[i + fIndexOfBunchStart[r]]>candidate.fQMax){
 120                 candidate.fQMax = bunchData[i + fIndexOfBunchStart[r]];
 121               }
 122               time++;
 123             }
 124             if(candidate.fTotalCharge>0){
 125               candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 126               candidate.fPad=candidate.fTotalCharge*pad;
 127               candidate.fPad2=candidate.fPad*pad;
 128               candidate.fLastMergedPad=pad;
 129               candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 130             }
 131             if(candidate.fTotalCharge>10 && candidate.fMean<924){
 132               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 133             }
 134           }
 135         }
 136       }
 137     }
 138   }
 139 }
 140
 141 Bool_t AliHLTTPCKryptonClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
 142   //Checking if we have a match on the next pad
 143   for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
 144     AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
 145     if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
 146       cluster->fMean=candidate->fMean;
 147       cluster->fTotalCharge+=candidate->fTotalCharge;
 148       cluster->fTime += candidate->fTime;
 149       cluster->fTime2 += candidate->fTime2;
 150       cluster->fPad+=candidate->fPad;
 151       cluster->fPad2=candidate->fPad2;
 152       cluster->fLastMergedPad=nextPad->GetPadNumber();
 153       if(candidate->fQMax>cluster->fQMax){
 154         cluster->fQMax = candidate->fQMax;
 155       }
 156
 157       //setting the matched pad to used
 158       nextPad->fUsedClusterCandidates[candidateNumber]=1;
 159       nextPadToRead++;
 160       if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
 161         nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
 162         ComparePads(nextPad,cluster,nextPadToRead);
 163       }
 164       else{
 165         return kFALSE;
 166       }
 167     }
 168     else{
 169       return kFALSE;
 170     }
 171   }
 172   return kFALSE;
 173 }
 174
 175 void AliHLTTPCKryptonClusterFinder::FindRowClusters()
 176 {
 177   // see header file for function documentation
 178
 179   AliHLTTPCClusters* tmpCandidate=NULL;
 180   for(UInt_t row=5;row<fNumberOfRows-5;row++){
 181     fRowOfFirstCandidate=row;
 182     for(UInt_t pad=5;pad<fNumberOfPadsInRow[row]-1-5;pad++){
 183       AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
 184       for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
 185         if(tmpPad->fUsedClusterCandidates[candidate]){
 186           continue;
 187         }
 188         if((Int_t)tmpPad->fClusterCandidates[candidate].fMean<100 || (Int_t)tmpPad->fClusterCandidates[candidate].fMean>AliHLTTPCTransform::GetNTimeBins()-100){
 189           continue;
 190         }
 191         tmpCandidate=&tmpPad->fClusterCandidates[candidate];
 192         tmpCandidate->fFirstPad=pad;
 193         UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
 194
 195         ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
 196         if(tmpCandidate->fTotalCharge>tmpTotalCharge){
 197           //we have a cluster
 198           tmpCandidate->fRowNumber= tmpCandidate->fRowNumber*tmpCandidate->fTotalCharge;
 199           //      tmpCandidate->fPad = tmpCandidate->fPad*tmpCandidate->fTotalCharge;
 200           fClusters.push_back(*tmpCandidate);
 201         }
 202       }
 203     }
 204   }
 205
 206   HLTDebug("Found %d normal clusters.",fClusters.size());
 207 }
 208
 209 void AliHLTTPCKryptonClusterFinder::FindKryptonClusters()
 210 {
 211   fMaxQOfCluster=0;
 212   if(fClusters.size()<2){
 213     return;
 214   }
 215   for(UInt_t i=0; i<fClusters.size();i++){
 216     AliHLTTPCClusters * tmpCluster=&fClusters[i];
 217     if(tmpCluster->fFlags==99){//quickfix to check if a cluster is used already
 218       continue;
 219     }
 220
 221     UInt_t prevRow=tmpCluster->fRowNumber/tmpCluster->fTotalCharge;
 222
 223     if((Int_t)tmpCluster->fQMax>fMaxQOfCluster){
 224       fMaxQOfCluster = tmpCluster->fQMax;
 225     }
 226
 227     //adds "normal" clusters belonging to the krypton cluster
 228     for(UInt_t j=i+1;j<fClusters.size();j++){
 229       AliHLTTPCClusters * nextCluster=&fClusters[j];
 230
 231       if(nextCluster->fFlags==99){//quickfix to check if a cluster is used already
 232         continue;
 233       }
 234       if(prevRow == (UInt_t)(nextCluster->fRowNumber/nextCluster->fTotalCharge)-1){//Checks if the row numbers are ok (next to eachother)
 235         if(abs((Int_t)(tmpCluster->fPad/tmpCluster->fTotalCharge) - (Int_t)(nextCluster->fPad/nextCluster->fTotalCharge))<3){ // checks if the pad numbers are ok
 236           if(abs((Int_t)tmpCluster->fMean-(Int_t)nextCluster->fMean)<2){
 237             tmpCluster->fMean=nextCluster->fMean;
 238             tmpCluster->fTotalCharge+=nextCluster->fTotalCharge;
 239             tmpCluster->fRowNumber+=nextCluster->fRowNumber;
 240             tmpCluster->fPad+=nextCluster->fPad;
 241             tmpCluster->fTime+=nextCluster->fTime;
 242             if((Int_t)nextCluster->fQMax>fMaxQOfCluster){
 243               fMaxQOfCluster = nextCluster->fQMax;
 244             }
 245
 246
 247
 248             if(tmpCluster->fFlags!=99){//means that this is the first time normal clusters match
 249               CheckForCandidateOnPreviousRow(tmpCluster);
 250               Int_t minFirst=0;
 251               Int_t maxFirst=0;
 252               if(tmpCluster->fFirstPad>1){
 253                 minFirst=2;
 254               }
 255               if(tmpCluster->fLastMergedPad+2<(UInt_t)AliHLTTPCTransform::GetNPads(prevRow)){
 256                 maxFirst=2;
 257               }
 258
 259               for(UInt_t ap = tmpCluster->fFirstPad -minFirst; ap<=tmpCluster->fLastMergedPad+maxFirst; ap++){
 260                 fHWAddressVector.push_back((AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(tmpCluster->fRowNumber/tmpCluster->fTotalCharge-AliHLTTPCTransform::GetFirstRow(fCurrentPatch),ap));
 261               }
 262             }
 263
 264             UInt_t minNext=0;
 265             UInt_t maxNext=0;
 266             if(nextCluster->fFirstPad>1){
 267               minNext=2;
 268             }
 269             if(nextCluster->fLastMergedPad+2<(UInt_t)AliHLTTPCTransform::GetNPads(prevRow+1)){
 270               maxNext=2;
 271             }
 272             for(UInt_t ap = nextCluster->fFirstPad-minNext; ap<=nextCluster->fLastMergedPad+maxNext; ap++){
 273               fHWAddressVector.push_back((AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(nextCluster->fRowNumber/nextCluster->fTotalCharge-AliHLTTPCTransform::GetFirstRow(fCurrentPatch),ap));
 274               HLTDebug("Pushing back hw address %d from row: %d and Pad: %d",fDigitReader->GetAltroBlockHWaddr(nextCluster->fRowNumber/nextCluster->fTotalCharge-AliHLTTPCTransform::GetFirstRow(fCurrentPatch),ap),nextCluster->fRowNumber/nextCluster->fTotalCharge,ap);
 275             }
 276
 277             prevRow=nextCluster->fRowNumber/nextCluster->fTotalCharge;
 278             nextCluster->fFlags=99;
 279             tmpCluster->fFlags=99;
 280             if(j!=fClusters.size()-1){
 281               continue;
 282             }
 283           }
 284         }
 285       }
 286
 287       if(tmpCluster->fFlags==99){
 288         //adds a clustercandidate on next row if present TODO
 289         /*      if(tmpCluster->fFlags==99){
 290           for(Int_t p=-1;p<2;p++){
 291           AliHLTTPCPad *padAfter=fRowPadVector[tmpCluster->fRowNumber+1][tmpCluster->fPad+p];
 292           for(UInt_t c=0;c<padAfter->fClusterCandidates.size();c++)
 293           if(abs((Int_t)tmpCluster->fMean - (Int_t)padAfter->fClusterCandidates[c].fMean)<2){
 294           tmpCluster->fTotalCharge+=padAfter->fClusterCandidates[c].fTotalCharge;
 295           }
 296           }//end add clustercandidate if present
 297           }*/
 298
 299
 300         //convert the (AliHLTClusters) cluster to spacepointdata and add it to the output array.
 301         if (fNKryptonClusters*sizeof(AliHLTTPCSpacePointData)>=fMaxOutputSize*sizeof(AliHLTUInt32_t)) {
 302           HLTWarning("Buffer too small too add more spacepoints: %d of %d byte(s) already used",fNKryptonClusters*sizeof(AliHLTTPCSpacePointData) , fMaxOutputSize*sizeof(AliHLTUInt32_t));
 303           return;
 304         }
 305         if(tmpCluster->fTotalCharge>10){
 306           Float_t xyz[3]={0,0,0};
 307           Int_t thisrow=-1;
 308           Int_t thissector=-1;
 309           AliHLTTPCTransform::Slice2Sector(fCurrentSlice, (Int_t)(tmpCluster->fRowNumber/tmpCluster->fTotalCharge), thissector, thisrow);
 310           AliHLTTPCTransform::Raw2Local(xyz, thissector, thisrow,(Float_t)(tmpCluster->fPad/tmpCluster->fTotalCharge),(Float_t)(tmpCluster->fTime/tmpCluster->fTotalCharge));
 311           fSpacePointData[fNKryptonClusters].fID = fCurrentSlice*10 +fCurrentPatch;
 312           fSpacePointData[fNKryptonClusters].fX = xyz[0];
 313           fSpacePointData[fNKryptonClusters].fY = xyz[1];
 314           fSpacePointData[fNKryptonClusters].fZ = xyz[2];
 315           fSpacePointData[fNKryptonClusters].fCharge = tmpCluster->fTotalCharge;
 316           fSpacePointData[fNKryptonClusters].fQMax = tmpCluster->fQMax;
 317           fSpacePointData[fNKryptonClusters].fPadRow = tmpCluster->fRowNumber/tmpCluster->fTotalCharge;
 318           HLTDebug("Krypton cluster found");
 319           HLTDebug("xyz=[%f,%f,%f]",fSpacePointData[fNKryptonClusters].fX,fSpacePointData[fNKryptonClusters].fY,fSpacePointData[fNKryptonClusters].fZ);
 320           HLTDebug("TotalCharge = %d    and   QMax = %d",fSpacePointData[fNKryptonClusters].fCharge,fSpacePointData[fNKryptonClusters].fQMax);
 321           fNKryptonClusters++;
 322           break;
 323         }
 324       }
 325     }
 326   }//end add "normal" clusters belonging to the krypton cluster
 327
 328   //resets the candidates for every pad and the fClusters(row clusters)
 329   for(UInt_t row=0;row<fNumberOfRows;row++){
 330     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
 331       fRowPadVector[row][pad]->fClusterCandidates.clear();
 332     }
 333   }
 334   fClusters.clear();
 335 }
 336
 337 void AliHLTTPCKryptonClusterFinder::CheckForCandidateOnPreviousRow(AliHLTTPCClusters *tmpCluster){
 338   if(tmpCluster->fRowNumber>1){
 339     for(Int_t p=-1;p<2;p++){
 340       if(tmpCluster->fPad+p>0 && tmpCluster->fPad+p<fNumberOfPadsInRow[tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1]){
 341         if(tmpCluster->fTotalCharge==0){
 342           HLTDebug("Charge of tmpCluster in AliHLTTPCKryptonClusterFinder::CheckForCandidateOnPreviousRow is 0");
 343           return;
 344         }
 345         if((Int_t)(tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch))<0){
 346           HLTDebug("AliHLTTPCKryptonClusterFinder::CheckForCandidateOnPreviousRow:    Rownumber is below 0");
 347           return;
 348         }
 349         if(tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch)>fNumberOfRows){
 350           HLTDebug("AliHLTTPCKryptonClusterFinder::CheckForCandidateOnPreviousRow:    Rownumber is too high");
 351           return;
 352         }
 353         AliHLTTPCPad *prevPad=fRowPadVector[tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch)][tmpCluster->fPad/tmpCluster->fTotalCharge+p];
 354         for(UInt_t i=0;i<prevPad->fClusterCandidates.size();i++){
 355           if(abs((Int_t)prevPad->fClusterCandidates[i].fMean - (Int_t)tmpCluster->fMean)<2 && prevPad->fUsedClusterCandidates[i]==0){
 356             tmpCluster->fTotalCharge += prevPad->fClusterCandidates[i].fTotalCharge;
 357             fHWAddressVector.push_back((AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch),tmpCluster->fPad/tmpCluster->fTotalCharge+p));
 358               //            fHWAddressVector.push_back((AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(prevPad->GetRowNumber(),prevPad->GetPadNumber()));
 359             HLTDebug("Pushing back hw address %d",fDigitReader->GetAltroBlockHWaddr(tmpCluster->fRowNumber/tmpCluster->fTotalCharge-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch),tmpCluster->fPad/tmpCluster->fTotalCharge+p));
 360           }
 361         }
 362       }
 363     }
 364   }
 365 }
 366
 367 void AliHLTTPCKryptonClusterFinder::SetSelection(Int_t minRow, Int_t maxRow){
 368   fSelectionMinRowNumber=minRow;
 369   fSelectionMaxRowNumber=maxRow;
 370 }
 371