UpdateRowPointer is now static function in memhandler.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx
CommitLineData
f80b98cb 1//Author: Anders Strand Vestbo
2//Last Modified: 28.6.01
3
162e2e5a 4#include <math.h>
5#include <TTree.h>
6#include <TFile.h>
7
8#include "GetGoodParticles.h"
9#include "AliL3TrackArray.h"
10#include "AliL3Logging.h"
99e7186b 11#include "AliL3HoughEval.h"
4de874d1 12#include "AliL3HoughTransformer.h"
99e7186b 13#include "AliL3DigitData.h"
4de874d1 14#include "AliL3HoughTrack.h"
99e7186b 15#include "AliL3Transform.h"
16#include "AliL3Histogram.h"
17#include "AliL3Defs.h"
4de874d1 18
19ClassImp(AliL3HoughEval)
20
4de874d1 21AliL3HoughEval::AliL3HoughEval()
22{
4fc9a6a4 23
99e7186b 24}
4de874d1 25
26AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
27{
4fc9a6a4 28
4de874d1 29 fHoughTransformer = transformer;
30 fTransform = new AliL3Transform();
99e7186b 31
32 fSlice = fHoughTransformer->GetSlice();
33 fPatch = fHoughTransformer->GetPatch();
34 fNrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
35 fNEtaSegments = fHoughTransformer->GetNEtaSegments();
36 fEtaMin = fHoughTransformer->GetEtaMin();
37 fEtaMax = fHoughTransformer->GetEtaMax();
38 fRemoveFoundTracks = kFALSE;
f80b98cb 39 fNumOfPadsToLook = 1;
99e7186b 40 fNumOfRowsToMiss = 1;
41 GenerateLUT();
4de874d1 42}
43
4de874d1 44AliL3HoughEval::~AliL3HoughEval()
45{
4de874d1 46 if(fTransform)
47 delete fTransform;
99e7186b 48 if(fRowPointers)
49 delete [] fRowPointers;
4de874d1 50}
51
99e7186b 52void AliL3HoughEval::GenerateLUT()
4de874d1 53{
99e7186b 54 //Generate a LUT, to limit the access to raw data
55
56 fRowPointers = new AliL3DigitRowData*[fNrows];
57
58 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
59 if(!tempPt)
60 printf("AliL3HoughEval::GenerateLUT : Zero data pointer\n");
f80b98cb 61
99e7186b 62 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
f80b98cb 63 {
99e7186b 64 Int_t prow = i - NRows[fPatch][0];
65 fRowPointers[prow] = tempPt;
4fc9a6a4 66 fHoughTransformer->UpdateDataPointer(tempPt);
f80b98cb 67 }
68
99e7186b 69}
f80b98cb 70
99e7186b 71Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
72{
73 //Look at rawdata along the road specified by the track candidates.
74 //If track is good, return true, if not return false.
162e2e5a 75
f80b98cb 76 Int_t sector,row;
99e7186b 77
f80b98cb 78 Int_t nrow=0,npixs=0;
79 Float_t xyz[3];
99e7186b 80
81 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
82 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
f80b98cb 83 {
99e7186b 84 Int_t prow = padrow - NRows[fPatch][0];
85
f80b98cb 86 if(!track->GetCrossingPoint(padrow,xyz))
87 {
88 printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
89 continue;
90 }
91
99e7186b 92 fTransform->Slice2Sector(fSlice,padrow,sector,row);
f80b98cb 93 fTransform->Local2Raw(xyz,sector,row);
94 npixs=0;
99e7186b 95
96 //Get the timebins for this pad
97 AliL3DigitRowData *tempPt = fRowPointers[prow];
98 if(!tempPt)
f80b98cb 99 {
99e7186b 100 printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
101 continue;
f80b98cb 102 }
103
99e7186b 104 //Look at both sides of the pad:
162e2e5a 105 for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
99e7186b 106 {
107 AliL3DigitData *digPt = tempPt->fDigitData;
108 for(UInt_t j=0; j<tempPt->fNDigit; j++)
109 {
110 UChar_t pad = digPt[j].fPad;
162e2e5a 111
99e7186b 112 if(pad < p) continue;
113 if(pad > p) break;
114 UShort_t time = digPt[j].fTime;
162e2e5a 115 Double_t eta = fTransform->GetEta(padrow,pad,time);
99e7186b 116 Int_t pixel_index = (Int_t)(eta/etaslice);
117 if(pixel_index > eta_index) continue;
118 if(pixel_index != eta_index) break;
119 if(remove)
120 digPt[j].fCharge = 0; //Delete the track from image
121 npixs++;
122 }
123 }
124
162e2e5a 125 if(npixs > 1)//At least 2 digits on this padrow
f80b98cb 126 {
127 nrow++;
128 }
129 }
99e7186b 130 if(remove)
131 return kTRUE;
132
133 if(nrow >= fNrows - fNumOfRowsToMiss)//this was a good track
f80b98cb 134 {
135 track->SetEtaIndex(eta_index);
99e7186b 136 if(fRemoveFoundTracks)
137 LookInsideRoad(track,eta_index,kTRUE);
f80b98cb 138 return kTRUE;
139 }
140 else
141 return kFALSE;
142}
143
99e7186b 144void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
f80b98cb 145{
99e7186b 146 //Display the current raw data inside the slice
4de874d1 147
99e7186b 148 if(!hist)
4de874d1 149 {
99e7186b 150 printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
4de874d1 151 return;
152 }
153
99e7186b 154 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
155 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
f80b98cb 156 {
99e7186b 157 Int_t prow = padrow - NRows[fPatch][0];
158
159 AliL3DigitRowData *tempPt = fRowPointers[prow];
160 if(!tempPt)
f80b98cb 161 {
99e7186b 162 printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
f80b98cb 163 continue;
164 }
165
99e7186b 166 AliL3DigitData *digPt = tempPt->fDigitData;
167 for(UInt_t j=0; j<tempPt->fNDigit; j++)
f80b98cb 168 {
99e7186b 169 UChar_t pad = digPt[j].fPad;
170 UChar_t charge = digPt[j].fCharge;
171 UShort_t time = digPt[j].fTime;
4fc9a6a4 172 if(charge < fHoughTransformer->GetThreshold()) continue;
99e7186b 173 Float_t xyz[3];
4de874d1 174 Int_t sector,row;
99e7186b 175 fTransform->Slice2Sector(fSlice,padrow,sector,row);
176 fTransform->Raw2Local(xyz,sector,row,pad,time);
177 Double_t eta = fTransform->GetEta(xyz);
178 Int_t pixel_index = (Int_t)(eta/etaslice);
179 if(pixel_index != eta_index) continue;
180 hist->Fill(xyz[0],xyz[1],charge);
f80b98cb 181 }
182 }
f80b98cb 183
4de874d1 184}
162e2e5a 185
186void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile)
187{
188
189 struct GoodTrack goodtracks[15000];
190 Int_t nt=0;
191 ifstream in(trackfile);
192 if(in)
193 {
194 printf("Reading good tracks from file %s\n",trackfile);
195 while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
196 goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
197 goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
198 {
199 nt++;
200 if (nt==15000)
201 {
202 cerr<<"Too many good tracks"<<endl;
203 break;
204 }
205 }
206 if (!in.eof())
207 {
208 LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
209 <<"Error in file reading"<<ENDLOG;
210 return;
211 }
212 }
213 else
214 {
215 LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
216 <<"No input trackfile "<<trackfile<<ENDLOG;
217 }
218
219 Int_t *particles = new Int_t[fNEtaSegments];
220 Int_t *ftracks = new Int_t[fNEtaSegments];
221 for(Int_t i=0; i<fNEtaSegments; i++)
222 {
223 particles[i]=0;
224 ftracks[i]=0;
225 }
226
227 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
228 for(Int_t i=0; i<tracks->GetNTracks(); i++)
229 {
230 AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
231 if(!tr) continue;
232 Int_t trackindex = tr->GetEtaIndex();
233 if(trackindex <0 || trackindex >= fNEtaSegments) continue;
234 ftracks[trackindex]++;
235 }
236 for(Int_t i=0; i<nt; i++)
237 {
238 if(goodtracks[i].nhits < 150) continue;
239 if(goodtracks[i].pt < 0.5) continue;
240 Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
241 if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
242 particles[particleindex]++;
243 }
244
245 Double_t found=0;
246 Double_t good =0;
247 for(Int_t i=0; i<fNEtaSegments; i++)
248 {
249 printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
250 found += ftracks[i];
251 good += particles[i];
252 }
253 printf("And the total efficiency was: %f\n",found/good);
254
255 delete [] particles;
256 delete [] ftracks;
257
258}
259