HLT/TPCLib/tracking/AliHLTTPCHoughEval.cxx

   1 // @(#) $Id$
   2 // origin: hough/AliL3HoughEval.cxx,v 1.28 Thu Jun 17 10:36:14 2004 UTC by cvetan
   3
   4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   5 //*-- Copyright &copy ALICE HLT Group
   6
   7 #include "AliHLTStdIncludes.h"
   8
   9 #include <TH1.h>
  10 #include <TFile.h>
  11
  12 #include "AliHLTTPCLogging.h"
  13 #include "AliHLTTPCHoughEval.h"
  14 #include "AliHLTTPCMemHandler.h"
  15 #include "AliHLTTPCTrackArray.h"
  16 #include "AliHLTTPCHoughTransformer.h"
  17 #include "AliHLTTPCDigitData.h"
  18 #include "AliHLTTPCHoughTrack.h"
  19 #include "AliHLTTPCTransform.h"
  20 #include "AliHLTTPCHistogram.h"
  21 #include "AliHLTTPCHistogram1D.h"
  22
  23 #if __GNUC__ == 3
  24 using namespace std;
  25 #endif
  26
  27 /** /class AliHLTTPCHoughEval
  28 //<pre>
  29 //_____________________________________________________________
  30 // AliHLTTPCHoughEval
  31 //
  32 // Evaluation class for tracklets produced by the Hough transform.
  33 //
  34 </pre>
  35 */
  36
  37 ClassImp(AliHLTTPCHoughEval)
  38
  39 AliHLTTPCHoughEval::AliHLTTPCHoughEval()
  40 {
  41   //default ctor
  42   fRemoveFoundTracks = kFALSE;
  43   fNumOfPadsToLook = 1;
  44   fNumOfRowsToMiss = 1;
  45   fEtaHistos=0;
  46   fRowPointers = 0;
  47 }
  48
  49
  50 AliHLTTPCHoughEval::~AliHLTTPCHoughEval()
  51 {
  52   //dtor
  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 AliHLTTPCHoughEval::InitTransformer(AliHLTTPCHoughTransformer *transformer)
  63 {
  64   //Init hough transformer
  65   fHoughTransformer = transformer;
  66   fSlice = fHoughTransformer->GetSlice();
  67   fPatch = fHoughTransformer->GetPatch();
  68   fNrows = AliHLTTPCTransform::GetLastRow(fPatch) - AliHLTTPCTransform::GetFirstRow(fPatch) + 1;
  69   fNEtaSegments = fHoughTransformer->GetNEtaSegments();
  70   fEtaMin = fHoughTransformer->GetEtaMin();
  71   fEtaMax = fHoughTransformer->GetEtaMax();
  72   fZVertex = fHoughTransformer->GetZVertex();
  73   GenerateLUT();
  74 }
  75
  76 void AliHLTTPCHoughEval::GenerateLUT()
  77 {
  78   //Generate a Look-up table, to limit the access to raw data
  79
  80   if(!fRowPointers)
  81     fRowPointers = new AliHLTTPCDigitRowData*[fNrows];
  82
  83   AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fHoughTransformer->GetDataPointer();
  84   if(!tempPt)
  85     printf("\nAliHLTTPCHoughEval::GenerateLUT : Zero data pointer\n");
  86
  87   for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
  88     {
  89       Int_t prow = i - AliHLTTPCTransform::GetFirstRow(fPatch);
  90       fRowPointers[prow] = tempPt;
  91       AliHLTTPCMemHandler::UpdateRowPointer(tempPt);
  92     }
  93
  94 }
  95
  96 Bool_t AliHLTTPCHoughEval::LookInsideRoad(AliHLTTPCHoughTrack *track,Int_t &nrowscrossed,Int_t *rowrange,Bool_t remove)
  97 {
  98   //Look at rawdata along the road specified by the track candidates.
  99   //If track is good, return true, if not return false.
 100
 101   Int_t sector,row;
 102
 103   Int_t nrow=0,npixs=0;//,rows_crossed=0;
 104   Float_t xyz[3];
 105
 106   Int_t totalcharge=0;//total charge along the road
 107
 108   //for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
 109   for(Int_t padrow = rowrange[0]; padrow<=rowrange[1]; padrow++)
 110     {
 111       Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
 112       if(track->IsHelix())
 113         {
 114           if(!track->GetCrossingPoint(padrow,xyz))
 115             {
 116               continue;
 117             }
 118         }
 119       else
 120         {
 121           track->GetLineCrossingPoint(padrow,xyz);
 122           xyz[0] += AliHLTTPCTransform::Row2X(track->GetFirstRow());
 123           Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 124           xyz[2] = r*track->GetTgl();
 125         }
 126
 127       AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
 128       AliHLTTPCTransform::Local2Raw(xyz,sector,row);
 129
 130       npixs=0;
 131
 132       //Get the timebins for this pad
 133       AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
 134       if(!tempPt)
 135         {
 136           printf("AliHLTTPCHoughEval::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           AliHLTTPCDigitData *digPt = tempPt->fDigitData;
 144           for(UInt_t j=0; j<tempPt->fNDigit; j++)
 145             {
 146               Int_t pad = digPt[j].fPad;
 147               Int_t charge = digPt[j].fCharge;
 148               if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
 149               if(pad < p) continue;
 150               if(pad > p) break;
 151               UShort_t time = digPt[j].fTime;
 152               Double_t eta = AliHLTTPCTransform::GetEta(fSlice,padrow,pad,time);
 153               Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);
 154               if(pixelindex != track->GetEtaIndex()) continue;
 155               totalcharge += 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   nrowscrossed += nrow; //Update the number of rows crossed.
 171
 172   if(nrow >= rowrange[1]-rowrange[0]+1 - fNumOfRowsToMiss)//this was a good track
 173     {
 174       if(fRemoveFoundTracks)
 175         {
 176           Int_t dummy=0;
 177           LookInsideRoad(track,dummy,rowrange,kTRUE);
 178         }
 179       return kTRUE;
 180     }
 181   else
 182     return kFALSE;
 183 }
 184
 185 void AliHLTTPCHoughEval::FindEta(AliHLTTPCTrackArray *tracks)
 186 {
 187   //Find the corresponding eta slice hough space
 188   Int_t sector,row;
 189   Float_t xyz[3];
 190
 191   Int_t ntracks = tracks->GetNTracks();
 192   fEtaHistos = new AliHLTTPCHistogram1D*[ntracks];
 193
 194   Char_t hname[100];
 195   for(Int_t i=0; i<ntracks; i++)
 196     {
 197       sprintf(hname,"etahist_%d",i);
 198       fEtaHistos[i] = new AliHLTTPCHistogram1D(hname,hname,100,0,1);
 199     }
 200   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
 201
 202   for(Int_t ntr=0; ntr<ntracks; ntr++)
 203     {
 204       AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(ntr);
 205       if(!track) continue;
 206       for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
 207         {
 208           Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
 209
 210           if(!track->GetCrossingPoint(padrow,xyz))
 211             {
 212               printf("AliHLTTPCHoughEval::LookInsideRoad : Track does not cross line!!\n");
 213               continue;
 214             }
 215
 216           AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
 217           AliHLTTPCTransform::Local2Raw(xyz,sector,row);
 218
 219           //Get the timebins for this pad
 220           AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
 221           if(!tempPt)
 222             {
 223               printf("AliHLTTPCHoughEval::LookInsideRoad : Zero data pointer\n");
 224               continue;
 225             }
 226
 227           //Look at both sides of the pad:
 228           for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
 229             {
 230               AliHLTTPCDigitData *digPt = tempPt->fDigitData;
 231               for(UInt_t j=0; j<tempPt->fNDigit; j++)
 232                 {
 233                   UChar_t pad = digPt[j].fPad;
 234                   Int_t charge = digPt[j].fCharge;
 235                   if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
 236                   if(pad < p) continue;
 237                   if(pad > p) break;
 238                   UShort_t time = digPt[j].fTime;
 239                   Double_t eta = AliHLTTPCTransform::GetEta(fSlice,padrow,pad,time);
 240                   Int_t pixelindex = (Int_t)(eta/etaslice);
 241                   if(pixelindex > track->GetEtaIndex()+1) continue;
 242                   if(pixelindex < 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       AliHLTTPCHistogram1D *hist = fEtaHistos[i];
 252       Int_t maxbin = hist->GetMaximumBin();
 253       Double_t maxvalue = hist->GetBinContent(maxbin);
 254       AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(i);
 255       if(!track) continue;
 256       if(hist->GetBinContent(maxbin-1)<maxvalue && hist->GetBinContent(maxbin+1)<maxvalue)
 257         {
 258           track->SetWeight((Int_t)maxvalue,kTRUE);
 259           track->SetEta(hist->GetBinCenter(maxbin));
 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 AliHLTTPCHoughEval::DisplayEtaSlice(Int_t etaindex,AliHLTTPCHistogram *hist)
 276 {
 277   //Display the current raw data inside the (slice,patch)
 278
 279   if(!hist)
 280     {
 281       printf("AliHLTTPCHoughEval::DisplayEtaSlice : No input histogram!\n");
 282       return;
 283     }
 284
 285   for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
 286     {
 287       Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
 288
 289       AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
 290       if(!tempPt)
 291         {
 292           printf("AliHLTTPCHoughEval::DisplayEtaSlice : Zero data pointer\n");
 293           continue;
 294         }
 295
 296       AliHLTTPCDigitData *digPt = tempPt->fDigitData;
 297       if((Int_t)tempPt->fRow != padrow)
 298         {
 299           printf("\nAliHLTTPCHoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
 300           return;
 301         }
 302       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 303         {
 304           UChar_t pad = digPt[j].fPad;
 305           UChar_t charge = digPt[j].fCharge;
 306           UShort_t time = digPt[j].fTime;
 307           if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
 308           Float_t xyz[3];
 309           Int_t sector,row;
 310           AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
 311           AliHLTTPCTransform::Raw2Local(xyz,sector,row,pad,time);
 312           xyz[2] -= fZVertex;
 313           Double_t eta = AliHLTTPCTransform::GetEta(xyz);
 314           Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 315           if(pixelindex != etaindex) continue;
 316           hist->Fill(xyz[0],xyz[1],charge);
 317         }
 318     }
 319
 320 }
 321
 322 void AliHLTTPCHoughEval::CompareMC(AliHLTTPCTrackArray */*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(AliHLTTPCLog::kError,"AliHLTTPCHoughEval::CompareMC","Input file")
 345             <<"Error in file reading"<<ENDLOG;
 346           return;
 347         }
 348     }
 349   else
 350     {
 351       LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughEval::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       AliHLTTPCHoughTrack *tr = (AliHLTTPCHoughTrack*)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 }