HLT/hough/AliL3HoughTransformerNew.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Logging.h"
   9 #include "AliL3HoughTransformerNew.h"
  10 #include "AliL3MemHandler.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3DigitData.h"
  13 #include "AliL3HistogramAdaptive.h"
  14 #include "AliL3TrackArray.h"
  15 #include "AliL3HoughTrack.h"
  16 #include "AliL3SpacePointData.h"
  17 #include "AliL3Vertex.h"
  18 #include "AliL3Fitter.h"
  19
  20 #if __GNUC__ >= 3
  21 using namespace std;
  22 #endif
  23
  24 //_____________________________________________________________
  25 // AliL3HoughTransformerNew
  26 //
  27 // Hough transformation class
  28 //
  29
  30 ClassImp(AliL3HoughTransformerNew)
  31
  32 AliL3HoughTransformerNew::AliL3HoughTransformerNew()
  33 {
  34   //default ctor
  35   fParamSpace3D=0;
  36 }
  37
  38 AliL3HoughTransformerNew::AliL3HoughTransformerNew(Int_t slice,Int_t patch,Int_t netasegments) : AliL3HoughTransformer(slice,patch,netasegments)
  39 {
  40   //normal ctor
  41   fParamSpace3D=0;
  42 }
  43 AliL3HoughTransformerNew::~AliL3HoughTransformerNew()
  44 {
  45   //dtor
  46   if(fParamSpace3D)
  47     delete fParamSpace3D;
  48 }
  49
  50 void AliL3HoughTransformerNew::CreateHistograms(Int_t nxbins,Float_t xlow,Float_t xup,
  51                                                 Int_t nybins,Float_t ylow,Float_t yup,
  52                                                 Int_t nzbins,Float_t zlow,Float_t zup)
  53 {
  54   //Create the histogram which contain the hough space
  55   Char_t name[1024];
  56   sprintf(name,"paramspace_%p",(void*)this);
  57   fParamSpace3D = new TH3F(name,"",nxbins,xlow,xup,nybins,ylow,yup,nzbins,zlow,zup);
  58 }
  59
  60 void AliL3HoughTransformerNew::Reset()
  61 {
  62   //Reset Hough space
  63   fParamSpace3D->Reset();
  64 }
  65
  66
  67
  68 void AliL3HoughTransformerNew::TransformLine(Int_t *rowrange,Float_t *phirange)
  69 {
  70   //Hough Transform
  71   AliL3DigitRowData *tempPt = GetDataPointer();
  72   if(!tempPt)
  73     {
  74       LOG(AliL3Log::kError,"AliL3HoughTransformerNew::TransformLine","Data")
  75         <<"No input data "<<ENDLOG;
  76       return;
  77     }
  78
  79   TH3 *hist = fParamSpace3D;
  80   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
  81     {
  82       AliL3DigitData *digPt = tempPt->fDigitData;
  83       if(i < rowrange[0])
  84         {
  85           AliL3MemHandler::UpdateRowPointer(tempPt);
  86           continue;
  87         }
  88       else if(i > rowrange[1])
  89         break;
  90       if(i != (Int_t)tempPt->fRow)
  91         {
  92           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
  93           continue;
  94         }
  95       for(UInt_t j=0; j<tempPt->fNDigit; j++)
  96         {
  97           UShort_t charge = digPt[j].fCharge;
  98           UChar_t pad = digPt[j].fPad;
  99           UShort_t time = digPt[j].fTime;
 100           if(charge < GetLowerThreshold() || charge > GetUpperThreshold()) continue;
 101
 102           Int_t sector,row;
 103           Float_t xyz[3];
 104           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 105           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 106
 107           Float_t phi = AliL3Transform::GetPhi(xyz);
 108           if(phi < phirange[0] || phi > phirange[1])
 109             continue;
 110
 111           xyz[0] = xyz[0] - AliL3Transform::Row2X(rowrange[0]);
 112           Float_t x = xyz[0] + AliL3Transform::Row2X(rowrange[0]);
 113           Float_t r = sqrt(x*x + xyz[1]*xyz[1]);
 114           Float_t delta = atan(xyz[2]/r);
 115           for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
 116             {
 117               Float_t theta = hist->GetXaxis()->GetBinCenter(xbin);
 118               Float_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 119               hist->Fill(theta,rho,delta,charge);
 120             }
 121         }
 122       AliL3MemHandler::UpdateRowPointer(tempPt);
 123     }
 124
 125 }
 126
 127 struct AliL3Digit {
 128   Int_t fRow;//padrow index
 129   Int_t fCharge;//digits charge
 130   Float_t fY;//Y coordinate of the digit
 131   Float_t fZ;//Z coordinate of the digit
 132   AliL3Digit *fNext;//pointer to the next digit
 133 };
 134
 135 struct AliL3EtaContainer {
 136   void *fFirst;//pointer to the first digit in a sequence
 137   void *fLast;//pointer to the last digit in a sequence
 138 };
 139
 140 void AliL3HoughTransformerNew::TransformLineC(Int_t *rowrange,Float_t *phirange)
 141 {
 142   //Hough Transform
 143   AliL3DigitRowData *tempPt = GetDataPointer();
 144   if(!tempPt)
 145     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 146       <<"No input data "<<ENDLOG;
 147
 148   Int_t counter=0;
 149   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 150     {
 151       counter += tempPt->fNDigit;
 152       AliL3MemHandler::UpdateRowPointer(tempPt);
 153     }
 154
 155   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
 156   AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
 157   memset(etaPt,0,bound*sizeof(AliL3EtaContainer));
 158
 159   AliL3Digit *digits = new AliL3Digit[counter];
 160   cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
 161   memset(digits,0,counter*sizeof(AliL3Digit));
 162
 163   Int_t sector,row;
 164   Float_t xyz[3];
 165
 166   counter=0;
 167   tempPt = GetDataPointer();
 168
 169   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
 170   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 171     {
 172       AliL3DigitData *digPt = tempPt->fDigitData;
 173       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 174         {
 175           Int_t charge = digPt[di].fCharge;
 176           Int_t pad = digPt[di].fPad;
 177           Int_t time = digPt[di].fTime;
 178           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 179           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 180           Double_t eta = AliL3Transform::GetEta(xyz);
 181
 182           Float_t phi = atan2(xyz[1],xyz[0]);
 183           if(phi < phirange[0] || phi > phirange[1]) continue;
 184
 185           digits[counter].fRow = i;
 186           digits[counter].fY = xyz[1];
 187           digits[counter].fZ = xyz[2];
 188           digits[counter].fCharge = charge;
 189
 190           Int_t etaindex = GetEtaIndex(eta);
 191           Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
 192
 193           if(index > 0 && index < bound)
 194             {
 195               if(etaPt[index].fFirst == 0)
 196                 etaPt[index].fFirst = (void*)(&digits[counter]);
 197               else
 198                 ((AliL3Digit*)(etaPt[index].fLast))->fNext = &digits[counter];
 199               etaPt[index].fLast = (void*)(&digits[counter]);
 200             }
 201           counter++;
 202         }
 203       AliL3MemHandler::UpdateRowPointer(tempPt);
 204     }
 205
 206   cout<<"Doing the combinatorics"<<endl;
 207
 208   AliL3Digit *dPt1,*dPt2;
 209   TH3 *hist = fParamSpace3D;
 210   for(Int_t e=0; e<GetNEtaSegments(); e++)
 211     {
 212       for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
 213         {
 214           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
 215
 216           for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
 217             {
 218               for(Int_t j=i+1; j<=rowrange[1]; j++)
 219                 {
 220                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
 221
 222                   for(dPt2 = (AliL3Digit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit*)dPt2->fNext)
 223                     {
 224                       if(dPt1->fRow == dPt2->fRow)
 225                         {
 226                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
 227                           exit(5);
 228                         }
 229
 230                       //Do the transform:
 231                       Float_t x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
 232                       Float_t x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
 233                       Float_t y1 = dPt1->fY;
 234                       Float_t y2 = dPt2->fY;
 235                       Float_t theta = atan2(x2-x1,y1-y2);
 236                       Float_t rho = x1*cos(theta)+y1*sin(theta);
 237                       Float_t r1 = sqrt(pow(AliL3Transform::Row2X(dPt1->fRow),2) + pow(y1,2));
 238                       Float_t delta1 = atan2(dPt1->fZ,r1);
 239                       Float_t r2 = sqrt(pow(AliL3Transform::Row2X(dPt2->fRow),2) + pow(y2,2));
 240                       Float_t delta2 = atan2(dPt2->fZ,r2);
 241                       Float_t delta = (dPt1->fCharge*delta1+dPt2->fCharge*delta2)/(dPt1->fCharge+dPt2->fCharge);
 242                       hist->Fill(theta,rho,delta,dPt1->fCharge+dPt2->fCharge);
 243                     }
 244                 }
 245             }
 246         }
 247     }
 248
 249   cout<<"done"<<endl;
 250   delete [] etaPt;
 251   delete [] digits;
 252
 253 }
 254