HLT/hough/AliL3HoughEval.cxx

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