Small bugfix concerning calculation of eta.
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Aug 2002 08:41:00 +0000 (08:41 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Aug 2002 08:41:00 +0000 (08:41 +0000)
HLT/hough/AliL3HoughTransformerLUT.cxx
HLT/hough/AliL3HoughTransformerLUT.h

index 7d7239c..ca28f4e 100644 (file)
@@ -33,6 +33,7 @@ AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer
   fNEtas=0;
   fNPhi0=0;
   fSector=0;
+  fSlice=0;
   fSectorRow=0;
   fZSign=0;
   fZLengthPlusOff=0;
@@ -50,7 +51,26 @@ AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer
 
 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
 {
-  AliL3HoughTransformerLUT();
+  fMinRow=0;
+  fMaxRow=0;
+
+  fNRows=0;
+  fNEtas=0;
+  fNPhi0=0;
+  fSector=0;
+  fSectorRow=0;
+  fZSign=0;
+  fZLengthPlusOff=0;
+  fTimeWidth=0;
+  fPadPitch=0;
+  fEtaSlice=0;
+
+  fLUTX=0;
+  fLUTY=0;
+  fLUTEta=0;
+  fLUTphi0=0;
+  fLUT2sinphi0=0;
+  fLUT2cosphi0=0;
 
   Init(slice,patch,n_eta_segments);
 }
@@ -62,7 +82,7 @@ AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
 #ifdef do_mc
   if(fTrackID)
     {
-      for(Int_t i=0; i<GetNEtaSegments(); i++)
+      for(Int_t i=0; i<fNEtas; i++)
        {
          if(!fTrackID[i]) continue;
          delete fTrackID[i];
@@ -76,6 +96,7 @@ AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
     delete[] fLUTY;
     fNRows=0;
   }
+
   if(fNEtas){
     delete[] fLUTEta;
     fNEtas=0;
@@ -95,7 +116,7 @@ void AliL3HoughTransformerLUT::DeleteHistograms()
   }
 
   if(!fParamSpace){
-    for(Int_t i=0; i<GetNEtaSegments(); i++)
+    for(Int_t i=0; i<fNEtas; i++)
       {
        if(!fParamSpace[i]) continue;
        delete fParamSpace[i];
@@ -124,9 +145,14 @@ void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments
   fMaxRow=AliL3Transform::GetLastRow(patch);;
   fNRows=AliL3Transform::GetNRows(patch);;
   fNEtas=n_eta_segments;
+  if(fNEtas!=GetNEtaSegments()) {
+    cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
+    exit(1);
+  }
 
   AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
   fZSign = slice < 18 ? 1:-1;
+  fSlice = slice;
   fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
   fTimeWidth=AliL3Transform::GetZWidth();
 
@@ -152,8 +178,9 @@ void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments
   //lookup tables for rz2s <=> etas
   fLUTEta=new Float_t[fNEtas];
   for(Int_t rr=0;rr<fNEtas;rr++){
-    fLUTEta[rr]=CalcRoverZ2(etamin_+(rr+1)*fEtaSlice);
-    //cout << rr << ": " << fLUTEta[rr] << endl;
+    Float_t eta=etamin_+(rr+1)*fEtaSlice;
+    fLUTEta[rr]=CalcRoverZ2(eta);
+    //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
   }
 }
 
@@ -189,9 +216,9 @@ void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double
     {
       AliL3Histogram *hist = fParamSpace[0];
       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
-      cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
-      fTrackID = new TrackIndex*[GetNEtaSegments()];
-      for(Int_t i=0; i<GetNEtaSegments(); i++)
+      cout<<"Allocating "<<fNEtas*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
+      fTrackID = new TrackIndex*[fNEtas];
+      for(Int_t i=0; i<fNEtas; i++)
        fTrackID[i] = new TrackIndex[ncells];
     }
 #endif
@@ -207,9 +234,9 @@ void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double
   for(Int_t i=0; i<fNPhi0; i++){
     phi0+=diff;
     fLUTphi0[i]=phi0;
-    fLUT2sinphi0[i]=Float_t(2*sin(phi0));
-    fLUT2cosphi0[i]=Float_t(2*cos(phi0));
-    cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
+    fLUT2sinphi0[i]=2.*sin(phi0);
+    fLUT2cosphi0[i]=2.*cos(phi0);
+    //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
   }  
 }
 
@@ -224,31 +251,40 @@ void AliL3HoughTransformerLUT::Reset()
       return;
     }
   
-  for(Int_t i=0; i<GetNEtaSegments(); i++)
+  for(Int_t i=0; i<fNEtas; i++)
     fParamSpace[i]->Reset();
 #ifdef do_mc
   if(fDoMC)
     {
       AliL3Histogram *hist = fParamSpace[0];
       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
-      for(Int_t i=0; i<GetNEtaSegments(); i++)
+      for(Int_t i=0; i<fNEtas; i++)
        memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
     }
 #endif
 }
 
