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
  31 #if __GNUC__ >= 3
  32 using namespace std;
  33 #endif
  34
  35 ClassImp(AliHLTTPCKryptonClusterFinder)
  36
  37 AliHLTTPCKryptonClusterFinder::AliHLTTPCKryptonClusterFinder()
  38   : AliHLTTPCClusterFinder(),
  39     fTimebinsInBunch(),
  40     fIndexOfBunchStart(),
  41     fHKryptonSpectrumFullPatch(NULL),
  42     fHKryptonSpectrumSelection(NULL),
  43     fHNumberOfKryptonClusters(NULL),
  44     fHNumberOfKryptonClustersSelection(NULL),
  45     fHMaxQofKryptonClusterLast1000(NULL),
  46     fHMaxQofKryptonClusterSelection(NULL),
  47     fStartBinKryptonSpectrum(0),
  48     fEndBinKryptonSpectrum(3000),
  49     fStartBinMaxQ(0),
  50     fEndBinMaxQ(1000),
  51     fStartBinNumberOfKryptonClusters(0),
  52     fEndBinNumberOfKryptonClusters(1000),
  53     fSelectionMinRowNumber(0),
  54     fSelectionMaxRowNumber(0),
  55     fMaxQOfCluster(0),
  56     fMaxQOfClusterBin(0),
  57     fNumberOfKryptonClusters(0),
  58     fNumberOfKryptonClustersBin(0),
  59     fHistogramsInitialized(kFALSE)
  60 {
  61   //constructor
  62 }
  63
  64 AliHLTTPCKryptonClusterFinder::~AliHLTTPCKryptonClusterFinder(){
  65   if(fHKryptonSpectrumFullPatch){
  66     delete fHKryptonSpectrumFullPatch;
  67     fHKryptonSpectrumFullPatch=NULL;
  68   }
  69   if(fHKryptonSpectrumSelection){
  70     delete fHKryptonSpectrumSelection;
  71     fHKryptonSpectrumSelection=NULL;
  72   }
  73   if(fHNumberOfKryptonClusters){
  74     delete fHNumberOfKryptonClusters;
  75     fHNumberOfKryptonClusters=NULL;
  76   }
  77   if(fHMaxQofKryptonClusterLast1000){
  78     delete fHMaxQofKryptonClusterLast1000;
  79     fHMaxQofKryptonClusterLast1000=NULL;
  80   }
  81   if(fHMaxQofKryptonClusterSelection){
  82     delete fHMaxQofKryptonClusterSelection;
  83     fHMaxQofKryptonClusterSelection=NULL;
  84   }
  85 }
  86
  87 void AliHLTTPCKryptonClusterFinder::ReBunch(const UInt_t *bunchData , Int_t bunchSize){
  88   fTimebinsInBunch.clear();
  89   fIndexOfBunchStart.clear();
  90   Bool_t newBunch=kTRUE;
  91   Int_t currentBunchNumber=-1;
  92   for(Int_t i=0;i<bunchSize;i++){
  93     if(newBunch){
  94       if(bunchData[i]>5){
  95         fIndexOfBunchStart.push_back(i);
  96         fTimebinsInBunch.push_back(1);
  97         currentBunchNumber++;
  98         newBunch=kFALSE;
  99       }
 100     }
 101     else{
 102       if(bunchData[i]>0){
 103         fTimebinsInBunch[currentBunchNumber]++;
 104       }
 105       else{
 106         newBunch=kTRUE;
 107       }
 108     }
 109   }
 110 }
 111
 112 void AliHLTTPCKryptonClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
 113 {
 114   //set input pointer
 115   fPtr = (UChar_t*)ptr;
 116   fSize = size;
 117
 118   if(!fVectorInitialized){
 119     InitializePadArray();
 120    }
 121
 122   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 123
 124   while(fDigitReader->NextChannel()){
 125     UInt_t row=fDigitReader->GetRow();
 126     UInt_t pad=fDigitReader->GetPad();
 127
 128     fRowPadVector[row][pad]->ClearCandidates();
 129
 130     while(fDigitReader->NextBunch()){
 131       if(fDigitReader->GetBunchSize()>1){
 132         const UInt_t *bunchData= fDigitReader->GetSignals();
 133
 134         ReBunch(bunchData,fDigitReader->GetBunchSize());
 135         Int_t rebunchCount=fIndexOfBunchStart.size();
 136
 137         for(Int_t r=0;r<rebunchCount;r++){
 138           UInt_t time = fDigitReader->GetTime()+fIndexOfBunchStart[r];
 139           AliHLTTPCClusters candidate;
 140           candidate.fTotalCharge=0;
 141           candidate.fChargeFalling=0;//PS This variable will store teh maximum charge of the candidate, used later. Keep this in mind
 142           if(fTimebinsInBunch[r]>2){
 143             for(Int_t i=0;i<fTimebinsInBunch[r];i++){
 144               candidate.fTotalCharge+=bunchData[i + fIndexOfBunchStart[r]];
 145               candidate.fTime += time*bunchData[i + fIndexOfBunchStart[r]];
 146               candidate.fTime2 += time*time*bunchData[i + fIndexOfBunchStart[r]];
 147               if(bunchData[i + fIndexOfBunchStart[r]]>candidate.fChargeFalling){
 148                 candidate.fChargeFalling = bunchData[i + fIndexOfBunchStart[r]];
 149               }
 150               time++;
 151             }
 152             if(candidate.fTotalCharge>0){
 153               candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 154               candidate.fPad=candidate.fTotalCharge*pad;
 155               candidate.fPad2=candidate.fPad*pad;
 156               candidate.fLastMergedPad=pad;
 157               candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 158             }
 159             if(candidate.fTotalCharge>10 && candidate.fMean<924){
 160               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 161             }
 162           }
 163         }
 164       }
 165     }
 166   }
 167 }
 168
 169 Bool_t AliHLTTPCKryptonClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
 170   //Checking if we have a match on the next pad
 171   for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
 172     AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
 173     if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
 174       cluster->fMean=candidate->fMean;
 175       cluster->fTotalCharge+=candidate->fTotalCharge;
 176       cluster->fTime += candidate->fTime;
 177       cluster->fTime2 += candidate->fTime2;
 178       cluster->fPad+=candidate->fPad;
 179       cluster->fPad2=candidate->fPad2;
 180       cluster->fLastMergedPad=candidate->fPad;
 181
 182       if(candidate->fChargeFalling>cluster->fChargeFalling){
 183         cluster->fChargeFalling = candidate->fChargeFalling;
 184       }
 185
 186       //setting the matched pad to used
 187       nextPad->fUsedClusterCandidates[candidateNumber]=1;
 188       nextPadToRead++;
 189       if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
 190         nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
 191         ComparePads(nextPad,cluster,nextPadToRead);
 192       }
 193       else{
 194         return kFALSE;
 195       }
 196     }
 197     else{
 198       return kFALSE;
 199     }
 200   }
 201   return kFALSE;
 202 }
 203
 204 void AliHLTTPCKryptonClusterFinder::FindRowClusters()
 205 {
 206   // see header file for function documentation
 207
 208   AliHLTTPCClusters* tmpCandidate=NULL;
 209   for(UInt_t row=5;row<fNumberOfRows-5;row++){
 210     fRowOfFirstCandidate=row;
 211     for(UInt_t pad=5;pad<fNumberOfPadsInRow[row]-1-5;pad++){
 212       AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
 213       for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
 214         if(tmpPad->fUsedClusterCandidates[candidate]){
 215           continue;
 216         }
 217         if((Int_t)tmpPad->fClusterCandidates[candidate].fMean<100 || (Int_t)tmpPad->fClusterCandidates[candidate].fMean>AliHLTTPCTransform::GetNTimeBins()-100){
 218           continue;
 219         }
 220         tmpCandidate=&tmpPad->fClusterCandidates[candidate];
 221         UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
 222
 223         ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
 224         if(tmpCandidate->fTotalCharge>tmpTotalCharge){
 225           //we have a cluster
 226           tmpCandidate->fPad=tmpCandidate->fPad/tmpCandidate->fTotalCharge;
 227           fClusters.push_back(*tmpCandidate);
 228         }
 229       }
 230     }
 231   }
 232
 233   HLTDebug("Found %d normal clusters.",fClusters.size());
 234 }
 235
 236 void AliHLTTPCKryptonClusterFinder::FindKryptonClusters()
 237 {
 238   if(fClusters.size()<2){
 239     return;
 240   }
 241   fMaxQOfCluster=0;
 242   for(UInt_t i=0; i<fClusters.size();i++){
 243     AliHLTTPCClusters * tmpCluster=&fClusters[i];
 244     if(tmpCluster->fFlags==99){//quickfix to check if a cluster is used already
 245       continue;
 246     }
 247     if((Int_t)tmpCluster->fChargeFalling>fMaxQOfCluster){
 248       fMaxQOfCluster = tmpCluster->fChargeFalling;
 249     }
 250
 251     //adds "normal" clusters belonging to the krypton cluster
 252     for(UInt_t j=i+1;j<fClusters.size();j++){
 253       AliHLTTPCClusters * nextCluster=&fClusters[j];
 254
 255       if(nextCluster->fFlags==99){//quickfix to check if a cluster is used already
 256         continue;
 257       }
 258
 259       if(tmpCluster->fRowNumber==nextCluster->fRowNumber-1){//Checks if the row numbers are ok (next to eachother)
 260         if(abs((Int_t)(tmpCluster->fPad) - (Int_t)(nextCluster->fPad))<3){ // checks if the pad numbers are ok
 261           if(abs((Int_t)tmpCluster->fMean-(Int_t)nextCluster->fMean)<2){
 262             tmpCluster->fMean=nextCluster->fMean;
 263             tmpCluster->fTotalCharge+=nextCluster->fTotalCharge;
 264             tmpCluster->fPad=nextCluster->fPad;
 265
 266             if((Int_t)nextCluster->fChargeFalling>fMaxQOfCluster){
 267               fMaxQOfCluster = nextCluster->fChargeFalling;
 268             }
 269
 270             if(tmpCluster->fFlags!=99){//means that this is the first time normal clusters match
 271               CheckForCandidateOnPreviousRow(tmpCluster);
 272             }
 273             tmpCluster->fRowNumber=nextCluster->fRowNumber;
 274             nextCluster->fFlags=99;
 275             tmpCluster->fFlags=99;
 276             if(j!=fClusters.size()-1){
 277               continue;
 278             }
 279
 280           }
 281         }
 282       }
 283
 284       if(tmpCluster->fFlags==99){
 285         //adds a clustercandidate on next row if present TODO
 286         /*      if(tmpCluster->fFlags==99){
 287           for(Int_t p=-1;p<2;p++){
 288           AliHLTTPCPad *padAfter=fRowPadVector[tmpCluster->fRowNumber+1][tmpCluster->fPad+p];
 289           for(UInt_t c=0;c<padAfter->fClusterCandidates.size();c++)
 290           if(abs((Int_t)tmpCluster->fMean - (Int_t)padAfter->fClusterCandidates[c].fMean)<2){
 291           tmpCluster->fTotalCharge+=padAfter->fClusterCandidates[c].fTotalCharge;
 292           }
 293           }//end add clustercandidate if present
 294           }*/
 295         fHMaxQofKryptonClusterLast1000->Fill(fMaxQOfClusterBin,fMaxQOfCluster);
 296         fMaxQOfClusterBin++;
 297         if(fMaxQOfClusterBin>fEndBinMaxQ){
 298           fMaxQOfClusterBin=0;
 299         }
 300         fHKryptonSpectrumFullPatch->Fill(tmpCluster->fTotalCharge);
 301         fHNumberOfKryptonClusters->Fill(fNumberOfKryptonClustersBin);
 302         if(fNumberOfKryptonClustersBin>fEndBinNumberOfKryptonClusters){
 303           fNumberOfKryptonClustersBin=0;
 304         }
 305         /*
 306         fHKryptonSpectrumSelection
 307         fHNumberOfKryptonClustersSelection
 308         fHMaxQofKryptonClusterSelection
 309         */
 310         HLTInfo("Krypton cluster found, charge: %d   in patch number: %d",tmpCluster->fTotalCharge,fCurrentPatch);
 311         break;
 312       }
 313
 314
 315     }
 316   }//end add "normal" clusters belonging to the krypton cluster
 317   fNumberOfKryptonClustersBin++;
 318
 319   //resets the candidates for every pad and the fClusters(row clusters)
 320   for(UInt_t row=5;row<fNumberOfRows-5;row++){
 321     for(UInt_t pad=5;pad<fNumberOfPadsInRow[row]-1-5;pad++){
 322       fRowPadVector[row][pad]->fClusterCandidates.clear();
 323     }
 324   }
 325   fClusters.clear();
 326 }
 327
 328 void AliHLTTPCKryptonClusterFinder::CheckForCandidateOnPreviousRow(AliHLTTPCClusters *tmpCluster){
 329   if(tmpCluster->fRowNumber>1){
 330     for(Int_t p=-1;p<2;p++){
 331       if(tmpCluster->fPad+p>0 && tmpCluster->fPad+p<fNumberOfPadsInRow[tmpCluster->fRowNumber-1]){
 332         AliHLTTPCPad *prevPad=fRowPadVector[tmpCluster->fRowNumber-1-AliHLTTPCTransform::GetFirstRow(fCurrentPatch)][tmpCluster->fPad+p];
 333         for(UInt_t i=0;i<prevPad->fClusterCandidates.size();i++){
 334           if(abs((Int_t)prevPad->fClusterCandidates[i].fMean - (Int_t)tmpCluster->fMean)<2 && prevPad->fUsedClusterCandidates[i]==0){
 335             tmpCluster->fTotalCharge += prevPad->fClusterCandidates[i].fTotalCharge;
 336           }
 337         }
 338       }
 339     }
 340   }
 341 }
 342
 343 void AliHLTTPCKryptonClusterFinder::InitializeHistograms(){
 344
 345
 346   if(fHistogramsInitialized){
 347     return;
 348   }
 349
 350   TString sliceStr("-Slice_");
 351   sliceStr+=fCurrentSlice;
 352   TString patchStr("-Patch_");
 353   patchStr+=fCurrentPatch;
 354
 355   TString namefKryptonSpectrumFullPatch("KryptonSpectrumFullPatch"+sliceStr+patchStr);
 356   TString namefKryptonSpectrumSelection = "KryptonSpectrumSelection"+sliceStr+patchStr;
 357   TString namefNumberOfKryptonClusters = "NumberOfKryptonClusters"+sliceStr+patchStr;
 358   TString namefNumberOfKryptonClustersSelection = "NumberOfKryptonClustersSelection"+sliceStr+patchStr;
 359   TString namefMaxQLast1000 = "MaxQ"+sliceStr+patchStr;
 360   TString namefMaxQSelection = "MaxQSelection"+sliceStr+patchStr;
 361
 362   fHKryptonSpectrumFullPatch = new TH1F(namefKryptonSpectrumFullPatch,namefKryptonSpectrumFullPatch,fEndBinKryptonSpectrum-fStartBinKryptonSpectrum,fStartBinKryptonSpectrum,fEndBinKryptonSpectrum);
 363
 364  fHKryptonSpectrumSelection = new TH1F(namefKryptonSpectrumSelection,namefKryptonSpectrumSelection,fEndBinKryptonSpectrum-fStartBinKryptonSpectrum,fStartBinKryptonSpectrum,fEndBinKryptonSpectrum);
 365
 366   fHNumberOfKryptonClusters = new TH1F(namefNumberOfKryptonClusters,namefNumberOfKryptonClusters,fEndBinNumberOfKryptonClusters-fStartBinNumberOfKryptonClusters,fStartBinNumberOfKryptonClusters,fEndBinNumberOfKryptonClusters);
 367
 368   fHNumberOfKryptonClustersSelection = new TH1F(namefNumberOfKryptonClustersSelection,namefNumberOfKryptonClustersSelection,fEndBinNumberOfKryptonClusters-fStartBinNumberOfKryptonClusters,fStartBinNumberOfKryptonClusters,fEndBinNumberOfKryptonClusters);
 369
 370   fHMaxQofKryptonClusterLast1000 = new TH1F(namefMaxQLast1000,namefMaxQLast1000,fEndBinMaxQ-fStartBinMaxQ,fStartBinMaxQ,fEndBinMaxQ);
 371
 372   fHMaxQofKryptonClusterSelection = new TH1F(namefMaxQSelection,namefMaxQSelection,fEndBinMaxQ-fStartBinMaxQ,fStartBinMaxQ,fEndBinMaxQ);
 373
 374   fHistogramsInitialized=kTRUE;
 375 }
 376
 377 void AliHLTTPCKryptonClusterFinder::ResetHistograms(TString histoName){
 378
 379   if (histoName.CompareTo("KryptonSpectrumFullPatch")==0){
 380     fHKryptonSpectrumFullPatch->Reset();
 381   }
 382   if (histoName.CompareTo("KryptonSpectrumSelection")==0){
 383     fHKryptonSpectrumSelection->Reset();
 384   }
 385   if (histoName.CompareTo("NumberOfKryptonClusters")==0){
 386     fHNumberOfKryptonClusters->Reset();
 387   }
 388   if (histoName.CompareTo("NumberOfKryptonClustersSelection")==0){
 389     fHNumberOfKryptonClustersSelection->Reset();
 390   }
 391   if (histoName.CompareTo("MaxQ")==0){
 392     fHMaxQofKryptonClusterLast1000->Reset();
 393   }
 394   if (histoName.CompareTo("MaxQSelection")==0){
 395     fHMaxQofKryptonClusterSelection->Reset();
 396   }
 397   if (histoName.CompareTo("All")==0){
 398     fHKryptonSpectrumFullPatch->Reset();
 399     fHKryptonSpectrumSelection->Reset();
 400     fHNumberOfKryptonClusters->Reset();
 401     fHNumberOfKryptonClustersSelection->Reset();
 402     fHMaxQofKryptonClusterLast1000->Reset();
 403     fHMaxQofKryptonClusterSelection->Reset();
 404   }
 405 }
 406
 407 void AliHLTTPCKryptonClusterFinder::SetSelection(Int_t minRow, Int_t maxRow){
 408   fSelectionMinRowNumber=minRow;
 409   fSelectionMaxRowNumber=maxRow;
 410 }
 411
 412 void AliHLTTPCKryptonClusterFinder::GetHistogramObjectArray(TObjArray & histos){
 413   //  TObjArray histos;
 414   histos.Clear();
 415   histos.Add(fHKryptonSpectrumFullPatch);
 416   histos.Add(fHKryptonSpectrumSelection);
 417   histos.Add(fHNumberOfKryptonClusters);
 418   histos.Add(fHNumberOfKryptonClustersSelection);
 419   histos.Add(fHMaxQofKryptonClusterLast1000);
 420   histos.Add(fHMaxQofKryptonClusterSelection);
 421   //  return histos;
 422 }
 423
 424 void AliHLTTPCKryptonClusterFinder::WriteHistograms(TString filename){
 425
 426   TFile file(filename,"recreate");
 427   fHKryptonSpectrumFullPatch->Write();
 428   fHKryptonSpectrumSelection->Write();
 429   fHNumberOfKryptonClusters->Write();
 430   fHNumberOfKryptonClustersSelection->Write();
 431   fHMaxQofKryptonClusterLast1000->Write();
 432   fHMaxQofKryptonClusterSelection->Write();
 433   file.Close();
 434 }