HLT/hough/AliL3HoughEval.cxx

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