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