3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliL3StandardIncludes.h"
13 #include "AliL3Logging.h"
14 #include "AliL3HoughEval.h"
15 #include "AliL3MemHandler.h"
16 #include "AliL3TrackArray.h"
17 #include "AliL3HoughBaseTransformer.h"
18 #include "AliL3DigitData.h"
19 #include "AliL3HoughTrack.h"
20 #include "AliL3Transform.h"
21 #include "AliL3Histogram.h"
22 #include "AliL3Histogram1D.h"
28 /** /class AliL3HoughEval
30 //_____________________________________________________________
33 // Evaluation class for tracklets produced by the Hough transform.
38 ClassImp(AliL3HoughEval)
40 AliL3HoughEval::AliL3HoughEval()
43 fRemoveFoundTracks = kFALSE;
51 AliL3HoughEval::~AliL3HoughEval()
54 fHoughTransformer = 0;
57 for(Int_t i=0; i<fNrows; i++)
59 delete [] fRowPointers;
63 void AliL3HoughEval::InitTransformer(AliL3HoughBaseTransformer *transformer)
65 //Init hough transformer
66 fHoughTransformer = transformer;
67 fSlice = fHoughTransformer->GetSlice();
68 fPatch = fHoughTransformer->GetPatch();
69 fNrows = AliL3Transform::GetLastRow(fPatch) - AliL3Transform::GetFirstRow(fPatch) + 1;
70 fNEtaSegments = fHoughTransformer->GetNEtaSegments();
71 fEtaMin = fHoughTransformer->GetEtaMin();
72 fEtaMax = fHoughTransformer->GetEtaMax();
73 fZVertex = fHoughTransformer->GetZVertex();
77 void AliL3HoughEval::GenerateLUT()
79 //Generate a Look-up table, to limit the access to raw data
82 fRowPointers = new AliL3DigitRowData*[fNrows];
84 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
86 printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
88 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
90 Int_t prow = i - AliL3Transform::GetFirstRow(fPatch);
91 fRowPointers[prow] = tempPt;
92 AliL3MemHandler::UpdateRowPointer(tempPt);
97 Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t &nrowscrossed,Int_t *rowrange,Bool_t remove)
99 //Look at rawdata along the road specified by the track candidates.
100 //If track is good, return true, if not return false.
104 Int_t nrow=0,npixs=0;//,rows_crossed=0;
107 Int_t totalcharge=0;//total charge along the road
109 //for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
110 for(Int_t padrow = rowrange[0]; padrow<=rowrange[1]; padrow++)
112 Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
115 if(!track->GetCrossingPoint(padrow,xyz))
122 track->GetLineCrossingPoint(padrow,xyz);
123 xyz[0] += AliL3Transform::Row2X(track->GetFirstRow());
124 Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
125 xyz[2] = r*track->GetTgl();
128 AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
129 AliL3Transform::Local2Raw(xyz,sector,row);
133 //Get the timebins for this pad
134 AliL3DigitRowData *tempPt = fRowPointers[prow];
137 printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
141 //Look at both sides of the pad:
142 for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
144 AliL3DigitData *digPt = tempPt->fDigitData;
145 for(UInt_t j=0; j<tempPt->fNDigit; j++)
147 Int_t pad = digPt[j].fPad;
148 Int_t charge = digPt[j].fCharge;
149 if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
150 if(pad < p) continue;
152 UShort_t time = digPt[j].fTime;
153 Double_t eta = AliL3Transform::GetEta(fSlice,padrow,pad,time);
154 Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);
155 if(pixelindex != track->GetEtaIndex()) continue;
156 totalcharge += digPt[j].fCharge;
158 digPt[j].fCharge = 0; //Erease the track from image
163 if(npixs > 1)//At least 2 digits on this padrow
171 nrowscrossed += nrow; //Update the number of rows crossed.
173 if(nrow >= rowrange[1]-rowrange[0]+1 - fNumOfRowsToMiss)//this was a good track
175 if(fRemoveFoundTracks)
178 LookInsideRoad(track,dummy,rowrange,kTRUE);
186 void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
188 //Find the corresponding eta slice hough space
192 Int_t ntracks = tracks->GetNTracks();
193 fEtaHistos = new AliL3Histogram1D*[ntracks];
196 for(Int_t i=0; i<ntracks; i++)
198 sprintf(hname,"etahist_%d",i);
199 fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
201 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
203 for(Int_t ntr=0; ntr<ntracks; ntr++)
205 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
207 for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
209 Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
211 if(!track->GetCrossingPoint(padrow,xyz))
213 printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
217 AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
218 AliL3Transform::Local2Raw(xyz,sector,row);
220 //Get the timebins for this pad
221 AliL3DigitRowData *tempPt = fRowPointers[prow];
224 printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
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++)
231 AliL3DigitData *digPt = tempPt->fDigitData;
232 for(UInt_t j=0; j<tempPt->fNDigit; j++)
234 UChar_t pad = digPt[j].fPad;
235 Int_t charge = digPt[j].fCharge;
236 if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
237 if(pad < p) continue;
239 UShort_t time = digPt[j].fTime;
240 Double_t eta = AliL3Transform::GetEta(fSlice,padrow,pad,time);
241 Int_t pixelindex = (Int_t)(eta/etaslice);
242 if(pixelindex > track->GetEtaIndex()+1) continue;
243 if(pixelindex < track->GetEtaIndex()-1) break;
244 fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
250 for(Int_t i=0; i<ntracks; i++)
252 AliL3Histogram1D *hist = fEtaHistos[i];
253 Int_t maxbin = hist->GetMaximumBin();
254 Double_t maxvalue = hist->GetBinContent(maxbin);
255 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
257 if(hist->GetBinContent(maxbin-1)<maxvalue && hist->GetBinContent(maxbin+1)<maxvalue)
259 track->SetWeight((Int_t)maxvalue,kTRUE);
260 track->SetEta(hist->GetBinCenter(maxbin));
261 track->SetNHits(track->GetWeight());
266 tracks->Remove(i); //remove this track, because it was not a peak
271 //for(Int_t i=0; i<ntracks; i++)
272 //delete fEtaHistos[i];
273 //delete [] fEtaHistos;
276 void AliL3HoughEval::DisplayEtaSlice(Int_t etaindex,AliL3Histogram *hist)
278 //Display the current raw data inside the (slice,patch)
282 printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
286 for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
288 Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
290 AliL3DigitRowData *tempPt = fRowPointers[prow];
293 printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
297 AliL3DigitData *digPt = tempPt->fDigitData;
298 if((Int_t)tempPt->fRow != padrow)
300 printf("\nAliL3HoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
303 for(UInt_t j=0; j<tempPt->fNDigit; j++)
305 UChar_t pad = digPt[j].fPad;
306 UChar_t charge = digPt[j].fCharge;
307 UShort_t time = digPt[j].fTime;
308 if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
311 AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
312 AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
314 Double_t eta = AliL3Transform::GetEta(xyz);
315 Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
316 if(pixelindex != etaindex) continue;
317 hist->Fill(xyz[0],xyz[1],charge);
324 void AliL3HoughEval::CompareMC(AliL3TrackArray */*tracks*/,Char_t */*trackfile*/,Int_t /*threshold*/)
327 struct GoodTrack goodtracks[15000];
329 ifstream in(trackfile);
332 printf("Reading good tracks from file %s\n",trackfile);
333 while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
334 goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
335 goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
340 cerr<<"Too many good tracks"<<endl;
346 LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
347 <<"Error in file reading"<<ENDLOG;
353 LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
354 <<"No input trackfile "<<trackfile<<ENDLOG;
357 Int_t *particles = new Int_t[fNEtaSegments];
358 Int_t *ftracks = new Int_t[fNEtaSegments];
359 for(Int_t i=0; i<fNEtaSegments; i++)
365 TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
366 TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
367 TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
368 TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
369 TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
370 TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
372 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
373 for(Int_t i=0; i<tracks->GetNTracks(); i++)
375 AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
377 if(tr->GetWeight()<threshold) continue;
378 Int_t trackindex = tr->GetEtaIndex();
379 if(trackindex <0 || trackindex >= fNEtaSegments) continue;
380 ftracks[trackindex]++;
381 ptfound->Fill(tr->GetPt());
382 etafound->Fill(tr->GetEta());
384 for(Int_t i=0; i<nt; i++)
386 if(goodtracks[i].nhits < 174) continue;
387 if(goodtracks[i].pt < 0.2) continue;
388 Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
389 if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
390 particles[particleindex]++;
391 ptgood->Fill(goodtracks[i].pt);
392 etagood->Fill(goodtracks[i].eta);
397 for(Int_t i=0; i<fNEtaSegments; i++)
399 //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
401 good += particles[i];
403 printf("And the total efficiency was: %f\n",found/good);
405 ptgood->Sumw2(); ptfound->Sumw2();
406 etagood->Sumw2(); etafound->Sumw2();
407 pteff->Divide(ptfound,ptgood,1,1,"b");
408 etaeff->Divide(etafound,etagood,1,1,"b");
409 TFile *file = TFile::Open("eff.root","RECREATE");