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   fParamSpace=0;
  35 }
  36
  37 AliL3HoughTransformerNew::AliL3HoughTransformerNew(Int_t slice,Int_t patch,Int_t netasegments) : AliL3HoughTransformer(slice,patch,netasegments)
  38 {
  39   fParamSpace=0;
  40 }
  41 AliL3HoughTransformerNew::~AliL3HoughTransformerNew()
  42 {
  43   if(fParamSpace)
  44     delete fParamSpace;
  45 }
  46
  47 void AliL3HoughTransformerNew::CreateHistograms(Int_t nxbins,Float_t xlow,Float_t xup,
  48                                                 Int_t nybins,Float_t ylow,Float_t yup,
  49                                                 Int_t nzbins,Float_t zlow,Float_t zup)
  50 {
  51   Char_t name[1024];
  52   sprintf(name,"paramspace_%d",(Int_t)this);
  53   fParamSpace = new TH3F(name,"",nxbins,xlow,xup,nybins,ylow,yup,nzbins,zlow,zup);
  54 }
  55
  56 void AliL3HoughTransformerNew::Reset()
  57 {
  58   fParamSpace->Reset();
  59 }
  60
  61
  62
  63 void AliL3HoughTransformerNew::TransformLine(Int_t *rowrange,Float_t *phirange)
  64 {
  65   AliL3DigitRowData *tempPt = GetDataPointer();
  66   if(!tempPt)
  67     {
  68       LOG(AliL3Log::kError,"AliL3HoughTransformerNew::TransformLine","Data")
  69         <<"No input data "<<ENDLOG;
  70       return;
  71     }
  72
  73   TH3 *hist = fParamSpace;
  74   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
  75     {
  76       AliL3DigitData *digPt = tempPt->fDigitData;
  77       if(i < rowrange[0])
  78         {
  79           AliL3MemHandler::UpdateRowPointer(tempPt);
  80           continue;
  81         }
  82       else if(i > rowrange[1])
  83         break;
  84       if(i != (Int_t)tempPt->fRow)
  85         {
  86           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
  87           continue;
  88         }
  89       for(UInt_t j=0; j<tempPt->fNDigit; j++)
  90         {
  91           UShort_t charge = digPt[j].fCharge;
  92           UChar_t pad = digPt[j].fPad;
  93           UShort_t time = digPt[j].fTime;
  94           if(charge < GetLowerThreshold() || charge > GetUpperThreshold()) continue;
  95
  96           Int_t sector,row;
  97           Float_t xyz[3];
  98           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
  99           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 100
 101           Float_t phi = AliL3Transform::GetPhi(xyz);
 102           if(phi < phirange[0] || phi > phirange[1])
 103             continue;
 104
 105           xyz[0] = xyz[0] - AliL3Transform::Row2X(rowrange[0]);
 106           Float_t x = xyz[0] + AliL3Transform::Row2X(rowrange[0]);
 107           Float_t R = sqrt(x*x + xyz[1]*xyz[1]);
 108           Float_t delta = atan(xyz[2]/R);
 109           for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
 110             {
 111               Float_t theta = hist->GetXaxis()->GetBinCenter(xbin);
 112               Float_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 113               hist->Fill(theta,rho,delta,charge);
 114             }
 115         }
 116       AliL3MemHandler::UpdateRowPointer(tempPt);
 117     }
 118
 119 }
 120
 121 struct Digit {
 122   Int_t row;
 123   Int_t charge;
 124   Float_t y;
 125   Float_t z;
 126   Digit *next;
 127 };
 128
 129 struct EtaContainer {
 130   void *first;
 131   void *last;
 132 };
 133
 134 void AliL3HoughTransformerNew::TransformLineC(Int_t *rowrange,Float_t *phirange)
 135 {
 136   AliL3DigitRowData *tempPt = GetDataPointer();
 137   if(!tempPt)
 138     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 139       <<"No input data "<<ENDLOG;
 140
 141   Int_t counter=0;
 142   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 143     {
 144       counter += tempPt->fNDigit;
 145       AliL3MemHandler::UpdateRowPointer(tempPt);
 146     }
 147
 148   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
 149   EtaContainer *etaPt = new EtaContainer[bound];
 150   memset(etaPt,0,bound*sizeof(EtaContainer));
 151
 152   Digit *digits = new Digit[counter];
 153   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 154   memset(digits,0,counter*sizeof(Digit));
 155
 156   Int_t sector,row;
 157   Float_t xyz[3];
 158
 159   counter=0;
 160   tempPt = GetDataPointer();
 161
 162   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
 163   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 164     {
 165       AliL3DigitData *digPt = tempPt->fDigitData;
 166       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 167         {
 168           Int_t charge = digPt[di].fCharge;
 169           Int_t pad = digPt[di].fPad;
 170           Int_t time = digPt[di].fTime;
 171           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 172           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 173           Double_t eta = AliL3Transform::GetEta(xyz);
 174
 175           Float_t phi = atan2(xyz[1],xyz[0]);
 176           if(phi < phirange[0] || phi > phirange[1]) continue;
 177
 178           digits[counter].row = i;
 179           digits[counter].y = xyz[1];
 180           digits[counter].z = xyz[2];
 181           digits[counter].charge = charge;
 182
 183           Int_t eta_index = GetEtaIndex(eta);
 184           Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index;
 185
 186           if(index > 0 && index < bound)
 187             {
 188               if(etaPt[index].first == 0)
 189                 etaPt[index].first = (void*)(&digits[counter]);
 190               else
 191                 ((Digit*)(etaPt[index].last))->next = &digits[counter];
 192               etaPt[index].last = (void*)(&digits[counter]);
 193             }
 194           counter++;
 195         }
 196       AliL3MemHandler::UpdateRowPointer(tempPt);
 197     }
 198
 199   cout<<"Doing the combinatorics"<<endl;
 200
 201   Digit *dPt1,*dPt2;
 202   TH3 *hist = fParamSpace;
 203   for(Int_t e=0; e<GetNEtaSegments(); e++)
 204     {
 205       for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
 206         {
 207           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
 208
 209           for(dPt1 = (Digit*)etaPt[index1].first; dPt1 != 0; dPt1 = (Digit*)dPt1->next)
 210             {
 211               for(Int_t j=i+1; j<=rowrange[1]; j++)
 212                 {
 213                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
 214
 215                   for(dPt2 = (Digit*)etaPt[index2].first; dPt2 != 0; dPt2 = (Digit*)dPt2->next)
 216                     {
 217                       if(dPt1->row == dPt2->row)
 218                         {
 219                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
 220                           exit(5);
 221                         }
 222
 223                       //Do the transform:
 224                       Float_t x1 = AliL3Transform::Row2X(dPt1->row) - AliL3Transform::Row2X(rowrange[0]);
 225                       Float_t x2 = AliL3Transform::Row2X(dPt2->row) - AliL3Transform::Row2X(rowrange[0]);
 226                       Float_t y1 = dPt1->y;
 227                       Float_t y2 = dPt2->y;
 228                       Float_t theta = atan2(x2-x1,y1-y2);
 229                       Float_t rho = x1*cos(theta)+y1*sin(theta);
 230                       Float_t R1 = sqrt(pow(AliL3Transform::Row2X(dPt1->row),2) + pow(y1,2));
 231                       Float_t delta1 = atan2(dPt1->z,R1);
 232                       Float_t R2 = sqrt(pow(AliL3Transform::Row2X(dPt2->row),2) + pow(y2,2));
 233                       Float_t delta2 = atan2(dPt2->z,R2);
 234                       Float_t delta = (dPt1->charge*delta1+dPt2->charge*delta2)/(dPt1->charge+dPt2->charge);
 235                       hist->Fill(theta,rho,delta,dPt1->charge+dPt2->charge);
 236                     }
 237                 }
 238             }
 239         }
 240     }
 241
 242   cout<<"done"<<endl;
 243   delete [] etaPt;
 244   delete [] digits;
 245
 246 }
 247