[u/mrichter/AliRoot.git] / HLT / comp / AliL3ClusterFitter.cxx

//$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;
}
Commit	Line	Data
6b8f6f6e	1	//$Id$
	2
	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	4	//*-- Copyright &copy ASV
	5
	6	#include "AliL3StandardIncludes.h"
	7	#include "AliL3ClusterFitter.h"
	8	#include "AliL3FitUtilities.h"
	9	#include "AliL3Transform.h"
	10	#include "AliL3DigitData.h"
	11	#include "AliL3ModelTrack.h"
	12	#include "AliL3TrackArray.h"
	13	#include "AliL3MemHandler.h"
	14
	15	#if GCCVERSION == 3
	16	using namespace std;
	17	#endif
	18
	19	//_____________________________________________________________
	20	//
	21	// AliL3ClusterFitter
	22	//
	23
	24
	25	ClassImp(AliL3ClusterFitter)
	26
	27	AliL3ClusterFitter::AliL3ClusterFitter()
	28	{
	29	fPadFitRange=2;
	30	fTimeFitRange=3;
	31	plane=0;
	32	fNmaxOverlaps = 3;
	33	fChiSqMax=12;
	34	}
	35
	36	AliL3ClusterFitter::~AliL3ClusterFitter()
	37	{
	38
	39	}
	40
	41	void AliL3ClusterFitter::FindClusters()
	42	{
	43	if(!fTracks)
	44	{
	45	cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
	46	return;
	47	}
	48	if(!fRowData)
	49	{
	50	cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
	51	return;
	52	}
	53
	54	AliL3DigitRowData *rowPt = fRowData;
	55	AliL3DigitData *digPt=0;
	56
	57	Int_t pad,time;
	58	Short_t charge;
	59
	60	for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
	61	{
	62	fCurrentPadRow = i;
	63	memset((void)fRow,0,(AliL3Transform::GetNTimeBins()+1)(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
	64	digPt = (AliL3DigitData*)rowPt->fDigitData;
65	//cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
66	for(UInt_t j=0; j<rowPt->fNDigit; j++)
67	{
68	pad = digPt[j].fPad;
69	time = digPt[j].fTime;
70	charge = digPt[j].fCharge;
71	fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
72	fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
73	//cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
74	}
75
76	for(Int_t k=0; k<fTracks->GetNTracks(); k++)
77	{
78	AliL3ModelTrack track = (AliL3ModelTrack)fTracks->GetCheckedTrack(k);
79	if(!track) continue;
80	if(track->GetPadHit(i) < 0 \|\| track->GetTimeHit(i) < 0)
81	{
82	track->SetCluster(i,0,0,0,0,0,0);
83	continue;
84	}
85	FitClusters(track);
86	}
87	FillZeros(rowPt);
88	AliL3MemHandler::UpdateRowPointer(rowPt);
89	}
90	}
91
92	void AliL3ClusterFitter::FitCluster(AliL3ModelTrack *track)
93	{
94	cout<<"Track had no overlaps"<<endl;
95	Int_t minpad = (Int_t)rint(track->GetPadHit(fCurrentPadRow)) - fPadFitRange;
96	Int_t maxpad = (Int_t)rint(track->GetPadHit(fCurrentPadRow)) + fPadFitRange;
97	Int_t mintime = (Int_t)rint(track->GetTimeHit(fCurrentPadRow)) - fTimeFitRange;
98	Int_t maxtime = (Int_t)rint(track->GetTimeHit(fCurrentPadRow)) + fTimeFitRange;
99
100	Int_t size = FIT_PTS;
101
102	if(minpad <= 0)
103	minpad=0;
104	if(mintime<=0)
105	mintime=0;
106	if(maxpad>=AliL3Transform::GetNPads(fCurrentPadRow)-1)
107	maxpad=AliL3Transform::GetNPads(fCurrentPadRow)-1;
108	if(maxtime>=AliL3Transform::GetNTimeBins()-1)
109	maxtime=AliL3Transform::GetNTimeBins()-1;
110
111	if(plane)
112	delete [] plane;
113	plane = new DPOINT[FIT_PTS];
114	memset(plane,0,FIT_PTS*sizeof(DPOINT));
115
116	Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
117
118	//Fill the fit parameters:
119	Double_t a[FIT_MAXPAR];
120	a[2] = track->GetPadHit(fCurrentPadRow);
121	a[4] = track->GetTimeHit(fCurrentPadRow);
122	a[1] = fRow[(AliL3Transform::GetNTimeBins()+1)*((Int_t)rint(a[2])) + (Int_t)rint(a[4])].fCharge;
123	a[3] = sqrt(track->GetParSigmaY2(fCurrentPadRow))*sqrt(2);
124	a[5] = sqrt(track->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
125	a[6] = sqrt(track->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
126
127	if(!a[1])
128	{
129	cout<<"No charge"<<endl;
130	if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
131	track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
132	return;
133	}
134
135	Int_t pad_num=0;
136	Int_t time_num_max=0;
137	Int_t ndata=0;
138	Int_t charge,tot_charge=0;
139
140	//Fill the proper arrays:
141	for(Int_t i=minpad; i<=maxpad; i++)
142	{
143	Int_t max_charge = 0;
144	Int_t time_num=0;
145	for(Int_t j=mintime; j<=maxtime; j++)
146	{
147	charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
148	if(charge > max_charge)
149	{
150	max_charge = charge;
151	time_num++;
152	}
153	if(charge <= 0) continue;
154	//cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
155	tot_charge += charge;
156	ndata++;
157	if(ndata >= size)
158	cerr<<"Too many points"<<endl;
159	plane[ndata].u = (Double_t)i;
160	plane[ndata].v = (Double_t)j;
161	x[ndata]=ndata;
162	y[ndata]=charge;
163	s[ndata]= 1 + sqrt((Double_t)charge);
164	}
165	if(max_charge) //there was charge on this pad
166	pad_num++;
167	if(time_num_max < time_num)
168	time_num_max = time_num;
169	}
170
171	if(pad_num <= 1 \|\| time_num_max <=1) //too few to do fit
172	{
173	//cout<<"To few pads for fit on row "<<fCurrentPadRow<<endl;
174	if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
175	track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
176
177	return;
178	}
179
180	Int_t npars = NUM_PARS;
181	Int_t lista[FIT_MAXPAR];
182	Double_t dev[FIT_MAXPAR],chisq_f;
183	lista[1] = 1;
184	lista[2] = 1;
185	lista[3] = 0;
186	lista[4] = 1;
187	lista[5] = 0;
188	lista[6] = 0;
189
190	//Declare a pointer to the C fitting function
191	void (funcs) ( double, double , double ,double ,int );
192	funcs = f2gauss5;
193
194	//cout<<"Doing fit with parameters "<<a[2]<<" "<<a[4]<<" "<<a[3]<<" "<<a[5]<<" "<<a[1]<<endl;
195
196	Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f,f2gauss5);
197	if(ret)
198	cerr<<"Fit error"<<endl;
199
200	tot_charge = (Int_t)(a[1] * a[3] * a[5]);
201	chisq_f /= (ndata - 3);
202	// cout<<"Chisq per degree of freedom : "<<chisq_f<<endl;
203	//cout<<"pad "<<a[2]<<" time "<<a[4]<<" adc "<<a[1]<<endl;
204	//cout<<"Setting cluster on row "<<fCurrentPadRow<<" pad "<<a[2]<<" time "<<a[4]<<" charge "<<tot_charge<<endl;
205
206	//Make sure that the cluster has not already been set before setting it:
207	if(track->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
208	track->SetCluster(fCurrentPadRow,a[2],a[4],tot_charge,0,0,pad_num);
209	}
210
211	void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track)
212	{
213	//Handle single and overlapping clusters
214
215	if(!track->IsPresent(fCurrentPadRow) && track->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
216	{
217	if(!track->IsPresent(fCurrentPadRow) &&
218	track->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch) + 1)//debug
219	{
220	cerr<<"AliL3ClusterFitter::FitClusters() : Mismatching clustercount "
221	<<track->GetNClusters()<<" "<<fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch) + 1<<endl;
222	exit(5);
223	}
224	return; //This cluster has been set before.
225	}
226
227
228	Int_t minpad,maxpad,mintime,maxtime;
229
230	Int_t size = FIT_PTS;
231	Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
232	if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
233	{
234	cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
235	return;
236	}
237	Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
238
239	//Check if at least one cluster is not already fitted
240	Bool_t all_fitted=kTRUE;
241	Int_t k=-1;
242	while(k < track->GetNOverlaps(fCurrentPadRow))
243	{
244	AliL3ModelTrack *tr=0;
245	if(k==-1)
246	tr = track;
247	else
248	tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
249	k++;
250	if(!tr) continue;
251	if(!tr->IsPresent(fCurrentPadRow))
252	{
253	all_fitted = kFALSE;
254	break;
255	}
256	}
257	if(all_fitted)
258	{
259	cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
260	return;
261	}
262
263	plane = new DPOINT[FIT_PTS];
264	memset(plane,0,FIT_PTS*sizeof(DPOINT));
265
266	Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
267
268	//Fill the fit parameters:
269	Double_t a[FIT_MAXPAR];
270	Int_t lista[FIT_MAXPAR];
271	Double_t dev[FIT_MAXPAR],chisq_f;
272
273	minpad = mintime = 999;
274	maxpad = maxtime = 0;
275	Int_t fit_pars=0;
276
277	Int_t ntracks = track->GetNOverlaps(fCurrentPadRow)+1;
278	Bool_t *do_fit = new Bool_t[ntracks];
279	for(Int_t i=0; i<ntracks; i++)
280	do_fit[i]=kTRUE;
281
282	Int_t n_overlaps=0;
283	k=-1;
284	//Fill the overlapping tracks:
285	while(k < track->GetNOverlaps(fCurrentPadRow))
286	{
287	AliL3ModelTrack *tr=0;
288	if(k==-1)
289	tr = track;
290	else
291	tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
292	k++;
293	if(!tr) continue;
294
295	Int_t hitpad = (Int_t)rint(tr->GetPadHit(fCurrentPadRow));
296	Int_t hittime = (Int_t)rint(tr->GetTimeHit(fCurrentPadRow));
297	Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
298	if(!charge) //There is not charge here, so the cluster is non-existing -- remove it.
299	{
300	if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
301	tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
302	do_fit[k] = kFALSE;
303	continue;
304	}
305
306	if(k==0)
307	LocateCluster(tr,minpad,maxpad,mintime,maxtime);
308
309	if(tr->GetPadHit(fCurrentPadRow) < (Double_t)minpad \|\| tr->GetPadHit(fCurrentPadRow) > (Double_t)maxpad \|\|
310	tr->GetTimeHit(fCurrentPadRow) < (Double_t)mintime \|\| tr->GetTimeHit(fCurrentPadRow) > (Double_t)maxtime)
311	{
312	do_fit[k] = kFALSE;//This cluster is outside the region already specified, so it will not be included in this fit.
313
314	//If this is the first: remove it, because it will not be checked again.
315	if(k==0)
316	if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
317	tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
318	continue;
319	}
320
321
322	cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
323	<<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
324	<<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
325
326	a[n_overlaps*NUM_PARS+2] = tr->GetPadHit(fCurrentPadRow);
327	a[n_overlaps*NUM_PARS+4] = tr->GetTimeHit(fCurrentPadRow);
328
329	if(!tr->IsPresent(fCurrentPadRow)) //Cluster is not fitted before
330	{
331	a[n_overlaps*NUM_PARS+1] = charge;
332	a[n_overlapsNUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow))sqrt(2);
333	a[n_overlapsNUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow))sqrt(2);
334	a[n_overlapsNUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow))sqrt(2);
335	lista[n_overlaps*NUM_PARS + 1] = 1;
336	lista[n_overlaps*NUM_PARS + 2] = 1;
337	lista[n_overlaps*NUM_PARS + 3] = 0;
338	lista[n_overlaps*NUM_PARS + 4] = 1;
339	lista[n_overlaps*NUM_PARS + 5] = 0;
340	lista[n_overlaps*NUM_PARS + 6] = 0;
341	fit_pars += 3;
342	}
343	else
344	{
345	Int_t charge;
346	Float_t xywidth,zwidth,pad,time;
347	tr->GetPad(fCurrentPadRow,pad);
348	tr->GetTime(fCurrentPadRow,time);
349	tr->GetClusterCharge(fCurrentPadRow,charge);
350	xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
351	zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
352	cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
353
354	a[n_overlaps*NUM_PARS+2] = pad;
355	a[n_overlaps*NUM_PARS+4] = time;
356	a[n_overlaps*NUM_PARS+1] = charge;
357	a[n_overlapsNUM_PARS+3] = sqrt(xywidth)sqrt(2);
358	a[n_overlapsNUM_PARS+5] = sqrt(zwidth)sqrt(2);
359	a[n_overlapsNUM_PARS+6] = sqrt(zwidth)sqrt(2);
360
361	lista[n_overlaps*NUM_PARS + 1] = 1;
362	lista[n_overlaps*NUM_PARS + 2] = 0;
363	lista[n_overlaps*NUM_PARS + 3] = 0;
364	lista[n_overlaps*NUM_PARS + 4] = 0;
365	lista[n_overlaps*NUM_PARS + 5] = 0;
366	lista[n_overlaps*NUM_PARS + 6] = 0;
367	fit_pars += 1;
368	}
369	n_overlaps++;
370	}
371
372	if(n_overlaps==0) //No clusters here
373	return;
374	cout<<"Setting init searchrange; pad "<<minpad<<" "<<maxpad<<" time "<<mintime<<" "<<maxtime<<endl;
375
376	if(minpad <= 0)
377	minpad=0;
378	if(mintime<=0)
379	mintime=0;
380	if(maxpad>=AliL3Transform::GetNPads(fCurrentPadRow)-1)
381	maxpad=AliL3Transform::GetNPads(fCurrentPadRow)-1;
382	if(maxtime>=AliL3Transform::GetNTimeBins()-1)
383	maxtime=AliL3Transform::GetNTimeBins()-1;
384
385	Int_t pad_num=0;
386	Int_t time_num_max=0;
387	Int_t ndata=0;
388	Int_t tot_charge=0;
389
390	for(Int_t i=minpad; i<=maxpad; i++)
391	{
392	Int_t max_charge = 0;
393	Int_t time_num=0;
394	for(Int_t j=mintime; j<=maxtime; j++)
395	{
396	Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
397
398	if(charge <= 0) continue;
399
400	if(charge > max_charge)
401	{
402	max_charge = charge;
403	time_num++;
404	}
405	cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
406	tot_charge += charge;
407	ndata++;
408	if(ndata >= size)
409	cerr<<"Too many points"<<endl;
410
411	//This digit will most likely be used:
412	fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
413	plane[ndata].u = (Double_t)i;
414	plane[ndata].v = (Double_t)j;
415	x[ndata]=ndata;
416	y[ndata]=charge;
417	s[ndata]= 1 + sqrt((Double_t)charge);
418	}
419	if(max_charge) //there was charge on this pad
420	pad_num++;
421	if(time_num_max < time_num)
422	time_num_max = time_num;
423	}
424
425	k=-1;
426	if(pad_num <= 1 \|\| time_num_max <=1 \|\| n_overlaps > fNmaxOverlaps) //too few to do fit
427	{
428	while(k < track->GetNOverlaps(fCurrentPadRow))
429	{
430	AliL3ModelTrack *tr=0;
431	if(k==-1)
432	tr = track;
433	else
434	tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
435	k++;
436	if(!tr) continue;
437	if(do_fit[k] == kFALSE) continue;
438
439	if(tr->GetNClusters() == fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch))
440	tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
441	}
442	cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<endl;
443	return;
444	}
445
446	Int_t npars = n_overlaps * NUM_PARS;
447	cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
448	Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
449	if(ret)
450	{
451	cerr<<"Fit error"<<endl;
452	exit(5);
453	}
454
455	chisq_f /= (ndata-fit_pars);
456	cout<<"Chisq "<<chisq_f<<endl;
457
458	k=-1;
459	n_overlaps=0;
460	while(k < track->GetNOverlaps(fCurrentPadRow))
461	{
462	AliL3ModelTrack *tr=0;
463	if(k==-1)
464	tr = track;
465	else
466	tr = (AliL3ModelTrack*)fTracks->GetCheckedTrack(overlaps[k]);
467	k++;
468	if(!tr) continue;
469	if(do_fit[k] == kFALSE) continue;
470
471	if(!tr->IsPresent(fCurrentPadRow))
472	{
473	if(tr->GetNClusters() != fCurrentPadRow - AliL3Transform::GetFirstRow(fPatch)) continue;//This cluster has been set before
474	if(chisq_f < fChiSqMax)//cluster fit is good enough
475	{
476	tot_charge = (Int_t)(a[n_overlapsNUM_PARS+1] a[n_overlapsNUM_PARS+3] a[n_overlaps*NUM_PARS+5]);
477	tr->SetCluster(fCurrentPadRow,a[n_overlapsNUM_PARS+2],a[n_overlapsNUM_PARS+4],tot_charge,0,0,pad_num);
478	cout<<"Setting cluster in pad "<<a[n_overlapsNUM_PARS+2]<<" time "<<a[n_overlapsNUM_PARS+4]<<" charge "<<tot_charge<<endl;
479	}
480	else //fit was too bad
481	tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
482	}
483	n_overlaps++;
484	}
485
486	delete [] plane;
487	delete [] do_fit;
488	}
489