HLT/hough/AliL3HoughEval.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <math.h>
   7 #include <iostream.h>
   8 #ifdef use_root
   9 #include <TH1.h>
  10 #include <TFile.h>
  11 #endif
  12 #include "AliL3MemHandler.h"
  13 #include "AliL3TrackArray.h"
  14 #include "AliL3Logging.h"
  15 #include "AliL3HoughEval.h"
  16 #include "AliL3HoughBaseTransformer.h"
  17 #include "AliL3DigitData.h"
  18 #include "AliL3HoughTrack.h"
  19 #include "AliL3Transform.h"
  20 #include "AliL3Histogram.h"
  21 #include "AliL3Histogram1D.h"
  22
  23 //_____________________________________________________________
  24 // AliL3HoughEval
  25 //
  26 // Evaluation class for tracklets produced by the Hough transform.
  27
  28 ClassImp(AliL3HoughEval)
  29
  30 AliL3HoughEval::AliL3HoughEval()
  31 {
  32
  33   fRemoveFoundTracks = kFALSE;
  34   fNumOfPadsToLook = 1;
  35   fNumOfRowsToMiss = 1;
  36   fEtaHistos=0;
  37   fRowPointers = 0;
  38 }
  39
  40
  41 AliL3HoughEval::~AliL3HoughEval()
  42 {
  43   fHoughTransformer = 0;
  44   if(fRowPointers)
  45     {
  46       for(Int_t i=0; i<fNrows; i++)
  47         fRowPointers[i] = 0;
  48       delete [] fRowPointers;
  49     }
  50 }
  51
  52 void AliL3HoughEval::InitTransformer(AliL3HoughBaseTransformer *transformer)
  53 {
  54   fHoughTransformer = transformer;
  55   fSlice = fHoughTransformer->GetSlice();
  56   fPatch = fHoughTransformer->GetPatch();
  57   fNrows = AliL3Transform::GetLastRow(fPatch) - AliL3Transform::GetFirstRow(fPatch) + 1;
  58   fNEtaSegments = fHoughTransformer->GetNEtaSegments();
  59   fEtaMin = fHoughTransformer->GetEtaMin();
  60   fEtaMax = fHoughTransformer->GetEtaMax();
  61   GenerateLUT();
  62 }
  63
  64 void AliL3HoughEval::GenerateLUT()
  65 {
  66   //Generate a Look-up table, to limit the access to raw data
  67
  68   if(!fRowPointers)
  69     fRowPointers = new AliL3DigitRowData*[fNrows];
  70
  71   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
  72   if(!tempPt)
  73     printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
  74
  75   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
  76     {
  77       Int_t prow = i - AliL3Transform::GetFirstRow(fPatch);
  78       fRowPointers[prow] = tempPt;
  79       AliL3MemHandler::UpdateRowPointer(tempPt);
  80     }
  81
  82 }
  83
  84 Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t &nrows_crossed,Bool_t remove)
  85 {
  86   //Look at rawdata along the road specified by the track candidates.
  87   //If track is good, return true, if not return false.
  88
  89   Int_t sector,row;
  90
  91   Int_t nrow=0,npixs=0;//,rows_crossed=0;
  92   Float_t xyz[3];
  93
  94   Int_t total_charge=0;//total charge along the road
  95
  96   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
  97     {
  98       //if(AliL3Transform::Row2X(padrow) > maxrow) break;//The track has left this slice
  99       //rows_crossed++;
 100
 101       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
 102       if(!track->GetCrossingPoint(padrow,xyz))
 103         {
 104           continue;
 105         }
 106
 107       AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
 108       AliL3Transform::Local2Raw(xyz,sector,row);
 109       npixs=0;
 110
 111       //Get the timebins for this pad
 112       AliL3DigitRowData *tempPt = fRowPointers[prow];
 113       if(!tempPt)
 114         {
 115           printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
 116           continue;
 117         }
 118
 119       //Look at both sides of the pad:
 120       for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
 121         {
 122           AliL3DigitData *digPt = tempPt->fDigitData;
 123           for(UInt_t j=0; j<tempPt->fNDigit; j++)
 124             {
 125               UChar_t pad = digPt[j].fPad;
 126               Int_t charge = digPt[j].fCharge;
 127               if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
 128               if(pad < p) continue;
 129               if(pad > p) break;
 130               UShort_t time = digPt[j].fTime;
 131               Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
 132               Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);
 133               if(pixel_index > track->GetEtaIndex()) continue;
 134               if(pixel_index != track->GetEtaIndex()) continue;
 135               total_charge += digPt[j].fCharge;
 136               if(remove)
 137                 digPt[j].fCharge = 0; //Erase the track from image
 138               npixs++;
 139             }
 140         }
 141
 142       if(npixs > 1)//At least 2 digits on this padrow
 143         {
 144           nrow++;
 145         }
 146     }
 147   if(remove)
 148     return kTRUE;
 149
 150   nrows_crossed += nrow; //Update the number of rows crossed.
 151
 152   if(nrow >= AliL3Transform::GetNRows(fPatch) - fNumOfRowsToMiss)//this was a good track
 153     {
 154       if(fRemoveFoundTracks)
 155         {
 156           Int_t dummy=0;
 157           LookInsideRoad(track,dummy,kTRUE);
 158         }
 159       return kTRUE;
 160     }
 161   else
 162     return kFALSE;
 163 }
 164
 165 void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
 166 {
 167
 168   Int_t sector,row;
 169   Float_t xyz[3];
 170
 171   Int_t ntracks = tracks->GetNTracks();
 172   fEtaHistos = new AliL3Histogram1D*[ntracks];
 173
 174   Char_t hname[100];
 175   for(Int_t i=0; i<ntracks; i++)
 176     {
 177       sprintf(hname,"etahist_%d",i);
 178       fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
 179     }
 180   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
 181
 182   for(Int_t ntr=0; ntr<ntracks; ntr++)
 183     {
 184       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
 185       if(!track) continue;
 186       for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
 187         {
 188           Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
 189
 190           if(!track->GetCrossingPoint(padrow,xyz))
 191             {
 192               printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
 193               continue;
 194             }
 195
 196           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
 197           AliL3Transform::Local2Raw(xyz,sector,row);
 198
 199           //Get the timebins for this pad
 200           AliL3DigitRowData *tempPt = fRowPointers[prow];
 201           if(!tempPt)
 202             {
 203               printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
 204               continue;
 205             }
 206
 207           //Look at both sides of the pad:
 208           for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
 209             {
 210               AliL3DigitData *digPt = tempPt->fDigitData;
 211               for(UInt_t j=0; j<tempPt->fNDigit; j++)
 212                 {
 213                   UChar_t pad = digPt[j].fPad;
 214                   Int_t charge = digPt[j].fCharge;
 215                   if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
 216                   if(pad < p) continue;
 217                   if(pad > p) break;
 218                   UShort_t time = digPt[j].fTime;
 219                   Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
 220                   Int_t pixel_index = (Int_t)(eta/etaslice);
 221                   if(pixel_index > track->GetEtaIndex()+1) continue;
 222                   if(pixel_index < track->GetEtaIndex()-1) break;
 223                   fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
 224                 }
 225             }
 226         }
 227     }
 228
 229   for(Int_t i=0; i<ntracks; i++)
 230     {
 231       AliL3Histogram1D *hist = fEtaHistos[i];
 232       Int_t max_bin = hist->GetMaximumBin();
 233       Double_t max_value = hist->GetBinContent(max_bin);
 234       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 235       if(!track) continue;
 236       if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
 237         {
 238           track->SetWeight((Int_t)max_value,kTRUE);
 239           track->SetEta(hist->GetBinCenter(max_bin));
 240           track->SetNHits(track->GetWeight());
 241         }
 242       else
 243         {
 244           track->SetWeight(0);
 245           tracks->Remove(i); //remove this track, because it was not a peak
 246         }
 247     }
 248   tracks->Compress();
 249
 250   //for(Int_t i=0; i<ntracks; i++)
 251   //delete fEtaHistos[i];
 252   //delete [] fEtaHistos;
 253 }
 254
 255 void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
 256 {
 257   //Display the current raw data inside the (slice,patch)
 258
 259   if(!hist)
 260     {
 261       printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
 262       return;
 263     }
 264
 265   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
 266     {
 267       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
 268
 269       AliL3DigitRowData *tempPt = fRowPointers[prow];
 270       if(!tempPt)
 271         {
 272           printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
 273           continue;
 274         }
 275
 276       AliL3DigitData *digPt = tempPt->fDigitData;
 277       if((Int_t)tempPt->fRow != padrow)
 278         {
 279           printf("\nAliL3HoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
 280           return;
 281         }
 282       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 283         {
 284           UChar_t pad = digPt[j].fPad;
 285           UChar_t charge = digPt[j].fCharge;
 286           UShort_t time = digPt[j].fTime;
 287           if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
 288           Float_t xyz[3];
 289           Int_t sector,row;
 290           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
 291           AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 292           Double_t eta = AliL3Transform::GetEta(xyz);
 293           Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 294           if(pixel_index != eta_index) continue;
 295           hist->Fill(xyz[0],xyz[1],charge);
 296         }
 297     }
 298
 299 }
 300
 301 #ifdef use_root
 302 void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile,Int_t threshold)
 303 {
 304   /*
 305   struct GoodTrack goodtracks[15000];
 306   Int_t nt=0;
 307   ifstream in(trackfile);
 308   if(in)
 309     {
 310       printf("Reading good tracks from file %s\n",trackfile);
 311       while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
 312              goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
 313              goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
 314         {
 315           nt++;
 316           if (nt==15000)
 317             {
 318               cerr<<"Too many good tracks"<<endl;
 319               break;
 320             }
 321         }
 322       if (!in.eof())
 323         {
 324           LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
 325             <<"Error in file reading"<<ENDLOG;
 326           return;
 327         }
 328     }
 329   else
 330     {
 331       LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
 332         <<"No input trackfile "<<trackfile<<ENDLOG;
 333     }
 334
 335   Int_t *particles = new Int_t[fNEtaSegments];
 336   Int_t *ftracks = new Int_t[fNEtaSegments];
 337   for(Int_t i=0; i<fNEtaSegments; i++)
 338     {
 339       particles[i]=0;
 340       ftracks[i]=0;
 341     }
 342
 343   TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
 344   TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
 345   TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
 346   TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
 347   TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
 348   TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
 349
 350   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
 351   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 352     {
 353       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 354       if(!tr) continue;
 355       if(tr->GetWeight()<threshold) continue;
 356       Int_t trackindex = tr->GetEtaIndex();
 357       if(trackindex <0 || trackindex >= fNEtaSegments) continue;
 358       ftracks[trackindex]++;
 359       ptfound->Fill(tr->GetPt());
 360       etafound->Fill(tr->GetEta());
 361     }
 362   for(Int_t i=0; i<nt; i++)
 363     {
 364       if(goodtracks[i].nhits < 174) continue;
 365       if(goodtracks[i].pt < 0.2) continue;
 366       Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
 367       if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
 368       particles[particleindex]++;
 369       ptgood->Fill(goodtracks[i].pt);
 370       etagood->Fill(goodtracks[i].eta);
 371     }
 372
 373   Double_t found=0;
 374   Double_t good =0;
 375   for(Int_t i=0; i<fNEtaSegments; i++)
 376     {
 377       //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
 378       found += ftracks[i];
 379       good += particles[i];
 380     }
 381   printf("And the total efficiency was: %f\n",found/good);
 382
 383   ptgood->Sumw2(); ptfound->Sumw2();
 384   etagood->Sumw2(); etafound->Sumw2();
 385   pteff->Divide(ptfound,ptgood,1,1,"b");
 386   etaeff->Divide(etafound,etagood,1,1,"b");
 387   TFile *file = TFile::Open("eff.root","RECREATE");
 388   ptgood->Write();
 389   ptfound->Write();
 390   pteff->Write();
 391   etafound->Write();
 392   etagood->Write();
 393   etaeff->Write();
 394   file->Close();
 395
 396   delete [] particles;
 397   delete [] ftracks;
 398   */
 399 }
 400
 401 #endif