Checking in the seeds of new cluster fitting code.
[u/mrichter/AliRoot.git] / HLT / comp / AliL3ClusterFitter.cxx
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_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow))*sqrt(2);
333           a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow))*sqrt(2);
334           a[n_overlaps*NUM_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_overlaps*NUM_PARS+3] = sqrt(xywidth)*sqrt(2);
358           a[n_overlaps*NUM_PARS+5] = sqrt(zwidth)*sqrt(2);
359           a[n_overlaps*NUM_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_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
477               tr->SetCluster(fCurrentPadRow,a[n_overlaps*NUM_PARS+2],a[n_overlaps*NUM_PARS+4],tot_charge,0,0,pad_num);
478               cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_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