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