-
 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
 {
   //Return the histogram index of the corresponding eta. 
-  Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
-  Double_t index = (eta-GetEtaMin())/etaslice;
-  return (Int_t)index;
+  
+  Float_t rz2=CalcRoverZ2(eta);
+  return FindIndex(rz2);
+}
+
+Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2)
+{
+  Int_t index=0; //could improve search through devide and conquere strategy
+  while((index<fNEtas)&&(rz2<=fLUTEta[index])){
+    index++;
+    //cout << index << ": " << rz2 << " " << fLUTEta[index] << endl;
+  }
+  return index;
 }
 
 inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
 {
-  if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
+  if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
     return 0;
   if(!fParamSpace[eta_index])
     return 0;
@@ -266,7 +302,7 @@ Float_t AliL3HoughTransformerLUT::CalcRoverZ2(Float_t eta)
 Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
 {
   Float_t rz=sqrt(roverz2);
-  if(fZSign>0) rz=-rz;
+  if(fZSign<0) rz=-rz;
   Float_t ret=(1+rz)/(rz-1);
   ret=0.5*log(ret);
   return ret;
@@ -292,43 +328,48 @@ Float_t AliL3HoughTransformerLUT::CalcZ(Int_t time)
 
 Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
 {
-  Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
-  Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
-  if(slice>17) eta*=-1;
-  return eta;
+  if(eta_index >= fNEtas || eta_index < 0){
+    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
+    return 0.;
+  }
+  if(slice != fSlice){
+    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
+    return 0.;
+  }
+
+  return (CalcEta(fLUTEta[eta_index])-0.5*fEtaSlice);
 }
 
 void AliL3HoughTransformerLUT::TransformCircle()
 {
-#if 0
   //Transform the input data with a circle HT.
-  //The function loops over all the data, and transforms each pixel with the equations:
+  //The function loops over all the data, and transforms each pixel with the equation:
   // 
-  //kappa = 2/R*sin(phi - phi0)
+  //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
   //
-  //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
+  //where R^2 = x^2 + y^2
   //
   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
   //which histogram in which the pixel should be transformed, the eta-value is calcluated
   //and the proper histogram index is found by GetEtaIndex(eta).
 
-
   AliL3DigitRowData *tempPt = GetDataPointer();
   if(!tempPt)
     {
-      LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
+      LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
        <<"No input data "<<ENDLOG;
       return;
     }
   
   //Loop over the padrows:
-  for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
+  for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
     {
       //Get the data on this padrow:
       AliL3DigitData *digPt = tempPt->fDigitData;
       if(i != (Int_t)tempPt->fRow)
        {
-         cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
+         LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
+           <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<endl;
          continue;
        }
 
@@ -336,68 +377,65 @@ void AliL3HoughTransformerLUT::TransformCircle()
       for(UInt_t j=0; j<tempPt->fNDigit; j++)
        {
          UShort_t charge = digPt[j].fCharge;
-         UChar_t pad = digPt[j].fPad;
-         UShort_t time = digPt[j].fTime;
+
+         //check threshold
          if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
            continue;
-         Int_t sector,row;
-         Float_t xyz[3];
-         
-         //Transform data to local cartesian coordinates:
-         AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
-         AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
-         
-         //Calculate the eta:
-         Double_t eta = AliL3Transform::GetEta(xyz);
-         
-         //Get the corresponding index, which determines which histogram to fill:
-         Int_t eta_index = GetEtaIndex(eta);
-         if(eta_index < 0 || eta_index >= GetNEtaSegments())
+
+         UChar_t pad = digPt[j].fPad;
+         UShort_t time = digPt[j].fTime;
+
+         Float_t x = CalcX(row);
+         Float_t y = CalcY(pad,row);
+         Float_t z = CalcZ(time);
+
+         //find eta slice
+         Float_t rz2 = 1 + (x*x+y*y)/(z*z);
+         Int_t eta_index = FindIndex(rz2);
+
+         if(eta_index < 0 || eta_index >= fNEtas){
+           //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
            continue;
-         
+         }       
+
          //Get the correct histogrampointer:
          AliL3Histogram *hist = fParamSpace[eta_index];
-         if(!hist)
-           {
-             printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
-             continue;
-           }
-
-         //Do the transformation:
-         Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
-         Float_t phi = AliL3Transform::GetPhi(xyz);
-         
+         if(!hist){
+           //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
+           continue;
+         }
+
          //Fill the histogram along the phirange
-         for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
-           {
-             Float_t phi0 = hist->GetBinCenterY(b);
-             Float_t kappa = 2*sin(phi - phi0)/R;
-             hist->Fill(kappa,phi0,charge);
+         Float_t R2=1/(x*x+y*y);
+         for(Int_t b=0; b<fNPhi0; b++){
+           Float_t kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
+
+           hist->Fill(kappa,fLUTphi0[b],charge);
+           //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
+
 #ifdef do_mc
-             if(fDoMC)
-               {
-                 Int_t bin = hist->FindBin(kappa,phi0);
-                 for(Int_t t=0; t<3; t++)
-                   {
-                     Int_t label = digPt[j].fTrackID[t];
-                     if(label < 0) break;
-                     UInt_t c;
-                     for(c=0; c<MaxTrack; c++)
-                       if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
-                         break;
-                     if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
-                     fTrackID[eta_index][bin].fLabel[c] = label;
-                     fTrackID[eta_index][bin].fNHits[c]++;
-                   }
-               }
+           if(fDoMC)
+             {
+               Int_t bin = hist->FindBin(kappa,phi0);
+               for(Int_t t=0; t<3; t++)
+                 {
+                   Int_t label = digPt[j].fTrackID[t];
+                   if(label < 0) break;
+                   UInt_t c;
+                   for(c=0; c<MaxTrack; c++)
+                     if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
+                       break;
+                   if(c == MaxTrack-1) LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","MCData") <<"Array reached maximum!! "<<c<<endl;
+                   fTrackID[eta_index][bin].fLabel[c] = label;
+                   fTrackID[eta_index][bin].fNHits[c]++;
+                 }
+             }
 #endif
-           }
+         }
        }
-      
       //Move the data pointer to the next padrow:
       AliL3MemHandler::UpdateRowPointer(tempPt);
     }
-#endif //to do later
 }
 
 Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
@@ -409,7 +447,7 @@ Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double
     }
   
 #ifdef do_mc
-  if(eta_index < 0 || eta_index > GetNEtaSegments())
+  if(eta_index < 0 || eta_index > fNEtas)
     {
       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
       return -1;
index affeb6e..2861f98 100644 (file)
@@ -22,6 +22,7 @@ class AliL3HoughTransformerLUT : public AliL3HoughBaseTransformer {
   Int_t fNRows;
   Int_t fNEtas;
   Int_t fNPhi0;
+  Int_t fSlice;
   Int_t fSector;
   Int_t fSectorRow;
   Int_t fZSign;
@@ -43,7 +44,7 @@ class AliL3HoughTransformerLUT : public AliL3HoughBaseTransformer {
   Float_t CalcY(Int_t pad, Int_t row);
   Float_t CalcZ(Int_t time);  
 
-  Int_t FindIndex(Double_t rz2);
+  Int_t FindIndex(Float_t rz2);
 
  public: