Checking in the seeds of new cluster fitting code.
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Nov 2002 12:33:46 +0000 (12:33 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Nov 2002 12:33:46 +0000 (12:33 +0000)
HLT/comp/AliL3ClusterFitter.cxx [new file with mode: 0644]
HLT/comp/AliL3ClusterFitter.h [new file with mode: 0644]
HLT/comp/AliL3FitUtilities.c [new file with mode: 0644]
HLT/comp/AliL3FitUtilities.h [new file with mode: 0644]

diff --git a/HLT/comp/AliL3ClusterFitter.cxx b/HLT/comp/AliL3ClusterFitter.cxx
new file mode 100644 (file)
index 0000000..1f35dc7
--- /dev/null
@@ -0,0 +1,489 @@
+//$Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV
+
+#include "AliL3StandardIncludes.h"
+#include "AliL3ClusterFitter.h"
+#include "AliL3FitUtilities.h"
+#include "AliL3Transform.h"
+#include "AliL3DigitData.h"
+#include "AliL3ModelTrack.h"
+#include "AliL3TrackArray.h"
+#include "AliL3MemHandler.h"
+
+#if GCCVERSION == 3
+using namespace std;
+#endif
+
+//_____________________________________________________________
+//
+//  AliL3ClusterFitter
+//
+
+
+ClassImp(AliL3ClusterFitter)
+
+AliL3ClusterFitter::AliL3ClusterFitter()
+{
+  fPadFitRange=2;
+  fTimeFitRange=3;
+  plane=0;
+  fNmaxOverlaps = 3;
+  fChiSqMax=12;
+}
+
+AliL3ClusterFitter::~AliL3ClusterFitter()
+{
+
+}
+
+void AliL3ClusterFitter::FindClusters()
+{
+  if(!fTracks)
+    {
+      cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
+      return;
+    }
+  if(!fRowData)
+    {
+      cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
+      return;
+    }
+  
+  AliL3DigitRowData *rowPt = fRowData;
+  AliL3DigitData *digPt=0;
+
+  Int_t pad,time;
+  Short_t charge;
+
+  for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
+    {
+      fCurrentPadRow = i;
+      memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
+      digPt = (AliL3DigitData*)rowPt->fDigitData;
+      //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
+      for(UInt_t j=0; j<rowPt->fNDigit; j++)
+       {
+         pad = digPt[j].fPad;
+         time = digPt[j].fTime;
+         charge = digPt[j].fCharge;
+         fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
+         fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
+         //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
+       }
+      
+      for(Int_t k=0; k<fTracks->GetNTracks(); k++)
+       {
+         AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
+         if(!track) continue;
+         if(track->GetPadHit(i) < 0 || track->GetTimeHit(i) < 0)
+           {
+             track->SetCluster(i,0,0,0,0,0,0);
+             continue;
+           }
+         FitClusters(track);
+       }
+      FillZeros(rowPt);
+      AliL3MemHandler::UpdateRowPointer(rowPt);
+    }
+}
+
+void AliL3ClusterFitter::FitCluster(AliL3ModelTrack *track)
+{
+  cout<<"Track had no overlaps"<<endl;
+  Int_t minpad = (Int_t)rint(track->GetPadHit(fCurrentPadRow)) - fPadFitRange;
+  Int_t maxpad = (Int_t)rint(track->GetPadHit(fCurrentPadRow)) + fPadFitRange;
+  Int_t mintime = (Int_t)rint(track->GetTimeHit(fCurrentPadRow)) - fTimeFitRange;
+  Int_t maxtime = (Int_t)rint(track->GetTimeHit(fCurrentPadRow)) + fTimeFitRange;
+  
+  Int_t size = FIT_PTS;
+
+  if(minpad <= 0)
+    minpad=0;
+  if(mintime<=0)
+    mintime=0;
+  if(maxpad>=AliL3Transform::GetNPads(fCurrentPadRow)-1)
+    maxpad=AliL3Transform::GetNPads(fCurrentPadRow)-1;
+  if(maxtime>=AliL3Transform::GetNTimeBins()-1)
+    maxtime=AliL3Transform::GetNTimeBins()-1;
+
+  if(plane)
+    delete [] plane;
+  plane = new DPOINT[FIT_PTS];
+  memset(plane,0,FIT_PTS*sizeof(DPOINT));
+
+  Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
+  
+  //Fill the fit parameters:
+  Double_t a[FIT_MAXPAR];
+  a[2] = track->GetPadHit(fCurrentPadRow);
+  a[4] = track->GetTimeHit(fCurrentPadRow);
+  a[1] = fRow[(AliL3Transform::GetNTimeBins()+1)*((Int_t)rint(a[2])) + (Int_t)rint(a[4])].fCharge;
+  a[3] = sqrt(track->GetParSigmaY2(fCurrentPadRow))*sqrt(2);
+  a[5] = sqrt(track->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
+  a[6] = sqrt(track->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
+  
+  if(!a[1])
+    {
+      cout<<"No charge"<<endl;
+      if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+       track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
+      return;
+    }
+  
+  Int_t pad_num=0;
+  Int_t time_num_max=0;
+  Int_t ndata=0;
+  Int_t charge,tot_charge=0;
+
+  //Fill the proper arrays:
+  for(Int_t i=minpad; i<=maxpad; i++)
+    {
+      Int_t max_charge = 0;
+      Int_t time_num=0;
+      for(Int_t j=mintime; j<=maxtime; j++)
+       {
+         charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
+         if(charge > max_charge)
+           {
+             max_charge = charge;
+             time_num++;
+            }
+         if(charge <= 0) continue;
+         //cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
+         tot_charge += charge;
+         ndata++;
+         if(ndata >= size)
+           cerr<<"Too many points"<<endl;
+         plane[ndata].u = (Double_t)i;
+         plane[ndata].v = (Double_t)j;
+         x[ndata]=ndata;
+         y[ndata]=charge;
+         s[ndata]= 1 + sqrt((Double_t)charge);
+       }
+      if(max_charge) //there was charge on this pad
+       pad_num++;
+      if(time_num_max < time_num)
+       time_num_max = time_num;
+    }
+
+  if(pad_num <= 1 || time_num_max <=1) //too few to do fit
+    {
+      //cout<<"To few pads for fit on row "<<fCurrentPadRow<<endl;
+      if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+       track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
+
+      return;
+    }
+  
+  Int_t npars = NUM_PARS;
+  Int_t lista[FIT_MAXPAR];
+  Double_t dev[FIT_MAXPAR],chisq_f;
+  lista[1] = 1;
+  lista[2] = 1;
+  lista[3] = 0;
+  lista[4] = 1;
+  lista[5] = 0;
+  lista[6] = 0;
+  
+  //Declare a pointer to the C fitting function
+  void (*funcs) ( double, double *, double *,double *,int );
+  funcs = f2gauss5;
+  
+  //cout<<"Doing fit with parameters "<<a[2]<<" "<<a[4]<<" "<<a[3]<<" "<<a[5]<<" "<<a[1]<<endl;
+  
+  Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f,f2gauss5);
+  if(ret)
+    cerr<<"Fit error"<<endl;
+  
+  tot_charge = (Int_t)(a[1] * a[3] * a[5]);
+  chisq_f /= (ndata - 3);
+  //  cout<<"Chisq per degree of freedom : "<<chisq_f<<endl;
+  //cout<<"pad "<<a[2]<<" time "<<a[4]<<" adc "<<a[1]<<endl;
+  //cout<<"Setting cluster on row "<<fCurrentPadRow<<" pad "<<a[2]<<" time "<<a[4]<<" charge "<<tot_charge<<endl;
+  
+  //Make sure that the cluster has not already been set before setting it:
+  if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+    track->SetCluster(fCurrentPadRow,a[2],a[4],tot_charge,0,0,pad_num);
+}
+
+void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track)
+{
+  //Handle single and overlapping clusters
+    
+  if(!track->IsPresent(fCurrentPadRow) && track->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+    {    
+      if(!track->IsPresent(fCurrentPadRow) && 
+        track->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch) + 1)//debug
+       {
+         cerr<<"AliL3ClusterFitter::FitClusters() : Mismatching clustercount "
+             <<track->GetNClusters()<<" "<<fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch) + 1<<endl;
+         exit(5);
+       }
+      return; //This cluster has been set before.
+    }
+  
+  
+  Int_t minpad,maxpad,mintime,maxtime;
+  
+  Int_t size = FIT_PTS;
+  Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
+  if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
+    {
+      cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
+      return;
+    }
+  Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
+  
+  //Check if at least one cluster is not already fitted
+  Bool_t all_fitted=kTRUE;
+  Int_t k=-1;
+  while(k < track->GetNOverlaps(fCurrentPadRow))
+    {
+      AliL3ModelTrack *tr=0;
+      if(k==-1)
+       tr = track;
+      else
+       tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
+      k++;
+      if(!tr) continue;
+      if(!tr->IsPresent(fCurrentPadRow))
+       {
+         all_fitted = kFALSE;
+         break;
+       }
+    }
+  if(all_fitted)
+    {
+      cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
+      return;
+    }
+  
+  plane = new DPOINT[FIT_PTS];
+  memset(plane,0,FIT_PTS*sizeof(DPOINT));
+
+  Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
+  
+  //Fill the fit parameters:
+  Double_t a[FIT_MAXPAR];
+  Int_t lista[FIT_MAXPAR];
+  Double_t dev[FIT_MAXPAR],chisq_f;
+  
+  minpad = mintime = 999;
+  maxpad = maxtime = 0;
+  Int_t fit_pars=0;
+  
+  Int_t ntracks = track->GetNOverlaps(fCurrentPadRow)+1;
+  Bool_t *do_fit = new Bool_t[ntracks];
+  for(Int_t i=0; i<ntracks; i++)
+    do_fit[i]=kTRUE;
+  
+  Int_t n_overlaps=0;
+  k=-1;
+  //Fill the overlapping tracks:
+  while(k < track->GetNOverlaps(fCurrentPadRow))
+    {
+      AliL3ModelTrack *tr=0;
+      if(k==-1)
+       tr = track;
+      else
+       tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
+      k++;
+      if(!tr) continue;
+      
+      Int_t hitpad = (Int_t)rint(tr->GetPadHit(fCurrentPadRow));
+      Int_t hittime = (Int_t)rint(tr->GetTimeHit(fCurrentPadRow));
+      Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
+      if(!charge) //There is not charge here, so the cluster is non-existing -- remove it.
+       {
+         if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+           tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);           
+         do_fit[k] = kFALSE;
+         continue;
+       }
+      
+      if(k==0)
+       LocateCluster(tr,minpad,maxpad,mintime,maxtime);
+
+      if(tr->GetPadHit(fCurrentPadRow) < (Double_t)minpad || tr->GetPadHit(fCurrentPadRow) > (Double_t)maxpad ||
+        tr->GetTimeHit(fCurrentPadRow) < (Double_t)mintime || tr->GetTimeHit(fCurrentPadRow) > (Double_t)maxtime)
+       {
+         do_fit[k] = kFALSE;//This cluster is outside the region already specified, so it will not be included in this fit.
+
+         //If this is the first: remove it, because it will not be checked again.
+         if(k==0)
+           if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+             tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);                 
+         continue;
+       }
+      
+      
+      cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
+         <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
+         <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
+      
+      a[n_overlaps*NUM_PARS+2] = tr->GetPadHit(fCurrentPadRow);
+      a[n_overlaps*NUM_PARS+4] = tr->GetTimeHit(fCurrentPadRow);
+      
+      if(!tr->IsPresent(fCurrentPadRow)) //Cluster is not fitted before
+       {
+         a[n_overlaps*NUM_PARS+1] = charge;
+         a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow))*sqrt(2);
+         a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
+         a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
+         lista[n_overlaps*NUM_PARS + 1] = 1;
+         lista[n_overlaps*NUM_PARS + 2] = 1;
+         lista[n_overlaps*NUM_PARS + 3] = 0;
+         lista[n_overlaps*NUM_PARS + 4] = 1;
+         lista[n_overlaps*NUM_PARS + 5] = 0;
+         lista[n_overlaps*NUM_PARS + 6] = 0;
+         fit_pars             += 3;
+       }
+      else
+       {
+         Int_t charge;
+         Float_t xywidth,zwidth,pad,time;
+         tr->GetPad(fCurrentPadRow,pad);
+         tr->GetTime(fCurrentPadRow,time);
+         tr->GetClusterCharge(fCurrentPadRow,charge);
+         xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
+         zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
+         cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
+         
+         a[n_overlaps*NUM_PARS+2] = pad;
+         a[n_overlaps*NUM_PARS+4] = time;
+         a[n_overlaps*NUM_PARS+1] = charge;
+         a[n_overlaps*NUM_PARS+3] = sqrt(xywidth)*sqrt(2);
+         a[n_overlaps*NUM_PARS+5] = sqrt(zwidth)*sqrt(2);
+         a[n_overlaps*NUM_PARS+6] = sqrt(zwidth)*sqrt(2);
+
+         lista[n_overlaps*NUM_PARS + 1] = 1;
+         lista[n_overlaps*NUM_PARS + 2] = 0;
+         lista[n_overlaps*NUM_PARS + 3] = 0;
+         lista[n_overlaps*NUM_PARS + 4] = 0;
+         lista[n_overlaps*NUM_PARS + 5] = 0;
+         lista[n_overlaps*NUM_PARS + 6] = 0;
+         fit_pars             += 1;
+       }
+      n_overlaps++;
+    }
+  
+  if(n_overlaps==0) //No clusters here
+    return;
+  cout<<"Setting init searchrange; pad "<<minpad<<" "<<maxpad<<" time "<<mintime<<" "<<maxtime<<endl;
+
+  if(minpad <= 0)
+    minpad=0;
+  if(mintime<=0)
+    mintime=0;
+  if(maxpad>=AliL3Transform::GetNPads(fCurrentPadRow)-1)
+    maxpad=AliL3Transform::GetNPads(fCurrentPadRow)-1;
+  if(maxtime>=AliL3Transform::GetNTimeBins()-1)
+    maxtime=AliL3Transform::GetNTimeBins()-1;
+  
+  Int_t pad_num=0;
+  Int_t time_num_max=0;
+  Int_t ndata=0;
+  Int_t tot_charge=0;
+
+  for(Int_t i=minpad; i<=maxpad; i++)
+    {
+      Int_t max_charge = 0;
+      Int_t time_num=0;
+      for(Int_t j=mintime; j<=maxtime; j++)
+       {
+         Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
+         
+         if(charge <= 0) continue;
+
+         if(charge > max_charge)
+           {
+             max_charge = charge;
+             time_num++;
+           }
+         cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
+         tot_charge += charge;
+         ndata++;
+         if(ndata >= size)
+           cerr<<"Too many points"<<endl;
+         
+         //This digit will most likely be used:
+         fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
+         plane[ndata].u = (Double_t)i;
+         plane[ndata].v = (Double_t)j;
+         x[ndata]=ndata;
+         y[ndata]=charge;
+         s[ndata]= 1 + sqrt((Double_t)charge);
+       }
+      if(max_charge) //there was charge on this pad
+       pad_num++;
+      if(time_num_max < time_num)
+       time_num_max = time_num;
+    }
+  
+  k=-1;
+  if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps) //too few to do fit
+    {
+      while(k < track->GetNOverlaps(fCurrentPadRow))
+       {
+         AliL3ModelTrack *tr=0;
+         if(k==-1)
+           tr = track;
+         else
+           tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
+         k++;
+         if(!tr) continue;
+         if(do_fit[k] == kFALSE) continue;
+         
+         if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
+           tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
+       }
+      cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<endl;
+      return;
+    }
+  
+  Int_t npars = n_overlaps * NUM_PARS;
+  cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
+  Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
+  if(ret)
+    {
+      cerr<<"Fit error"<<endl;
+      exit(5);
+    }
+
+  chisq_f /= (ndata-fit_pars);
+  cout<<"Chisq "<<chisq_f<<endl;
+
+  k=-1;
+  n_overlaps=0;
+  while(k < track->GetNOverlaps(fCurrentPadRow))
+    {
+      AliL3ModelTrack *tr=0;
+      if(k==-1)
+       tr = track;
+      else
+       tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
+      k++;
+      if(!tr) continue;
+      if(do_fit[k] == kFALSE) continue;
+      
+      if(!tr->IsPresent(fCurrentPadRow))
+       {
+         if(tr->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch)) continue;//This cluster has been set before
+         if(chisq_f < fChiSqMax)//cluster fit is good enough
+           {
+             tot_charge = (Int_t)(a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
+             tr->SetCluster(fCurrentPadRow,a[n_overlaps*NUM_PARS+2],a[n_overlaps*NUM_PARS+4],tot_charge,0,0,pad_num);
+             cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]<<" charge "<<tot_charge<<endl;
+           }
+         else //fit was too bad
+           tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
+       }
+      n_overlaps++;
+    }
+  
+  delete [] plane;
+  delete [] do_fit;
+}
+
diff --git a/HLT/comp/AliL3ClusterFitter.h b/HLT/comp/AliL3ClusterFitter.h
new file mode 100644 (file)
index 0000000..e7ae5b1
--- /dev/null
@@ -0,0 +1,33 @@
+#ifndef AliL3_ClusterFitter
+#define AliL3_ClusterFitter
+
+#include "AliL3RootTypes.h"
+#include "AliL3Modeller.h"
+
+class AliL3ModelTrack;
+
+class AliL3ClusterFitter : public AliL3Modeller {
+  
+ private:
+  Int_t fPadFitRange;
+  Int_t fTimeFitRange;
+  Int_t fNmaxOverlaps;
+  Float_t fChiSqMax;
+  
+  void FitCluster(AliL3ModelTrack *track);
+  void FitClusters(AliL3ModelTrack *track);
+  
+ public:
+  AliL3ClusterFitter();
+  virtual ~AliL3ClusterFitter();
+  
+  void FindClusters();
+
+  void SetFitRange(Int_t p,Int_t t) {fPadFitRange=p; fTimeFitRange=t;}
+  void SetNmaxOverlaps(Int_t i) {fNmaxOverlaps=i;}
+
+  ClassDef(AliL3ClusterFitter,1) 
+
+};
+
+#endif
diff --git a/HLT/comp/AliL3FitUtilities.c b/HLT/comp/AliL3FitUtilities.c
new file mode 100644 (file)
index 0000000..1a829dc
--- /dev/null
@@ -0,0 +1,446 @@
+
+#include "AliL3FitUtilities.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <setjmp.h>
+
+
+jmp_buf env;
+
+DPOINT *plane;
+
+#ifndef PB
+int build_exp_table()
+/******************************************************************************/
+{
+       exp_header_t tab_h;
+       int i;
+       FLOAT_SIZE *exp_val;
+       FILE  *fp;
+       
+       if( !(fp = fopen( "exp_table", "wb" )))
+         printf("build_exp_table: I/O error\n");
+               
+       tab_h.float_size = sizeof( FLOAT_SIZE );
+       tab_h.steps     = 100000;
+       tab_h.dmin      = -5;
+       tab_h.dmax      = 1;
+       tab_h.step_size = (tab_h.dmax - tab_h.dmin)/tab_h.steps;
+       if( !(exp_val = (FLOAT_SIZE *)malloc( tab_h.steps * sizeof( FLOAT_SIZE ) ) ) )
+         printf("build_exp_table: malloc error\n");
+       fwrite((const void *)&tab_h, (size_t)sizeof(exp_header_t), (size_t)1, fp );
+       for( i=0; i<tab_h.steps; ++i ) {
+               exp_val[i] = exp(tab_h.dmin + i*tab_h.step_size);
+       }
+       fwrite((const void *)exp_val, (size_t)sizeof(FLOAT_SIZE), (size_t)tab_h.steps, fp );
+       
+       free( (void *)exp_val );
+       return 0;
+}
+
+#ifdef TEST_TABLE
+test_table()
+/**********************************************************/
+{
+       double i, delta, dmin, dmax, exp_tab();
+       FILE   *fp;
+       
+       fp = fopen( "exp_test", "w" );
+       
+       delta = 6.0/10000.;
+       dmin  = -5.1;
+       dmax  =  1.1;
+       for( i=dmin; i<dmax; i+=delta ){
+               fprintf( fp, "%lf\t%lf\t%lf\n", i, exp(i), exp_tab(i) );
+       }
+       return 0;
+}
+#endif
+double exp_tab(double in)
+{
+       return exp( in );
+}
+#else
+double exp_tab(in)
+/**********************************************************/
+double in;
+{
+       static exp_header_t tab_h;
+       static FLOAT_SIZE   *exp_val=NULL, slope;
+       FILE                *fp=NULL;
+       
+#ifdef HP
+       if( isnan(in) ) {
+               if( (num_nan /100)*100 == num_nan )
+                       fprintf( stderr, "exp_tab: NaN %d\n", num_nan );
+               ++num_nan;
+               return 1;
+       }
+#endif
+#ifdef MAC
+               return exp( in );
+#else
+       if( !exp_val ) {
+               if( !(fp = fopen( "exp_table", "rb" ))) {
+                       build_exp_table();
+               }
+               if( !(fp = fopen( "exp_table", "rb" )))
+                       printf("exp_tab: I/O error\n");
+               fread(&tab_h, (size_t)sizeof(exp_header_t), (size_t)1, fp );
+               if( tab_h.float_size != sizeof( FLOAT_SIZE ) )
+                       build_exp_table();
+
+               if( !(exp_val = (FLOAT_SIZE *)malloc( tab_h.steps * sizeof( FLOAT_SIZE ) ) ) )
+                       printf("exp_tab: malloc error\n");
+               fread(exp_val, (size_t)sizeof(FLOAT_SIZE), (size_t)tab_h.steps, fp );
+               slope = tab_h.steps / (tab_h.dmax - tab_h.dmin);
+       }
+       if( in < tab_h.dmin || in > tab_h.dmax )
+               return exp( in );
+       else
+               return exp_val[(int)((in-tab_h.dmin) * slope)];
+#endif
+}
+#endif
+
+void f2gauss5(double x,double a[],double *y,double dyda[],int na)
+{
+  int i,index;
+  double fac,fac1,fac2,ex,ex1,ex2,arg1,arg2,u,v;
+  
+       
+       index = nint(x);
+       if( index < 0 || index >=FIT_PTS )
+       {
+               fprintf( stderr, "ff2gauss: wrong index %ld\n", index );
+               return;
+       }
+       u     = plane[index].u;
+       v     = plane[index].v;
+       
+       *y=0.0;
+       for (i=1;i<=na-1;i+=5)
+       {
+               arg1      = (u-a[i+1])/a[i+2];
+               arg2      = (v-a[i+3])/a[i+4];
+               ex1       = exp_tab(-arg1*arg1);
+               ex2       = exp_tab(-arg2*arg2);
+               ex        = ex1*ex2;
+               fac       = a[i]*ex*2.0;
+               fac1      = fac * arg1;
+               fac2      = fac * arg2;
+               *y       += a[i] * ex;
+               dyda[i]   = ex;
+               dyda[i+1] = fac1/a[i+2];
+               dyda[i+2] = fac1*arg1/a[i+2];
+               dyda[i+3] = fac2/a[i+4];
+               dyda[i+4] = fac2*arg2/a[i+4];
+       }
+}
+
+void nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+  printf("%s\n",error_text);
+  exit(1);
+}
+
+void free_vector(double *v, long nl, long nh)
+/* free a double vector allocated with vector() */
+{
+       free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+       free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by matrix() */
+{
+       free((FREE_ARG) (m[nrl]+ncl-NR_END));
+       free((FREE_ARG) (m+nrl-NR_END));
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+       int *v;
+
+       v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+       if (!v) nrerror("allocation failure in ivector()");
+       return v-nl+NR_END;
+}
+
+double *vector(long nl,long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+       double *v;
+
+       v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
+       if (!v) nrerror("allocation failure in vector()");
+       return v-nl+NR_END;
+}
+
+double **matrix(long nrl,long nrh,long ncl,long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+       long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+       double **m;
+
+       /* allocate pointers to rows */
+       m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
+       if (!m) nrerror("allocation failure 1 in matrix()");
+       m += NR_END;
+       m -= nrl;
+
+       /* allocate rows and set pointers to them */
+       m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
+       if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+       m[nrl] += NR_END;
+       m[nrl] -= ncl;
+
+       for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+       /* return pointer to array of pointers to rows */
+       return m;
+}
+
+void gaussj(double **a, int n, double **b, int m)
+{
+       int *indxc,*indxr,*ipiv;
+       int i,icol,irow,j,k,l,ll;
+       double big,dum,pivinv,swap;
+
+       indxc=ivector(1,n);
+       indxr=ivector(1,n);
+       ipiv=ivector(1,n);
+       for (j=1;j<=n;j++) ipiv[j]=0;
+       for (i=1;i<=n;i++) {
+               big=0.0;
+               for (j=1;j<=n;j++)
+                       if (ipiv[j] != 1)
+                               for (k=1;k<=n;k++) {
+                                       if (ipiv[k] == 0) {
+                                               if (fabs(a[j][k]) >= big) {
+                                                       big=fabs(a[j][k]);
+                                                       irow=j;
+                                                       icol=k;
+                                               }
+                                       } else if (ipiv[k] > 1) {
+                                               free_ivector(ipiv,1,n);
+                                               free_ivector(indxr,1,n);
+                                               free_ivector(indxc,1,n);
+                                               nrerror("gaussj: Singular Matrix-1");
+                                       }
+                               }
+               ++(ipiv[icol]);
+               if (irow != icol) {
+                       for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
+                       for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
+               }
+               indxr[i]=irow;
+               indxc[i]=icol;
+               if (a[icol][icol] == 0.0){
+                       free_ivector(ipiv,1,n);
+                       free_ivector(indxr,1,n);
+                       free_ivector(indxc,1,n);
+                       nrerror("gaussj: Singular Matrix-2");
+               }
+               if (mabs(a[icol][icol]) < EPSILON ) 
+                       nrerror("gaussj: a[icol][icol] == 0");
+               pivinv=1.0/a[icol][icol];
+
+               pivinv=1.0/a[icol][icol];
+               a[icol][icol]=1.0;
+               for (l=1;l<=n;l++) a[icol][l] *= pivinv;
+               for (l=1;l<=m;l++) b[icol][l] *= pivinv;
+               for (ll=1;ll<=n;ll++)
+                       if (ll != icol) {
+                               dum=a[ll][icol];
+                               a[ll][icol]=0.0;
+                               for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
+                               for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
+                       }
+       }
+       for (l=n;l>=1;l--) {
+               if (indxr[l] != indxc[l])
+                       for (k=1;k<=n;k++)
+                               SWAP(a[k][indxr[l]],a[k][indxc[l]]);
+       }
+       free_ivector(ipiv,1,n);
+       free_ivector(indxr,1,n);
+       free_ivector(indxc,1,n);
+}
+
+void covsrt(double **covar, int ma, int ia[], int mfit)
+{
+       int i,j,k;
+       double swap;
+
+       for (i=mfit+1;i<=ma;i++)
+               for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
+       k=mfit;
+       for (j=ma;j>=1;j--) {
+               if (ia[j]) {
+                       for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
+                       for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
+                       k--;
+               }
+       }
+}
+
+void mrqcof(double x[], double y[], double sig[], int ndata, double a[], int ia[],
+       int ma, double **alpha, double beta[], double *chisq,
+       void (*funcs)(double, double [], double *, double [], int))
+{
+       int i,j,k,l,m,mfit=0;
+       double ymod,wt,sig2i,dy,*dyda;
+
+       dyda=vector(1,ma);
+
+       for (j=1;j<=ma;j++)
+               if (ia[j]) mfit++;
+       for (j=1;j<=mfit;j++) {
+               for (k=1;k<=j;k++) alpha[j][k]=0.0;
+               beta[j]=0.0;
+       }
+       *chisq=0.0;
+       for (i=1;i<=ndata;i++) {
+               (*funcs)(x[i],a,&ymod,dyda,ma);
+               sig2i=1.0/(sig[i]*sig[i]);
+               dy=y[i]-ymod;
+               for (j=0,l=1;l<=ma;l++) {
+                       if (ia[l]) {
+                               wt=dyda[l]*sig2i;
+                               for (j++,k=0,m=1;m<=l;m++)
+                                       if (ia[m]) alpha[j][++k] += wt*dyda[m];
+                               beta[j] += dy*wt;
+                       }
+               }
+               *chisq += dy*dy*sig2i;
+       }
+       for (j=2;j<=mfit;j++)
+               for (k=1;k<j;k++) alpha[k][j]=alpha[j][k];
+       free_vector(dyda,1,ma);
+}
+
+void mrqmin(double x[], double y[], double sig[], int ndata, double a[], int ia[],
+       int ma, double **covar, double **alpha, double *chisq,
+       void (*funcs)(double, double [], double *, double [], int), double *alamda)
+{
+       void covsrt(double **covar, int ma, int ia[], int mfit);
+       void gaussj(double **a, int n, double **b, int m);
+       void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
+               int ia[], int ma, double **alpha, double beta[], double *chisq,
+               void (*funcs)(double, double [], double *, double [], int));
+       int j,k,l,m;
+       static int mfit;
+       static double ochisq,*atry,*beta,*da,**oneda;
+
+       if (*alamda < 0.0) {
+               atry=vector(1,ma);
+               beta=vector(1,ma);
+               da=vector(1,ma);
+               for (mfit=0,j=1;j<=ma;j++)
+                       if (ia[j]) mfit++;
+               oneda=matrix(1,mfit,1,1);
+               *alamda=0.001;
+               mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,chisq,funcs);
+               ochisq=(*chisq);
+               for (j=1;j<=ma;j++) atry[j]=a[j];
+       }
+       for (j=0,l=1;l<=ma;l++) {
+               if (ia[l]) {
+                       for (j++,k=0,m=1;m<=ma;m++) {
+                               if (ia[m]) {
+                                       k++;
+                                       covar[j][k]=alpha[j][k];
+                               }
+                       }
+                       covar[j][j]=alpha[j][j]*(1.0+(*alamda));
+                       oneda[j][1]=beta[j];
+               }
+       }
+       gaussj(covar,mfit,oneda,1);
+       for (j=1;j<=mfit;j++) da[j]=oneda[j][1];
+       if (*alamda == 0.0) {
+               covsrt(covar,ma,ia,mfit);
+               free_matrix(oneda,1,mfit,1,1);
+               free_vector(da,1,ma);
+               free_vector(beta,1,ma);
+               free_vector(atry,1,ma);
+               return;
+       }
+       for (j=0,l=1;l<=ma;l++)
+               if (ia[l]) atry[l]=a[l]+da[++j];
+       mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq,funcs);
+       if (*chisq < ochisq) {
+               *alamda *= 0.1;
+               ochisq=(*chisq);
+               for (j=0,l=1;l<=ma;l++) {
+                       if (ia[l]) {
+                               for (j++,k=0,m=1;m<=ma;m++) {
+                                       if (ia[m]) {
+                                               k++;
+                                               alpha[j][k]=covar[j][k];
+                                       }
+                               }
+                               beta[j]=da[j];
+                               a[l]=atry[l];
+                       }
+               }
+       } else {
+               *alamda *= 10.0;
+               *chisq=ochisq;
+       }
+}
+
+
+int lev_marq_fit( double x[], double y[], double sig[], int NPT, double a[], int ia[], double dev[], int MA,
+                  double *chisq_p, void (*funcs)(double, double [], double *, double [], int) )
+{
+       int     i,itst,k;
+       double   alamda,chisq,ochisq,**covar,**alpha;
+       
+       covar=matrix(1,MA,1,MA);
+       alpha=matrix(1,MA,1,MA);
+
+       if( setjmp(env) == 0 ) {        
+               alamda = -1;
+               mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
+               k=1;
+               itst=0;
+               for (;;) {
+                       k++;
+                       ochisq=chisq;
+                       mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
+                       if (chisq > ochisq)
+                               itst=0;
+                       else if (fabs(ochisq-chisq) < 0.1)
+                               itst++;
+                       if (itst < 4) continue;
+                       alamda=0.0;
+                       mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
+                       *chisq_p = chisq;
+                       for (i=1;i<=MA;i++) 
+                               dev[i] = sqrt(covar[i][i]);
+                       break;
+               }
+               free_matrix(alpha,1,MA,1,MA);
+               free_matrix(covar,1,MA,1,MA);
+               return 0;
+       }
+       else {
+         //if( control_g.print_fit_errors==2 )
+         fprintf( stderr, " runtime error\n" );
+
+               free_matrix(alpha,1,MA,1,MA);
+               free_matrix(covar,1,MA,1,MA);
+               return -1;
+       }
+}
+
diff --git a/HLT/comp/AliL3FitUtilities.h b/HLT/comp/AliL3FitUtilities.h
new file mode 100644 (file)
index 0000000..4db6851
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef AliL3FitUtilities
+#define AliL3FitUtilities
+
+void   f2gauss5( double, double *, double *,double *,int );
+int lev_marq_fit( double x[], double y[], double sig[], int NPT, double a[], int ia[], double dev[], int MA,
+                  double *chisq_p, void (*funcs)(double, double [], double *, double [], int) );
+
+#define max(x, y) (((x) > (y)) ? (x) : (y))
+#define min(x, y) (((x) < (y)) ? (x) : (y))
+#define mabs(x)   (((x) > 0) ? (x) : (-(x)))
+#define nint(x)   ((int)((x) < 0 ? (x)-0.5 : (x)+0.5))
+#define DBL(x)    ((double)(x))
+#define veclen2(x,y)  (DBL(x)*DBL(x) + DBL(y)*DBL(y))
+#define samesign(x,y)   ((((x)>=0 && (y)>=0) || ((x)<0&&(y)<0)) ? TRUE : FALSE )
+
+#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
+#define SQR(x)       ((x)*(x))
+
+#define NR_END 1
+#define FREE_ARG char*
+#define EPSILON             1.0E-12
+#define TRUE 1
+#define FALSE 1
+#define FIT_PTS     2000
+#define  FIT_MAXPAR   41
+#define NUM_PARS 5
+
+/*--- fitting 2-dimensional cluster --------------------------*/
+struct DPOINT {
+  double u;
+  double v;
+};
+typedef struct DPOINT DPOINT;
+
+extern  DPOINT *plane; //!
+
+typedef struct { 
+                                       long   float_size;
+                                       long   steps;
+                                       double dmin;
+                                       double dmax;
+                                       double step_size;
+                               } exp_header_t;
+typedef float FLOAT_SIZE;
+
+#endif