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