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