HLT/hough/AliHLTHoughEval.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliHLTStandardIncludes.h"
   7
   8 #ifdef use_root
   9 #include <TH1.h>
  10 #include <TFile.h>
  11 #endif
  12
  13 #include "AliHLTLogging.h"
  14 #include "AliHLTHoughEval.h"
  15 #include "AliHLTMemHandler.h"
  16 #include "AliHLTTrackArray.h"
  17 #include "AliHLTHoughBaseTransformer.h"
  18 #include "AliHLTDigitData.h"
  19 #include "AliHLTHoughTrack.h"
  20 #include "AliHLTTransform.h"
  21 #include "AliHLTHistogram.h"
  22 #include "AliHLTHistogram1D.h"
  23
  24 #if __GNUC__ == 3
  25 using namespace std;
  26 #endif
  27
  28 /** /class AliHLTHoughEval
  29 //<pre>
  30 //_____________________________________________________________
  31 // AliHLTHoughEval
  32 //
  33 // Evaluation class for tracklets produced by the Hough transform.
  34 //
  35 </pre>
  36 */
  37
  38 ClassImp(AliHLTHoughEval)
  39
  40 AliHLTHoughEval::AliHLTHoughEval()
  41 {
  42   //default ctor
  43   fRemoveFoundTracks = kFALSE;
  44   fNumOfPadsToLook = 1;
  45   fNumOfRowsToMiss = 1;
  46   fEtaHistos=0;
  47   fRowPointers = 0;
  48 }
  49
  50
  51 AliHLTHoughEval::~AliHLTHoughEval()
  52 {
  53   //dtor
  54   fHoughTransformer = 0;
  55   if(fRowPointers)
  56     {
  57       for(Int_t i=0; i<fNrows; i++)
  58         fRowPointers[i] = 0;
  59       delete [] fRowPointers;
  60     }
  61 }
  62
  63 void AliHLTHoughEval::InitTransformer(AliHLTHoughBaseTransformer *transformer)
  64 {
  65   //Init hough transformer
  66   fHoughTransformer = transformer;
  67   fSlice = fHoughTransformer->GetSlice();
  68   fPatch = fHoughTransformer->GetPatch();
  69   fNrows = AliHLTTransform::GetLastRow(fPatch) - AliHLTTransform::GetFirstRow(fPatch) + 1;
  70   fNEtaSegments = fHoughTransformer->GetNEtaSegments();
  71   fEtaMin = fHoughTransformer->GetEtaMin();
  72   fEtaMax = fHoughTransformer->GetEtaMax();
  73   fZVertex = fHoughTransformer->GetZVertex();
  74   GenerateLUT();
  75 }
  76
  77 void AliHLTHoughEval::GenerateLUT()
  78 {
  79   //Generate a Look-up table, to limit the access to raw data
  80
  81   if(!fRowPointers)
  82     fRowPointers = new AliHLTDigitRowData*[fNrows];
  83
  84   AliHLTDigitRowData *tempPt = (AliHLTDigitRowData*)fHoughTransformer->GetDataPointer();
  85   if(!tempPt)
  86     printf("\nAliHLTHoughEval::GenerateLUT : Zero data pointer\n");
  87
  88   for(Int_t i=AliHLTTransform::GetFirstRow(fPatch); i<=AliHLTTransform::GetLastRow(fPatch); i++)
  89     {
  90       Int_t prow = i - AliHLTTransform::GetFirstRow(fPatch);
  91       fRowPointers[prow] = tempPt;
  92       AliHLTMemHandler::UpdateRowPointer(tempPt);
  93     }
  94
  95 }
  96
  97 Bool_t AliHLTHoughEval::LookInsideRoad(AliHLTHoughTrack *track,Int_t &nrowscrossed,Int_t *rowrange,Bool_t remove)
  98 {
  99   //Look at rawdata along the road specified by the track candidates.
 100   //If track is good, return true, if not return false.
 101
 102   Int_t sector,row;
 103
 104   Int_t nrow=0,npixs=0;//,rows_crossed=0;
 105   Float_t xyz[3];
 106
 107   Int_t totalcharge=0;//total charge along the road
 108
 109   //for(Int_t padrow = AliHLTTransform::GetFirstRow(fPatch); padrow <= AliHLTTransform::GetLastRow(fPatch); padrow++)
 110   for(Int_t padrow = rowrange[0]; padrow<=rowrange[1]; padrow++)
 111     {
 112       Int_t prow = padrow - AliHLTTransform::GetFirstRow(fPatch);
 113       if(track->IsHelix())
 114         {
 115           if(!track->GetCrossingPoint(padrow,xyz))
 116             {
 117               continue;
 118             }
 119         }
 120       else
 121         {
 122           track->GetLineCrossingPoint(padrow,xyz);
 123           xyz[0] += AliHLTTransform::Row2X(track->GetFirstRow());
 124           Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 125           xyz[2] = r*track->GetTgl();
 126         }
 127
 128       AliHLTTransform::Slice2Sector(fSlice,padrow,sector,row);
 129       AliHLTTransform::Local2Raw(xyz,sector,row);
 130
 131       npixs=0;
 132
 133       //Get the timebins for this pad
 134       AliHLTDigitRowData *tempPt = fRowPointers[prow];
 135       if(!tempPt)
 136         {
 137           printf("AliHLTHoughEval::LookInsideRoad : Zero data pointer\n");
 138           continue;
 139         }
 140
 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++)
 143         {
 144           AliHLTDigitData *digPt = tempPt->fDigitData;
 145           for(UInt_t j=0; j<tempPt->fNDigit; j++)
 146             {
 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;
 151               if(pad > p) break;
 152               UShort_t time = digPt[j].fTime;
 153               Double_t eta = AliHLTTransform::GetEta(fSlice,padrow,pad,time);
 154               Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);
 155               if(pixelindex != track->GetEtaIndex()) continue;
 156               totalcharge += digPt[j].fCharge;
 157               if(remove)
 158                 digPt[j].fCharge = 0; //Erease the track from image
 159               npixs++;
 160             }
 161         }
 162
 163       if(npixs > 1)//At least 2 digits on this padrow
 164         {
 165           nrow++;
 166         }
 167     }
 168   if(remove)
 169     return kTRUE;
 170
 171   nrowscrossed += nrow; //Update the number of rows crossed.
 172
 173   if(nrow >= rowrange[1]-rowrange[0]+1 - fNumOfRowsToMiss)//this was a good track
 174     {
 175       if(fRemoveFoundTracks)
 176         {
 177           Int_t dummy=0;
 178           LookInsideRoad(track,dummy,rowrange,kTRUE);
 179         }
 180       return kTRUE;
 181     }
 182   else
 183     return kFALSE;
 184 }
 185
 186 void AliHLTHoughEval::FindEta(AliHLTTrackArray *tracks)
 187 {
 188   //Find the corresponding eta slice hough space
 189   Int_t sector,row;
 190   Float_t xyz[3];
 191
 192   Int_t ntracks = tracks->GetNTracks();
 193   fEtaHistos = new AliHLTHistogram1D*[ntracks];
 194
 195   Char_t hname[100];
 196   for(Int_t i=0; i<ntracks; i++)
 197     {
 198       sprintf(hname,"etahist_%d",i);
 199       fEtaHistos[i] = new AliHLTHistogram1D(hname,hname,100,0,1);
 200     }
 201   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
 202
 203   for(Int_t ntr=0; ntr<ntracks; ntr++)
 204     {
 205       AliHLTHoughTrack *track = (AliHLTHoughTrack*)tracks->GetCheckedTrack(ntr);
 206       if(!track) continue;
 207       for(Int_t padrow = AliHLTTransform::GetFirstRow(fPatch); padrow <= AliHLTTransform::GetLastRow(fPatch); padrow++)
 208         {
 209           Int_t prow = padrow - AliHLTTransform::GetFirstRow(fPatch);
 210
 211           if(!track->GetCrossingPoint(padrow,xyz))
 212             {
 213               printf("AliHLTHoughEval::LookInsideRoad : Track does not cross line!!\n");
 214               continue;
 215             }
 216
 217           AliHLTTransform::Slice2Sector(fSlice,padrow,sector,row);
 218           AliHLTTransform::Local2Raw(xyz,sector,row);
 219
 220           //Get the timebins for this pad
 221           AliHLTDigitRowData *tempPt = fRowPointers[prow];
 222           if(!tempPt)
 223             {
 224               printf("AliHLTHoughEval::LookInsideRoad : Zero data pointer\n");
 225               continue;
 226             }
 227
 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++)
 230             {
 231               AliHLTDigitData *digPt = tempPt->fDigitData;
 232               for(UInt_t j=0; j<tempPt->fNDigit; j++)
 233                 {
 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;
 238                   if(pad > p) break;
 239                   UShort_t time = digPt[j].fTime;
 240                   Double_t eta = AliHLTTransform::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);
 245                 }
 246             }
 247         }
 248     }
 249
 250   for(Int_t i=0; i<ntracks; i++)
 251     {
 252       AliHLTHistogram1D *hist = fEtaHistos[i];
 253       Int_t maxbin = hist->GetMaximumBin();
 254       Double_t maxvalue = hist->GetBinContent(maxbin);
 255       AliHLTHoughTrack *track = (AliHLTHoughTrack*)tracks->GetCheckedTrack(i);
 256       if(!track) continue;
 257       if(hist->GetBinContent(maxbin-1)<maxvalue && hist->GetBinContent(maxbin+1)<maxvalue)
 258         {
 259           track->SetWeight((Int_t)maxvalue,kTRUE);
 260           track->SetEta(hist->GetBinCenter(maxbin));
 261           track->SetNHits(track->GetWeight());
 262         }
 263       else
 264         {
 265           track->SetWeight(0);
 266           tracks->Remove(i); //remove this track, because it was not a peak
 267         }
 268     }
 269   tracks->Compress();
 270
 271   //for(Int_t i=0; i<ntracks; i++)
 272   //delete fEtaHistos[i];
 273   //delete [] fEtaHistos;
 274 }
 275
 276 void AliHLTHoughEval::DisplayEtaSlice(Int_t etaindex,AliHLTHistogram *hist)
 277 {
 278   //Display the current raw data inside the (slice,patch)
 279
 280   if(!hist)
 281     {
 282       printf("AliHLTHoughEval::DisplayEtaSlice : No input histogram!\n");
 283       return;
 284     }
 285
 286   for(Int_t padrow = AliHLTTransform::GetFirstRow(fPatch); padrow <= AliHLTTransform::GetLastRow(fPatch); padrow++)
 287     {
 288       Int_t prow = padrow - AliHLTTransform::GetFirstRow(fPatch);
 289
 290       AliHLTDigitRowData *tempPt = fRowPointers[prow];
 291       if(!tempPt)
 292         {
 293           printf("AliHLTHoughEval::DisplayEtaSlice : Zero data pointer\n");
 294           continue;
 295         }
 296
 297       AliHLTDigitData *digPt = tempPt->fDigitData;
 298       if((Int_t)tempPt->fRow != padrow)
 299         {
 300           printf("\nAliHLTHoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
 301           return;
 302         }
 303       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 304         {
 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;
 309           Float_t xyz[3];
 310           Int_t sector,row;
 311           AliHLTTransform::Slice2Sector(fSlice,padrow,sector,row);
 312           AliHLTTransform::Raw2Local(xyz,sector,row,pad,time);
 313           xyz[2] -= fZVertex;
 314           Double_t eta = AliHLTTransform::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);
 318         }
 319     }
 320
 321 }
 322
 323 #ifdef use_root
 324 void AliHLTHoughEval::CompareMC(AliHLTTrackArray */*tracks*/,Char_t */*trackfile*/,Int_t /*threshold*/)
 325 {
 326   /*
 327   struct GoodTrack goodtracks[15000];
 328   Int_t nt=0;
 329   ifstream in(trackfile);
 330   if(in)
 331     {
 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)
 336         {
 337           nt++;
 338           if (nt==15000)
 339             {
 340               cerr<<"Too many good tracks"<<endl;
 341               break;
 342             }
 343         }
 344       if (!in.eof())
 345         {
 346           LOG(AliHLTLog::kError,"AliHLTHoughEval::CompareMC","Input file")
 347             <<"Error in file reading"<<ENDLOG;
 348           return;
 349         }
 350     }
 351   else
 352     {
 353       LOG(AliHLTLog::kError,"AliHLTHoughEval::CompareMC","Input")
 354         <<"No input trackfile "<<trackfile<<ENDLOG;
 355     }
 356
 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++)
 360     {
 361       particles[i]=0;
 362       ftracks[i]=0;
 363     }
 364
 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);
 371
 372   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
 373   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 374     {
 375       AliHLTHoughTrack *tr = (AliHLTHoughTrack*)tracks->GetCheckedTrack(i);
 376       if(!tr) continue;
 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());
 383     }
 384   for(Int_t i=0; i<nt; i++)
 385     {
 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);
 393     }
 394
 395   Double_t found=0;
 396   Double_t good =0;
 397   for(Int_t i=0; i<fNEtaSegments; i++)
 398     {
 399       //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
 400       found += ftracks[i];
 401       good += particles[i];
 402     }
 403   printf("And the total efficiency was: %f\n",found/good);
 404
 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");
 410   ptgood->Write();
 411   ptfound->Write();
 412   pteff->Write();
 413   etafound->Write();
 414   etagood->Write();
 415   etaeff->Write();
 416   file->Close();
 417
 418   delete [] particles;
 419   delete [] ftracks;
 420   */
 421 }
 422
 423 #endif