HLT/hough/AliL3HoughTransformer.cxx

   1 #include <math.h>
   2
   3 #include <TFile.h>
   4 #include <TTree.h>
   5 #include <TH2.h>
   6
   7 #include "AliTPCParam.h"
   8 #include "AliSimDigits.h"
   9
  10 #include "AliL3Defs.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3HoughPixel.h"
  13 #include "AliL3HoughTransformer.h"
  14 #include "AliL3HoughTrack.h"
  15 #include "AliL3TrackArray.h"
  16
  17 ClassImp(AliL3HoughTransformer)
  18
  19 AliL3HoughTransformer::AliL3HoughTransformer()
  20 {
  21   //Default constructor
  22   fTransform = 0;
  23   fEtaMin = 0;
  24   fEtaMax = 0;
  25   fSlice = 0;
  26   fPatch = 0;
  27   fRowContainer = 0;
  28   fPhiRowContainer = 0;
  29   fNumOfPadRows=0;
  30   fContainerBounds=0;
  31   fNDigits=0;
  32   fIndex = 0;
  33 }
  34
  35
  36 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
  37 {
  38   //Constructor
  39
  40   fTransform = new AliL3Transform();
  41
  42
  43   fEtaMin = etarange[0];
  44   fEtaMax = etarange[1];
  45   fSlice = slice;
  46   fPatch = patch;
  47   fNPhiSegments = phi_segments;
  48   fNumOfPadRows=NRowsSlice;
  49 }
  50
  51
  52 AliL3HoughTransformer::~AliL3HoughTransformer()
  53 {
  54   //Destructor
  55   if(fRowContainer)
  56     delete [] fRowContainer;
  57   if(fTransform)
  58     delete fTransform;
  59   if(fPhiRowContainer)
  60     delete [] fPhiRowContainer;
  61   if(fIndex)
  62     {
  63       for(Int_t i=0; i<fNDigits; i++)
  64         delete [] fIndex[i];
  65       delete [] fIndex;
  66     }
  67 }
  68
  69 void AliL3HoughTransformer::InitTemplates(TH2F *hist)
  70 {
  71
  72   AliL3Digits *pixel;
  73
  74   Int_t ymin = hist->GetYaxis()->GetFirst();
  75   Int_t ymax = hist->GetYaxis()->GetLast();
  76   Int_t nbinsy = hist->GetNbinsY();
  77
  78   fIndex = new Int_t*[fNDigits];
  79   for(Int_t i=0; i<fNDigits; i++)
  80     fIndex[i] = new Int_t[nbinsy+1];
  81
  82   Int_t sector,row;
  83   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
  84     {
  85
  86       for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
  87         {
  88           Float_t xyz[3];
  89           fTransform->Slice2Sector(fSlice,padrow,sector,row);
  90           fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
  91
  92           Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  93
  94           Double_t phi_pix = fTransform->GetPhi(xyz);
  95           Int_t index = pixel->fIndex;
  96           if(index >= fNDigits)
  97             printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
  98           for(Int_t p=ymin; p<=ymax; p++)
  99             {
 100
 101               Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
 102               Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
 103
 104               Int_t bin = hist->FindBin(kappa,phi0);
 105               if(fIndex[index][p]!=0)
 106                 printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
 107               fIndex[index][p] = bin;
 108             }
 109         }
 110     }
 111
 112 }
 113
 114
 115 void AliL3HoughTransformer::CountBins()
 116 {
 117
 118   Int_t middle_row = 87; //middle of the slice
 119
 120   Double_t r_in_bundle = fTransform->Row2X(middle_row);
 121   //  Double_t phi_min = (fSlice*20 - 10)*ToRad;
 122   //Double_t phi_max = (fSlice*20 + 10)*ToRad;
 123   Double_t phi_min = -15*ToRad;
 124   Double_t phi_max = 15*ToRad;
 125
 126   Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
 127   Double_t min_phi0 = 10000;
 128   Double_t max_phi0 = 0;
 129   Double_t min_kappa = 100000;
 130   Double_t max_kappa = 0;
 131
 132   Int_t xbin = 60;
 133   Int_t ybin = 60;
 134   Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
 135   Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
 136
 137   TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
 138
 139   for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
 140     {
 141       for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
 142         {
 143           for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
 144             {
 145               Float_t xyz[3];
 146               Int_t sector,row;
 147               fTransform->Slice2Sector(fSlice,padrow,sector,row);
 148               fTransform->Raw2Global(xyz,sector,row,pad,time);
 149               Double_t eta = fTransform->GetEta(xyz);
 150               if(eta < fEtaMin || eta > fEtaMax) continue;
 151               fTransform->Global2Local(xyz,sector);
 152               Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 153               Double_t phi_pix = fTransform->GetPhi(xyz);
 154
 155               for(Int_t p=0; p<=fNPhiSegments; p++)
 156                 {
 157                   Double_t phi_in_bundle = phi_min + p*phi_slice;
 158
 159                   Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
 160
 161                   Double_t phi0 = atan(tanPhi0);
 162                   //  if(phi0 < 0.55 || phi0 > 0.88) continue;
 163
 164                   //if(phi0 < 0) phi0 = phi0 +2*Pi;
 165                   //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
 166
 167                   Double_t angle = phi_pix - phi0;
 168                   Double_t kappa = 2*sin(angle)/r_pix;
 169                   histo->Fill(kappa,phi0,1);
 170                   if(phi0 < min_phi0)
 171                     min_phi0 = phi0;
 172                   if(phi0 > max_phi0)
 173                     max_phi0 = phi0;
 174                   if(kappa < min_kappa)
 175                     min_kappa = kappa;
 176                   if(kappa > max_kappa)
 177                     max_kappa = kappa;
 178
 179                 }
 180
 181             }
 182
 183         }
 184
 185     }
 186   Int_t count=0,bi=0;
 187
 188   Int_t xmin = histo->GetXaxis()->GetFirst();
 189   Int_t xmax = histo->GetXaxis()->GetLast();
 190   Int_t ymin = histo->GetYaxis()->GetFirst();
 191   Int_t ymax = histo->GetYaxis()->GetLast();
 192
 193
 194   for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
 195     {
 196       for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
 197         {
 198           bi++;
 199           Int_t bin = histo->GetBin(xbin,ybin);
 200           if(histo->GetBinContent(bin)>0)
 201             count++;
 202         }
 203     }
 204
 205
 206   printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
 207   printf("Phi, min %f max %f\n",min_phi0,max_phi0);
 208   printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
 209   histo->Draw("box");
 210 }
 211
 212
 213 void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
 214 {
 215   //Transformation is done with respect to local coordinates in slice.
 216
 217
 218   AliL3Digits *pix1;
 219   Int_t sector,row;
 220
 221   //Define a common point
 222
 223   /*
 224     Double_t rowdist1 = fTransform->Row2X(middle_row-1);
 225     Double_t rowdist2 = fTransform->Row2X(middle_row);
 226     Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
 227   */
 228
 229   Double_t r_in_bundle = fTransform->Row2X(middle_row);
 230   //Make overlap between slices
 231   Double_t phi_min = -15*ToRad;
 232   Double_t phi_max = 15*ToRad;
 233
 234   Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
 235
 236   for(Int_t p=0; p <= fNPhiSegments; p++)
 237     {
 238       Double_t phi_in_bundle = phi_min + p*phi_slice;
 239       //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
 240
 241       for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
 242         //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
 243         //for(Int_t padrow = 0; padrow <= middle_row; padrow++)
 244         {
 245
 246           for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
 247             {
 248
 249               Float_t xyz[3];
 250               fTransform->Slice2Sector(fSlice,padrow,sector,row);
 251               fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
 252               //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime);
 253
 254               Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 255
 256               Double_t phi_pix = fTransform->GetPhi(xyz);
 257
 258               //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
 259               Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
 260
 261               Double_t phi0 = atan(tanPhi0);
 262               //if(padrow > middle_row)
 263               //phi0 = -phi0;
 264               //if(phi0 < 0.55 || phi0 > 0.88) continue;
 265
 266               Double_t angle = phi_pix - phi0;
 267               Double_t kappa = 2*sin(angle)/r_pix;
 268
 269               //Double_t angle = phi_in_bundle - phi0;
 270               //Double_t kappa = 2*sin(angle)/r_in_bundle;
 271
 272               //if(kappa < -0.006 || kappa > 0.006) continue;
 273
 274               Short_t signal = pix1->fCharge;
 275
 276               hist->Fill(kappa,phi0,signal);
 277
 278
 279             }
 280
 281         }
 282     }
 283
 284 }
 285
 286
 287 void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
 288 {
 289   //Transformation is done with respect to local coordinates in slice.
 290   //Transform every pixel into whole phirange, using parametrisation:
 291   //kappa = 2*sin(phi-phi0)/R
 292   //Assumes you run InitTemplates firstly!!!!
 293
 294   AliL3Digits *pix1;
 295
 296   Int_t nbinsy = hist->GetNbinsY();
 297
 298   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
 299     {
 300
 301       for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
 302         {
 303
 304           Short_t signal = pix1->fCharge;
 305           Int_t index = pix1->fIndex;
 306
 307           for(Int_t p=0; p <= nbinsy; p++)
 308             hist->AddBinContent(fIndex[index][p],signal);
 309
 310         }
 311
 312     }
 313
 314
 315 }
 316
 317 /*
 318   void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
 319   {
 320   //Transformation is done with respect to local coordinates in slice.
 321   //Transform every pixel into whole phirange, using parametrisation:
 322   //kappa = 2*sin(phi-phi0)/R
 323
 324   printf("Transforming 1 pixel only\n");
 325
 326   AliL3Digits *pix1;
 327   Int_t sector,row;
 328
 329   Int_t ymin = hist->GetYaxis()->GetFirst();
 330   Int_t ymax = hist->GetYaxis()->GetLast();
 331
 332   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
 333   {
 334
 335   for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
 336   {
 337
 338   Float_t xyz[3];
 339   fTransform->Slice2Sector(fSlice,padrow,sector,row);
 340   fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
 341
 342   Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 343
 344   Double_t phi_pix = fTransform->GetPhi(xyz);
 345   Short_t signal = pix1->fCharge;
 346
 347   for(Int_t p=ymin+1; p<=ymax; p++)
 348   {
 349   Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
 350   Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
 351   hist->Fill(kappa,phi0,signal);
 352   }
 353
 354   }
 355
 356   }
 357
 358
 359   }
 360 */
 361
 362 void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
 363 {
 364
 365   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 366     {
 367       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 368       if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
 369
 370       Double_t xmin = fTransform->Row2X(track->GetFirstRow());
 371       Double_t xmax = fTransform->Row2X(track->GetLastRow());
 372
 373       Double_t a = -1./tan(track->GetPsiLine());
 374       Double_t b = track->GetDLine()/sin(track->GetPsiLine());
 375
 376       Double_t ymin = a*xmin + b;
 377       Double_t ymax = a*xmax + b;
 378
 379       Double_t middle_x = xmin + (xmax-xmin)/2;
 380       Double_t middle_y = ymin + (ymax-ymin)/2;
 381
 382       Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
 383       Double_t phi = atan2(middle_y,middle_x);
 384       Double_t phi0 = 2*phi - track->GetPsiLine();
 385       Double_t kappa = 2*sin(phi-phi0)/r_middle;
 386       hist->Fill(kappa,phi0,track->GetWeight());
 387
 388     }
 389
 390
 391 }
 392
 393 void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
 394 {
 395   //Transform every pixel into histogram, using parametrisation:
 396   //D = x*cos(psi) + y*sin(psi)
 397
 398   //  printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
 399
 400   AliL3Digits *pix1;
 401
 402   Int_t xmin = hist->GetXaxis()->GetFirst();
 403   Int_t xmax = hist->GetXaxis()->GetLast();
 404
 405   Double_t x0 = fTransform->Row2X(ref_row);
 406   Double_t y0 = 0;
 407
 408   Int_t sector,row;
 409   //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
 410
 411   Double_t phi_min = -10*ToRad;
 412   Double_t phi_max = 10*ToRad;
 413   Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
 414
 415
 416   Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
 417   Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
 418
 419
 420   Int_t index;
 421
 422   for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
 423     {
 424       for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
 425         {
 426
 427           index = phi*fNumOfPadRows + padrow;
 428           //printf("Looping index %d\n",index);
 429           if(index > fContainerBounds || index < 0)
 430             {
 431               printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
 432               return;
 433             }
 434           //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
 435           for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
 436             {
 437               //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
 438               Float_t xyz[3];
 439               fTransform->Slice2Sector(fSlice,padrow,sector,row);
 440               fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
 441
 442               if(raw)
 443                 raw->Fill(xyz[0],xyz[1],pix1->fCharge);
 444
 445               xyz[0] = xyz[0]-x0;
 446               xyz[1] = xyz[1]-y0;
 447
 448               //printf("calculating...");
 449               for(Int_t d=xmin+1; d<=xmax; d++)
 450                 {
 451                   Double_t psi = hist->GetXaxis()->GetBinCenter(d);
 452                   Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
 453
 454                   Short_t signal = pix1->fCharge;
 455                   hist->Fill(psi,D,signal);
 456                   //printf("Filling histo, psi %f D %f\n",psi,D);
 457                 }
 458               //printf("done\n");
 459             }
 460           //printf(" done\n");
 461         }
 462
 463     }
 464
 465 }
 466
 467 void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
 468 {
 469
 470   TFile *file = new TFile(rootfile);
 471   file->cd();
 472
 473   AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
 474   TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");
 475
 476   AliSimDigits da, *digarr=&da;
 477   t->GetBranch("Segment")->SetAddress(&digarr);
 478   Stat_t num_of_entries=t->GetEntries();
 479
 480   Int_t digit_counter=0;
 481   Float_t xyz[3];
 482   Double_t eta,phi;
 483
 484   Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
 485   printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
 486
 487   if(fRowContainer)
 488     delete [] fRowContainer;
 489   fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
 490   memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
 491
 492
 493   fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
 494   printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
 495   if(fPhiRowContainer)
 496     delete [] fPhiRowContainer;
 497   fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
 498   memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
 499
 500   Double_t phi_min = -10*ToRad;
 501   Double_t phi_max = 10*ToRad;
 502   Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
 503   Int_t index;
 504   digit_counter=0;
 505
 506   for (Int_t i=0; i<num_of_entries; i++)
 507     {
 508       t->GetEvent(i);
 509       Int_t sector;
 510       Int_t row;
 511       param->AdjustSectorRow(digarr->GetID(),sector,row);
 512       Int_t slice,padrow;
 513       fTransform->Sector2Slice(slice,padrow,sector,row);
 514       if(slice != fSlice) continue;
 515       if(padrow < NRows[fPatch][0]) continue;
 516       if(padrow > NRows[fPatch][1]) break;
 517       digarr->First();
 518       do {
 519         Int_t time=digarr->CurrentRow();
 520         Int_t pad=digarr->CurrentColumn();
 521         Short_t signal=digarr->CurrentDigit();
 522         if(time < param->GetMaxTBin()-1 && time > 0)
 523           if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
 524              && digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
 525             continue;
 526
 527
 528         fTransform->Raw2Global(xyz,sector,row,pad,time);
 529         eta = fTransform->GetEta(xyz);
 530         if(eta < fEtaMin || eta > fEtaMax) continue;
 531         fTransform->Global2Local(xyz,sector);
 532
 533
 534         phi = fTransform->GetPhi(xyz);
 535         if(hist)
 536           hist->Fill(xyz[0],xyz[1],signal);
 537
 538         AliL3Digits *dig = new AliL3Digits;
 539         dig->fIndex = digit_counter;
 540         digit_counter++;
 541         dig->fCharge = signal;
 542         dig->fPad = pad;
 543         dig->fTime = time;
 544
 545         if(fRowContainer[padrow].first == NULL)
 546           fRowContainer[padrow].first = (void*)dig;
 547         else
 548           ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
 549         fRowContainer[padrow].last = (void*)dig;
 550
 551         Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
 552         index = phi_index*fNumOfPadRows + padrow;
 553         if(phi_index > fContainerBounds || phi_index < 0)
 554           {
 555             printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
 556             continue;
 557           }
 558
 559         if(fPhiRowContainer[index].first == NULL)
 560           fPhiRowContainer[index].first = (void*)dig;
 561         else
 562           ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
 563         fPhiRowContainer[index].last=(void*)dig;
 564
 565       }while (digarr->Next());
 566
 567     }
 568
 569   fNDigits = digit_counter;
 570   printf("digitcounter %d\n",digit_counter);
 571   printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
 572   file->Close();
 573   delete file;
 574
 575 }