Some additional changes related to the previous changes. AliL3Transform
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx
CommitLineData
95a00d93 1//$Id$
2
b1886074 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
f80b98cb 5
162e2e5a 6#include <math.h>
95a00d93 7#include <iostream.h>
8#ifdef use_root
b46b21a5 9#include <TH1.h>
162e2e5a 10#include <TFile.h>
95a00d93 11#endif
b46b21a5 12#include "AliL3MemHandler.h"
162e2e5a 13#include "GetGoodParticles.h"
14#include "AliL3TrackArray.h"
15#include "AliL3Logging.h"
99e7186b 16#include "AliL3HoughEval.h"
237d3f5c 17#include "AliL3HoughBaseTransformer.h"
99e7186b 18#include "AliL3DigitData.h"
4de874d1 19#include "AliL3HoughTrack.h"
99e7186b 20#include "AliL3Transform.h"
21#include "AliL3Histogram.h"
b46b21a5 22#include "AliL3Histogram1D.h"
99e7186b 23#include "AliL3Defs.h"
4de874d1 24
b1886074 25//_____________________________________________________________
26// AliL3HoughEval
27//
28// Evaluation class for tracklets produced by the Hough transform.
29
4de874d1 30ClassImp(AliL3HoughEval)
31
4de874d1 32AliL3HoughEval::AliL3HoughEval()
4de874d1 33{
4fc9a6a4 34
4de874d1 35 fTransform = new AliL3Transform();
99e7186b 36 fRemoveFoundTracks = kFALSE;
f80b98cb 37 fNumOfPadsToLook = 1;
99e7186b 38 fNumOfRowsToMiss = 1;
b46b21a5 39 fEtaHistos=0;
b1886074 40 fRowPointers = 0;
4de874d1 41}
42
b46b21a5 43
4de874d1 44AliL3HoughEval::~AliL3HoughEval()
45{
b46b21a5 46 fHoughTransformer = 0;
4de874d1 47 if(fTransform)
48 delete fTransform;
99e7186b 49 if(fRowPointers)
b1886074 50 {
51 for(Int_t i=0; i<fNrows; i++)
52 fRowPointers[i] = 0;
53 delete [] fRowPointers;
54 }
4de874d1 55}
56
237d3f5c 57void AliL3HoughEval::InitTransformer(AliL3HoughBaseTransformer *transformer)
b46b21a5 58{
59 fHoughTransformer = transformer;
60 fSlice = fHoughTransformer->GetSlice();
61 fPatch = fHoughTransformer->GetPatch();
62 fNrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
63 fNEtaSegments = fHoughTransformer->GetNEtaSegments();
64 fEtaMin = fHoughTransformer->GetEtaMin();
65 fEtaMax = fHoughTransformer->GetEtaMax();
66 GenerateLUT();
67}
68
99e7186b 69void AliL3HoughEval::GenerateLUT()
4de874d1 70{
b46b21a5 71 //Generate a Look-up table, to limit the access to raw data
99e7186b 72
b1886074 73 if(!fRowPointers)
74 fRowPointers = new AliL3DigitRowData*[fNrows];
99e7186b 75
76 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
77 if(!tempPt)
b46b21a5 78 printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
f80b98cb 79
99e7186b 80 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
f80b98cb 81 {
99e7186b 82 Int_t prow = i - NRows[fPatch][0];
83 fRowPointers[prow] = tempPt;
b46b21a5 84 AliL3MemHandler::UpdateRowPointer(tempPt);
f80b98cb 85 }
86
99e7186b 87}
f80b98cb 88
99e7186b 89Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
90{
91 //Look at rawdata along the road specified by the track candidates.
92 //If track is good, return true, if not return false.
162e2e5a 93
f80b98cb 94 Int_t sector,row;
99e7186b 95
5de3a7f2 96 Int_t nrow=0,npixs=0,rows_crossed=0;
f80b98cb 97 Float_t xyz[3];
99e7186b 98
b46b21a5 99 Int_t total_charge=0;//total charge along the road
99e7186b 100 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
95a00d93 101
102
103 //Check if the track is leaving the sector at some point
5de3a7f2 104 Float_t maxrow=300;
105 Double_t angle=Pi/18;
106 track->CalculateEdgePoint(angle);
107 if(!track->IsPoint())
108 {
109 track->CalculateEdgePoint(-1.*angle);
110 if(track->IsPoint())
111 maxrow = track->GetPointX();
112 }
113 else
114 maxrow = track->GetPointX();
115
99e7186b 116 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
f80b98cb 117 {
5de3a7f2 118 if(fTransform->Row2X(padrow) > maxrow) break;//The track has left this slice
119 rows_crossed++;
99e7186b 120 Int_t prow = padrow - NRows[fPatch][0];
f80b98cb 121 if(!track->GetCrossingPoint(padrow,xyz))
122 {
5de3a7f2 123 //printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!; pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
f80b98cb 124 continue;
125 }
126
99e7186b 127 fTransform->Slice2Sector(fSlice,padrow,sector,row);
f80b98cb 128 fTransform->Local2Raw(xyz,sector,row);
129 npixs=0;
99e7186b 130
95a00d93 131
99e7186b 132 //Get the timebins for this pad
133 AliL3DigitRowData *tempPt = fRowPointers[prow];
134 if(!tempPt)
f80b98cb 135 {
99e7186b 136 printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
137 continue;
f80b98cb 138 }
139
99e7186b 140 //Look at both sides of the pad:
162e2e5a 141 for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
99e7186b 142 {
143 AliL3DigitData *digPt = tempPt->fDigitData;
144 for(UInt_t j=0; j<tempPt->fNDigit; j++)
145 {
b1886074 146 //if(digPt->fCharge <= fHoughTransformer->GetThreshold()) continue;
99e7186b 147 UChar_t pad = digPt[j].fPad;
148 if(pad < p) continue;
149 if(pad > p) break;
150 UShort_t time = digPt[j].fTime;
162e2e5a 151 Double_t eta = fTransform->GetEta(padrow,pad,time);
99e7186b 152 Int_t pixel_index = (Int_t)(eta/etaslice);
153 if(pixel_index > eta_index) continue;
154 if(pixel_index != eta_index) break;
b46b21a5 155 total_charge += digPt[j].fCharge;
99e7186b 156 if(remove)
b46b21a5 157 digPt[j].fCharge = 0; //Erease the track from image
99e7186b 158 npixs++;
159 }
160 }
161
162e2e5a 162 if(npixs > 1)//At least 2 digits on this padrow
f80b98cb 163 {
164 nrow++;
165 }
166 }
99e7186b 167 if(remove)
168 return kTRUE;
169
5de3a7f2 170 if(nrow >= rows_crossed - fNumOfRowsToMiss)//this was a good track
f80b98cb 171 {
b46b21a5 172 Double_t eta_track = (Double_t)eta_index*etaslice;
f80b98cb 173 track->SetEtaIndex(eta_index);
b46b21a5 174 track->SetWeight(total_charge,kTRUE);
175 track->SetEta(eta_track);
b1886074 176 track->SetRowRange(NRows[fPatch][0],NRows[fPatch][1]);
177 track->SetSlice(fSlice);
99e7186b 178 if(fRemoveFoundTracks)
179 LookInsideRoad(track,eta_index,kTRUE);
f80b98cb 180 return kTRUE;
181 }
182 else
183 return kFALSE;
184}
185
b46b21a5 186void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
187{
188
189 Int_t sector,row;
190 Float_t xyz[3];
191
192 Int_t ntracks = tracks->GetNTracks();
193 fEtaHistos = new AliL3Histogram1D*[ntracks];
194
195 Char_t hname[100];
196 for(Int_t i=0; i<ntracks; i++)
197 {
198 sprintf(hname,"etahist_%d",i);
199 fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
200 }
201 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
202
203 for(Int_t ntr=0; ntr<ntracks; ntr++)
204 {
205 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
206 if(!track) continue;
207 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
208 {
209 Int_t prow = padrow - NRows[fPatch][0];
210
211 if(!track->GetCrossingPoint(padrow,xyz))
212 {
213 printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
214 continue;
215 }
216
217 fTransform->Slice2Sector(fSlice,padrow,sector,row);
218 fTransform->Local2Raw(xyz,sector,row);
219
220 //Get the timebins for this pad
221 AliL3DigitRowData *tempPt = fRowPointers[prow];
222 if(!tempPt)
223 {
224 printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
225 continue;
226 }
227
228 //Look at both sides of the pad:
229 for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
230 {
231 AliL3DigitData *digPt = tempPt->fDigitData;
232 for(UInt_t j=0; j<tempPt->fNDigit; j++)
233 {
234 UChar_t pad = digPt[j].fPad;
235
236 if(pad < p) continue;
237 if(pad > p) break;
238 UShort_t time = digPt[j].fTime;
239 Double_t eta = fTransform->GetEta(padrow,pad,time);
240 Int_t pixel_index = (Int_t)(eta/etaslice);
241 if(pixel_index > track->GetEtaIndex()+1) continue;
242 if(pixel_index < track->GetEtaIndex()-1) break;
243 fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
244 }
245 }
246 }
247 }
248
249 for(Int_t i=0; i<ntracks; i++)
250 {
251 AliL3Histogram1D *hist = fEtaHistos[i];
252 Int_t max_bin = hist->GetMaximumBin();
253 Double_t max_value = hist->GetBinContent(max_bin);
254 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
255 if(!track) continue;
256 if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
257 {
258 track->SetWeight((Int_t)max_value,kTRUE);
259 track->SetEta(hist->GetBinCenter(max_bin));
260 track->SetNHits(track->GetWeight());
261 }
262 else
263 {
264 track->SetWeight(0);
265 tracks->Remove(i); //remove this track, because it was not a peak
266 }
267 }
268 tracks->Compress();
269
270 //for(Int_t i=0; i<ntracks; i++)
271 //delete fEtaHistos[i];
272