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