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