HLT/hough/AliL3HoughEval.cxx

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