HLT/hough/AliL3HoughEval.cxx

   1 //Author:        Anders Strand Vestbo
   2 //Last Modified: 28.6.01
   3
   4 #include <TClonesArray.h>
   5 #include <TH2.h>
   6 #include <TH1.h>
   7 #include <TFile.h>
   8 #include <AliRun.h>
   9 #include <TParticle.h>
  10 #include <TTree.h>
  11
  12 #include "AliTPCParam.h"
  13 #include "AliSimDigits.h"
  14 #include "AliL3TrackArray.h"
  15 #include "AliL3Transform.h"
  16 #include "AliL3HoughTransformer.h"
  17 #include "AliL3Defs.h"
  18 #include "AliL3HoughTrack.h"
  19 #include "AliL3HoughEval.h"
  20
  21 ClassImp(AliL3HoughEval)
  22
  23
  24 AliL3HoughEval::AliL3HoughEval()
  25 {
  26   //Default constructor
  27   fTransform = new AliL3Transform();
  28   fHoughTransformer = 0;
  29   fNumOfRowsToMiss = 1;
  30   fNumOfPadsToLook = 1;
  31 }
  32
  33
  34 AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
  35 {
  36   //Constructor
  37   fHoughTransformer = transformer;
  38   fTransform = new AliL3Transform();
  39   fNumOfRowsToMiss = 1;
  40   fNumOfPadsToLook = 1;
  41 }
  42
  43
  44 AliL3HoughEval::~AliL3HoughEval()
  45 {
  46   //Destructor
  47   if(fTransform)
  48     delete fTransform;
  49   if(fMcTrackTable)
  50     delete [] fMcTrackTable;
  51 }
  52
  53
  54 Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
  55 {
  56   //Look at rawdata along the road specified by the track candidates.
  57
  58   if(!fHoughTransformer)
  59     {
  60       printf("\nAliL3HoughEval: No transformer object\n");
  61       return kFALSE;
  62     }
  63
  64   Int_t patch = fHoughTransformer->fPatch;
  65   Int_t slice = fHoughTransformer->fSlice;
  66
  67
  68   Int_t sector,row;
  69   Int_t lut_index;
  70   Char_t **track_lut = fHoughTransformer->fTrackTable;
  71   Int_t nrow=0,npixs=0;
  72   Float_t xyz[3];
  73
  74   for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
  75     {
  76       Int_t prow = padrow - NRows[patch][0];
  77       if(!track->GetCrossingPoint(padrow,xyz))
  78         {
  79           printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
  80           continue;
  81         }
  82
  83       fTransform->Slice2Sector(slice,padrow,sector,row);
  84       fTransform->Local2Raw(xyz,sector,row);
  85       npixs=0;
  86       //Look at both sides of the crossing point:
  87       for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
  88         {
  89           if(p<0 || p>fTransform->GetNPads(padrow)) continue;
  90           lut_index = (prow<<8) + p;
  91           if(track_lut[eta_index][lut_index]>0) //There was a signal here
  92             npixs++;
  93         }
  94
  95       if(npixs > 0)
  96         {
  97           nrow++;
  98         }
  99     }
 100   if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
 101     {
 102       track->SetEtaIndex(eta_index);
 103       if(remove)
 104         RemoveTrackFromImage(track,eta_index);
 105       return kTRUE;
 106     }
 107   else
 108     return kFALSE;
 109 }
 110
 111 void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
 112 {
 113   //Evalaute the track candidates by looking along the trajectory.
 114   //If remove is on, the pixels along the track will be removed.
 115
 116
 117   if(!fHoughTransformer)
 118     {
 119       printf("AliL3HoughEval: No transformer object\n");
 120       return;
 121     }
 122
 123   Int_t patch = fHoughTransformer->fPatch;
 124   Int_t slice = fHoughTransformer->fSlice;
 125
 126
 127   Int_t sector,row;
 128   Int_t lut_index;
 129   Char_t **track_lut = fHoughTransformer->fTrackTable;
 130   Int_t nrow,npixs;
 131   Float_t xyz[3];
 132
 133   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 134     {
 135       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 136       if(!track) {printf("No track\n"); return;}
 137
 138       nrow = 0;
 139
 140       for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
 141         {
 142           Int_t prow = padrow - NRows[patch][0];
 143           if(!track->GetCrossingPoint(padrow,xyz))
 144             {
 145               printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
 146               continue;
 147             }
 148
 149           fTransform->Slice2Sector(slice,padrow,sector,row);
 150           fTransform->Local2Raw(xyz,sector,row);
 151           npixs=0;
 152           //Look at both sides of the crossing point:
 153           for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
 154             {
 155               if(p<0 || p>fTransform->GetNPads(padrow)) continue;
 156               lut_index = (prow<<8) + p;
 157               if(track_lut[eta_index][lut_index]>0) //There was a signal here
 158                 npixs++;
 159             }
 160
 161           if(npixs > 0)
 162             {
 163               nrow++;
 164             }
 165         }
 166       if(nrow < NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)
 167         tracks->Remove(i); //this was not a good enough track
 168       else if(remove)
 169         RemoveTrackFromImage(track,eta_index); //this was a good track, so remove it from the image.
 170     }
 171
 172   tracks->Compress();
 173 }
 174
 175 void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
 176 {
 177   //Remove the pixels along the track in the image.
 178
 179   Int_t patch = fHoughTransformer->fPatch;
 180   Int_t slice = fHoughTransformer->fSlice;
 181
 182   Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
 183
 184   Int_t lut_index;
 185   Char_t **track_lut = fHoughTransformer->fTrackTable;
 186   Int_t sector,row;
 187   Float_t xyz[3];
 188   for(Int_t padrow = NRows[patch][0]; padrow<=NRows[patch][1]; padrow++)
 189     {
 190       Int_t prow = padrow - NRows[patch][0];
 191       if(!track->GetCrossingPoint(padrow,xyz))
 192         {
 193           printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
 194           continue;
 195         }
 196
 197       fTransform->Slice2Sector(slice,padrow,sector,row);
 198       fTransform->Local2Raw(xyz,sector,row);
 199
 200       for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
 201         {
 202           if(p<0 || p>fTransform->GetNPads(padrow)) continue;
 203           lut_index = (prow<<8) + p;
 204           track_lut[eta_index][lut_index] = -1; //remove it!!
 205         }
 206     }
 207
 208
 209 }
 210
 211 void AliL3HoughEval::DisplaySlice(TH2F *hist)
 212 {
 213
 214   AliL3Digits *pixel;
 215
 216   for(Int_t padrow = 0; padrow < 174; padrow++)
 217     {
 218
 219       for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
 220         {
 221
 222           Int_t sector,row;
 223           Float_t xyz_pix[3];
 224           fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
 225           fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
 226
 227           if(hist)
 228             {
 229               hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
 230             }
 231
 232         }
 233     }
 234
 235 }
 236
 237 void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
 238 {
 239   //define the particles that produce good enough signals to be recognized in the transform
 240
 241   Int_t num_eta_segments = fHoughTransformer->fNumEtaSegments;
 242   Double_t eta_max = fHoughTransformer->fEtaMax;
 243   fMcTrackTable = new Int_t[num_eta_segments];
 244   for(Int_t i=0; i<num_eta_segments; i++)
 245     fMcTrackTable[i]=0;
 246   Double_t etaslice = eta_max/num_eta_segments;
 247
 248   Int_t patch = fHoughTransformer->fPatch;
 249   Int_t slice = fHoughTransformer->fSlice;
 250
 251   TFile *file = new TFile(rootfile);
 252   file->cd();
 253
 254   AliRun *gAlice = (AliRun*)file->Get("gAlice");
 255   gAlice->GetEvent(0);
 256
 257   TClonesArray *particles=gAlice->Particles();
 258   Int_t np = particles->GetEntriesFast();
 259
 260   Int_t row_to_miss = fNumOfRowsToMiss;
 261   AliTPCParam *param = (AliTPCParam*)file->Get("75x40_100x60");
 262   Int_t zero = param->GetZeroSup();
 263   TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
 264   AliSimDigits da, *digits=&da;
 265   TD->GetBranch("Segment")->SetAddress(&digits); //Return pointer to branch segment.
 266   Int_t *good=new Int_t[np];
 267   Int_t *count = new Int_t[np]; //np number of particles.
 268   Int_t good_number = NRows[patch][1]-NRows[patch][0]-row_to_miss;
 269   Int_t i;
 270   for (i=0; i<np; i++) count[i]=0;
 271   Int_t sectors_by_rows=(Int_t)TD->GetEntries();
 272   for (i=0; i<sectors_by_rows; i++)
 273     {
 274       if (!TD->GetEvent(i)) continue;
 275       Int_t sec,row;
 276       param->AdjustSectorRow(digits->GetID(),sec,row);
 277       Int_t ss,sr;
 278       fTransform->Sector2Slice(ss,sr,sec,row);
 279       if(ss!=slice) continue;
 280       if(sr < NRows[patch][0]) continue;
 281       if(sr > NRows[patch][1]) break;
 282       digits->First();
 283       while (digits->Next()) {
 284         Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
 285         Short_t dig = digits->GetDigit(it,ip);
 286         Int_t idx0=digits->GetTrackID(it,ip,0);
 287         Int_t idx1=digits->GetTrackID(it,ip,1);
 288         Int_t idx2=digits->GetTrackID(it,ip,2);
 289         if (idx0>=0 && dig>=zero) count[idx0]+=1;
 290         if (idx1>=0 && dig>=zero) count[idx1]+=1;
 291         if (idx2>=0 && dig>=zero) count[idx2]+=1;
 292       }
 293       for (Int_t j=0; j<np; j++)
 294         {
 295           if (count[j]>1) {//at least two digits at this padrow
 296             good[j]++;
 297           }
 298           count[j]=0;
 299         }
 300     }
 301   delete[] count;
 302
 303   Int_t good_one=0;
 304   //TObjArray *part = new TObjArray(0,0);
 305   for(i=0; i<np; i++)
 306    {
 307      TParticle *p = (TParticle*)particles->UncheckedAt(i);
 308      if(p->GetFirstMother()>0) continue; //secondary particle
 309      if(good[i] < good_number) continue; //too few padrows
 310      Double_t ptg=p->Pt();
 311      if(ptg<pet) continue;
 312      Int_t eta_index = (Int_t)(p->Eta()/etaslice);
 313      if(eta_index < 0 || eta_index > num_eta_segments)
 314        continue;
 315      //fMcTrackTable[eta_index]++;
 316      //part->AddLast(p);
 317      good_one++;
 318    }
 319
 320   printf("nparticles %d\n",good_one);
 321   file->Close();
 322   delete [] good;
 323   delete file;
 324   //return part;
 325 }
 326
 327
 328 void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
 329 {
 330
 331   Int_t slice = fHoughTransformer->fSlice;
 332
 333   TFile *file = new TFile(rootfile);
 334   file->cd();
 335
 336   AliRun *gAlice = (AliRun*)file->Get("gAlice");
 337   gAlice->GetEvent(0);
 338
 339   TClonesArray *particles=gAlice->Particles();
 340   Int_t n=particles->GetEntriesFast();
 341   Float_t torad=TMath::Pi()/180;
 342   Float_t phi_min = slice*20 - 10;
 343   Float_t phi_max = slice*20 + 10;
 344
 345   for (Int_t j=0; j<n; j++) {
 346     TParticle *p=(TParticle*)particles->UncheckedAt(j);
 347     if (p->GetFirstMother()>=0) continue;  //secondary particle
 348     if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
 349     Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
 350     Double_t phi_part = TMath::ATan2(pyg,pxg);
 351     if (phi_part < 0) phi_part += 2*TMath::Pi();
 352
 353     if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 354     if(ptg<0.100) continue;
 355
 356     AliL3HoughTrack *sel_track=0;
 357
 358     Double_t min_dist = 10000;
 359
 360     for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
 361       {
 362         AliL3HoughTrack *track = (AliL3HoughTrack*)merged_tracks->GetCheckedTrack(t);
 363         if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
 364         Float_t phi0[2] = {track->GetPhi0(),0};
 365         fTransform->Local2GlobalAngle(phi0,slice);
 366         Double_t dPt = ptg - track->GetPt();
 367         Double_t dPhi0 = phi_part - phi0[0];
 368         Double_t distance = pow(dPt,2) + pow(dPhi0,2);
 369         if(distance < min_dist)
 370           {
 371             min_dist = distance;
 372             sel_track = track;
 373           }
 374
 375       }
 376     if(sel_track)
 377       sel_track->SetBestMCid(j,min_dist);
 378
 379     Double_t dpt = fabs(ptg-sel_track->GetPt());
 380     Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
 381     //printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
 382   }
 383   file->Close();
 384   delete file;
 385
 386 }