When Pt is bad defined (ex. no field), the multiple scattering effect is calculated...
[u/mrichter/AliRoot.git] / 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