]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/comp/AliL3ClusterFitter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merged HLT tag v1-2 with ALIROOT tag v3-09-Release.
[u/mrichter/AliRoot.git] / HLT / comp / AliL3ClusterFitter.cxx
CommitLineData
3e87ef69 1// @(#) $Id$
6b8f6f6e 2
3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
3e87ef69 4//*-- Copyright &copy ALICE HLT Group
6b8f6f6e 5
6#include "AliL3StandardIncludes.h"
3e87ef69 7
8#include "AliL3Logging.h"
6b8f6f6e 9#include "AliL3ClusterFitter.h"
10#include "AliL3FitUtilities.h"
6b8f6f6e 11#include "AliL3DigitData.h"
12#include "AliL3ModelTrack.h"
13#include "AliL3TrackArray.h"
14#include "AliL3MemHandler.h"
3e87ef69 15#include "AliL3HoughTrack.h"
16#include "AliL3SpacePointData.h"
17#include "AliL3Compress.h"
6b8f6f6e 18
19#if GCCVERSION == 3
20using namespace std;
21#endif
22
23//_____________________________________________________________
24//
25// AliL3ClusterFitter
26//
27
28
29ClassImp(AliL3ClusterFitter)
30
3e87ef69 31Int_t AliL3ClusterFitter::fBadFitError=0;
32Int_t AliL3ClusterFitter::fFitError=0;
33Int_t AliL3ClusterFitter::fResultError=0;
34Int_t AliL3ClusterFitter::fFitRangeError=0;
35
6b8f6f6e 36AliL3ClusterFitter::AliL3ClusterFitter()
37{
6b8f6f6e 38 plane=0;
39 fNmaxOverlaps = 3;
40 fChiSqMax=12;
3e87ef69 41 fRowMin=-1;
42 fRowMax=-1;
43 fFitted=0;
44 fFailed=0;
45 fYInnerWidthFactor=1;
46 fZInnerWidthFactor=1;
47 fYOuterWidthFactor=1;
48 fZOuterWidthFactor=1;
49 fSeeds=0;
50 fProcessTracks=0;
51 fClusters=0;
52 fNMaxClusters=0;
53 fNClusters=0;
54}
55
56AliL3ClusterFitter::AliL3ClusterFitter(Char_t *path)
57{
58 strcpy(fPath,path);
59 plane=0;
60 fNmaxOverlaps = 3;
61 fChiSqMax=12;
62 fRowMin=-1;
63 fRowMax=-1;
64 fFitted=0;
65 fFailed=0;
66 fYInnerWidthFactor=1;
67 fZInnerWidthFactor=1;
68 fYOuterWidthFactor=1;
69 fZOuterWidthFactor=1;
70 fSeeds=0;
71 fProcessTracks=0;
72 fNMaxClusters=100000;
73 fClusters=0;
74 fNClusters=0;
6b8f6f6e 75}
76
77AliL3ClusterFitter::~AliL3ClusterFitter()
78{
3e87ef69 79 if(fSeeds)
80 delete fSeeds;
81 if(fClusters)
82 delete [] fClusters;
83}
84
85void AliL3ClusterFitter::Init(Int_t slice,Int_t patch,Int_t *rowrange,AliL3TrackArray *tracks)
86{
87 //Assuming tracklets found by the line transform
88
89 fSlice=slice;
90 fPatch=patch;
91
92 if(rowrange[0] > AliL3Transform::GetLastRow(patch) || rowrange[1] < AliL3Transform::GetFirstRow(patch))
93 cerr<<"AliL3ClusterFitter::Init : Wrong rows "<<rowrange[0]<<" "<<rowrange[1]<<endl;
94 fRowMin=rowrange[0];
95 fRowMax=rowrange[1];
96
97 if(fRowMin < 0)
98 fRowMin = 0;
99 if(fRowMax > AliL3Transform::GetLastRow(fPatch))
100 fRowMax = AliL3Transform::GetLastRow(fPatch);
101
102 fFitted=fFailed=0;
103
104 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
105 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
106 Int_t bounds = ntimes*npads;
107 if(fRow)
108 delete [] fRow;
109 fRow = new Digit[bounds];
110 if(fTracks)
111 delete fTracks;
112
113 fTracks = new AliL3TrackArray("AliL3ModelTrack");
114
115 for(Int_t i=0; i<tracks->GetNTracks(); i++)
116 {
117 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
118 if(!track) continue;
119 AliL3ModelTrack *mtrack = (AliL3ModelTrack*)fTracks->NextTrack();
120 mtrack->Init(slice,patch);
121 mtrack->SetTgl(track->GetTgl());
122 mtrack->SetRowRange(rowrange[0],rowrange[1]);
123 for(Int_t j=fRowMin; j<=fRowMax; j++)
124 {
125 Float_t hit[3];
126 track->GetLineCrossingPoint(j,hit);
127 hit[0] += AliL3Transform::Row2X(track->GetFirstRow());
128 Float_t R = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
129 hit[2] = R*track->GetTgl();
130 Int_t se,ro;
131 AliL3Transform::Slice2Sector(slice,j,se,ro);
132 AliL3Transform::Local2Raw(hit,se,ro);
133 if(hit[1]<0 || hit[1]>=AliL3Transform::GetNPads(j) || hit[2]<0 || hit[2]>=AliL3Transform::GetNTimeBins())
134 {
135 mtrack->SetPadHit(j,-1);
136 mtrack->SetTimeHit(j,-1);
137 continue;
138 }
139 mtrack->SetPadHit(j,hit[1]);
140 mtrack->SetTimeHit(j,hit[2]);
141 mtrack->SetCrossingAngleLUT(j,fabs(track->GetPsiLine() - AliL3Transform::Pi()/2));
142 //if(mtrack->GetCrossingAngleLUT(j) > AliL3Transform::Deg2Rad(20))
143 // cout<<"Angle "<<mtrack->GetCrossingAngleLUT(j)<<" psiline "<<track->GetPsiLine()*180/3.1415<<endl;
144 mtrack->CalculateClusterWidths(j);
145 }
146 }
147 // cout<<"Copied "<<fTracks->GetNTracks()<<" tracks "<<endl;
148}
149
150void AliL3ClusterFitter::Init(Int_t slice,Int_t patch)
151{
152 fSlice=slice;
153 fPatch=patch;
154
155 fRowMin=AliL3Transform::GetFirstRow(patch);
156 fRowMax=AliL3Transform::GetLastRow(patch);
157
158 fFitted=fFailed=0;
159
160 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
161 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
162 Int_t bounds = ntimes*npads;
163 if(fRow)
164 delete [] fRow;
165 fRow = new Digit[bounds];
166 if(fTracks)
167 delete fTracks;
168 fTracks = new AliL3TrackArray("AliL3ModelTrack");
169
170}
171
6b8f6f6e 172
3e87ef69 173void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline)
174{
175
176 cout<<"Loading the seeds"<<endl;
177 Char_t fname[1024];
178
179 if(offline)
180 sprintf(fname,"%s/offline/tracks_%d.raw",fPath,0);
181 else
182 sprintf(fname,"%s/hough/tracks_%d.raw",fPath,0);
183
184 cout<<"AliL3ClusterFitter::LoadSeeds : Loading input tracks from "<<fname<<endl;
185
186 AliL3MemHandler tfile;
187 tfile.SetBinaryInput(fname);
188
189 if(fSeeds)
190 delete fSeeds;
191 fSeeds = new AliL3TrackArray("AliL3ModelTrack");
192 tfile.Binary2TrackArray(fSeeds);
193 tfile.CloseBinaryInput();
194
195 //if(!offline)
196 //fSeeds->QSort();
197
198 Int_t clustercount=0;
199 for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
200 {
201 AliL3ModelTrack *track = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
202 if(!track) continue;
203
204 if(!offline)
205 {
206 /*
207 if(track->GetNHits() < 10 || track->GetPt() < 0.08)
208 {
209 fSeeds->Remove(i);
210 continue;
211 }
212 */
213 }
214 clustercount += track->GetNHits();
215 track->CalculateHelix();
216
217 Int_t nhits = track->GetNHits();
218 UInt_t *hitids = track->GetHitNumbers();
219
220 Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
221
222 track->Init(origslice,-1);
223 Int_t slice = origslice;
224
225 //for(Int_t j=rowrange[1]; j>=rowrange[0]; j--)
226 for(Int_t j=rowrange[0]; j<=rowrange[1]; j++)
227 {
228
229 //Calculate the crossing point between track and padrow
230 Float_t angle = 0; //Perpendicular to padrow in local coordinates
231 AliL3Transform::Local2GlobalAngle(&angle,slice);
232 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
233 {
234 cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
235 continue;
236 //track->Print();
237 //exit(5);
238 }
239 Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
240
241 Int_t sector,row;
242 AliL3Transform::Slice2Sector(slice,j,sector,row);
243 AliL3Transform::Global2Raw(xyz_cross,sector,row);
244 //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
245 if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
246 {
247 newslice:
248
249 Int_t tslice=slice;
250 Float_t lastcross=xyz_cross[1];
251 if(xyz_cross[1] > 0)
252 {
253 if(slice == 17)
254 slice=0;
255 else if(slice == 35)
256 slice = 18;
257 else
258 slice += 1;
259 }
260 else
261 {
262 if(slice == 0)
263 slice = 17;
264 else if(slice==18)
265 slice = 35;
266 else
267 slice -= 1;
268 }
269 if(slice < 0 || slice>35)
270 {
271 cerr<<"Wrong slice "<<slice<<" on row "<<j<<endl;
272 exit(5);
273 }
274 //cout<<"Track leaving, trying slice "<<slice<<endl;
275 angle=0;
276 AliL3Transform::Local2GlobalAngle(&angle,slice);
277 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
278 {
279 cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
280 continue;
281 //track->Print();
282 //exit(5);
283 }
284 xyz_cross[0] = track->GetPointX();
285 xyz_cross[1] = track->GetPointY();
286 xyz_cross[2] = track->GetPointZ();
287 Int_t sector,row;
288 AliL3Transform::Slice2Sector(slice,j,sector,row);
289 AliL3Transform::Global2Raw(xyz_cross,sector,row);
290 if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
291 {
292 if(xyz_cross[1] > 0 && lastcross > 0 || xyz_cross[1] < 0 && lastcross < 0)
293 goto newslice;
294 else
295 {
296 slice = tslice;//Track is on the border of two slices
297 continue;
298 }
299 }
300 }
301
302 if(xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
303 continue;
304
305 if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j))
306 {
307 cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
308 track->Print();
309 exit(5);
310 }
311
312 track->SetPadHit(j,xyz_cross[1]);
313 track->SetTimeHit(j,xyz_cross[2]);
314 angle=0;
315 AliL3Transform::Local2GlobalAngle(&angle,slice);
316 Float_t crossingangle = track->GetCrossingAngle(j,slice);
317 track->SetCrossingAngleLUT(j,crossingangle);
318
319 track->CalculateClusterWidths(j);
320
321 track->GetClusterModel(j)->fSlice = slice;
322
323 }
324 memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));
325 track->SetNHits(0);
326 }
327 fSeeds->Compress();
328
329 AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
330 c->WriteFile(fSeeds,"tracks_before.raw");
331 delete c;
332
333 cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
6b8f6f6e 334}
335
336void AliL3ClusterFitter::FindClusters()
337{
338 if(!fTracks)
339 {
340 cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
341 return;
342 }
343 if(!fRowData)
344 {
345 cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
346 return;
347 }
348
349 AliL3DigitRowData *rowPt = fRowData;
350 AliL3DigitData *digPt=0;
351
352 Int_t pad,time;
353 Short_t charge;
3e87ef69 354
355 if(fRowMin < 0)
356 {
357 fRowMin = AliL3Transform::GetFirstRow(fPatch);
358 fRowMax = AliL3Transform::GetLastRow(fPatch);
359 }
6b8f6f6e 360 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
361 {
3e87ef69 362 if((Int_t)rowPt->fRow < fRowMin)
363 {
364 AliL3MemHandler::UpdateRowPointer(rowPt);
365 continue;
366 }
367 else if((Int_t)rowPt->fRow > fRowMax)
368 break;
369 else if((Int_t)rowPt->fRow != i)
370 {
371 cerr<<"AliL3ClusterFitter::FindClusters : Mismatching row numbering "<<i<<" "<<rowPt->fRow<<endl;
372 exit(5);
373 }
6b8f6f6e 374 fCurrentPadRow = i;
375 memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
376 digPt = (AliL3DigitData*)rowPt->fDigitData;
3e87ef69 377
6b8f6f6e 378 for(UInt_t j=0; j<rowPt->fNDigit; j++)
379 {
380 pad = digPt[j].fPad;
381 time = digPt[j].fTime;
382 charge = digPt[j].fCharge;
3e87ef69 383 if(charge > 1024)
384 charge -= 1024;
6b8f6f6e 385 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
386 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
387 //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
388 }
389
3e87ef69 390 for(Int_t it=0; it<2; it++)
6b8f6f6e 391 {
3e87ef69 392 if(it==0)
6b8f6f6e 393 {
3e87ef69 394 fProcessTracks = fSeeds;
395 fSeeding = kTRUE;
396 }
397 else
398 {
399 fProcessTracks = fTracks;
400 fSeeding = kFALSE;
401 }
402 if(!fProcessTracks)
403 continue;
404
405 for(Int_t k=0; k<fProcessTracks->GetNTracks(); k++)
406 {
407 AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(k);
408 if(!track) continue;
409
410 if(fSeeding)
411 if(track->GetClusterModel(i)->fSlice != fSlice) continue;
412
413 if(track->GetPadHit(i) < 0 || track->GetPadHit(i) > AliL3Transform::GetNPads(i)-1 ||
414 track->GetTimeHit(i) < 0 || track->GetTimeHit(i) > AliL3Transform::GetNTimeBins()-1)
415 {
416 track->SetCluster(i,0,0,0,0,0,0);
417 continue;
418 }
419
420 if(CheckCluster(k) == kFALSE)
421 fFailed++;
6b8f6f6e 422 }
6b8f6f6e 423 }
3e87ef69 424 //FillZeros(rowPt);
6b8f6f6e 425 AliL3MemHandler::UpdateRowPointer(rowPt);
426 }
3e87ef69 427
428 fSeeding = kTRUE;
429 AddClusters();
430 fSeeding = kFALSE;
431 AddClusters();
432
433 cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
434 cout<<"Distribution:"<<endl;
435 cout<<"Bad fit "<<fBadFitError<<endl;
436 cout<<"Fit error "<<fFitError<<endl;
437 cout<<"Result error "<<fResultError<<endl;
438 cout<<"Fit range error "<<fFitRangeError<<endl;
439
6b8f6f6e 440}
441
3e87ef69 442Bool_t AliL3ClusterFitter::CheckCluster(Int_t trackindex)
6b8f6f6e 443{
3e87ef69 444 //Check if this is a single or overlapping cluster
6b8f6f6e 445
3e87ef69 446 AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(trackindex);
6b8f6f6e 447
3e87ef69 448 Int_t row = fCurrentPadRow;
449
450 if(track->IsSet(row)) //A fit has already be performed on this one
451 return kTRUE;
452
453 //Define the cluster region of this hit:
454 Int_t padr[2]={999,-1};
455 Int_t timer[2]={999,-1};
456
457 if(!SetFitRange(track,padr,timer))
6b8f6f6e 458 {
3e87ef69 459 track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
460
461 if(fDebug)
462 cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
463 <<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
464 <<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
465 fFitRangeError++;
466 return kFALSE;
6b8f6f6e 467 }
3e87ef69 468
469 //Check if any other track contributes to this cluster:
470 //This is done by checking if the tracks are overlapping within
471 //the range defined by the track parameters
472 for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
473 {
474 AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
475 if(!tr) continue;
476 if(fSeeding)
477 if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
478 Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))) + 1;
479 Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))) + 1;
480 if(
481 (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
482 tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
483
484 (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
485 tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
486
487 (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
488 tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) ||
489
490 (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
491 tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1])
492 )
493 {
494 if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
495 track->SetOverlap(row,t); //Set overlap
496 }
497 }
498
499 if(fDebug)
500 cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
501 FitClusters(track,padr,timer);
502 return kTRUE;
503}
504
505Bool_t AliL3ClusterFitter::SetFitRange(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
506{
507 Int_t row = fCurrentPadRow;
508 Int_t nt = AliL3Transform::GetNTimeBins()+1;
6b8f6f6e 509
3e87ef69 510 Int_t nsearchbins=0;
511 if(row < 63)
512 nsearchbins=25;
513 else
514 nsearchbins=49;
515 Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
516 ,0,1,2,3,3,3,3,3,3,3
517 ,2,1,0,-1,-2,-3
518 ,-3,-3,-3,-3,-3,-3,-1,-1};
519 Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
520 ,-3,-3,-3,-3,-2,-1,0,1,2,3
521 ,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
522
523 Int_t padhit = (Int_t)rint(track->GetPadHit(row));
524 Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
525 Int_t padmax=-1;
526 Int_t timemax=-1;
6b8f6f6e 527
3e87ef69 528 for(Int_t index=0; index<nsearchbins; index++)
6b8f6f6e 529 {
3e87ef69 530 if(IsMaximum(padhit + padloop[index],timehit + timeloop[index]))
6b8f6f6e 531 {
3e87ef69 532 padmax = padhit + padloop[index];
533 timemax = timehit + timeloop[index];
534 break;
6b8f6f6e 535 }
6b8f6f6e 536 }
537
3e87ef69 538 /*
539 if(IsMaximum(padhit,timehit))
540 {
541 padmax = padhit;
542 timemax = timehit;
543 }
544 */
545 //Define the cluster region of this hit:
546 //The region we look for, is centered at the local maxima
547 //and expanded around using the parametrized cluster width
548 //according to track parameters.
549 Int_t xyw = (Int_t)ceil(sqrt(track->GetParSigmaY2(row)))+1;
550 Int_t zw = (Int_t)ceil(sqrt(track->GetParSigmaZ2(row)))+1;
551 if(padmax>=0 && timemax>=0)
6b8f6f6e 552 {
3e87ef69 553 if(fDebug)
554 {
555 cout<<"Expanding cluster range using expected cluster widths: "<<xyw<<" "<<zw
556 <<" and setting local maxima pad "<<padmax<<" time "<<timemax<<endl;
557 if(xyw > 10 || zw > 10)
558 track->Print();
559 }
560
561 //Set the hit to the local maxima of the cluster.
562 //Store the maxima in the cluster model structure,
563 //-only temporary, it will be overwritten when calling SetCluster.
564
565 track->GetClusterModel(row)->fDPad = padmax;
566 track->GetClusterModel(row)->fDTime = timemax;
6b8f6f6e 567
3e87ef69 568 for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
569 {
570 for(Int_t j=timemax-zw; j<=timemax+zw; j++)
571 {
572 if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
573 if(fRow[nt*i+j].fCharge)
574 {
575 if(i < padrange[0]) padrange[0]=i;
576 if(i > padrange[1]) padrange[1]=i;
577 if(j < timerange[0]) timerange[0]=j;
578 if(j > timerange[1]) timerange[1]=j;
579 }
580 }
581 }
582 if(fDebug)
583 cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
584 return kTRUE;
6b8f6f6e 585 }
3e87ef69 586 return kFALSE;
587}
588
589Bool_t AliL3ClusterFitter::IsMaximum(Int_t pad,Int_t time)
590{
591 if(pad<0 || pad >= AliL3Transform::GetNPads(fCurrentPadRow) ||
592 time<0 || time >= AliL3Transform::GetNTimeBins())
593 return kFALSE;
594 Int_t nt = AliL3Transform::GetNTimeBins()+1;
595 if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
596 Int_t charge = fRow[nt*pad+time].fCharge;
597 if(charge == 1023 || charge==0) return kFALSE;
6b8f6f6e 598
3e87ef69 599 fRow[nt*pad+time].fUsed = kTRUE;
600 return kTRUE;
601
602 //if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
603 if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
604 //if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
605 if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
606 if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
607 //if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
608 if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
609 //if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
610 fRow[nt*pad+time].fUsed = kTRUE;
611 return kTRUE;
6b8f6f6e 612}
613
3e87ef69 614void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
6b8f6f6e 615{
616 //Handle single and overlapping clusters
617
3e87ef69 618 //Check whether this cluster has been set before:
6b8f6f6e 619
620 Int_t size = FIT_PTS;
621 Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
622 if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
623 {
624 cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
625 return;
626 }
627 Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
628
629 //Check if at least one cluster is not already fitted
630 Bool_t all_fitted=kTRUE;
3e87ef69 631
6b8f6f6e 632 Int_t k=-1;
633 while(k < track->GetNOverlaps(fCurrentPadRow))
634 {
635 AliL3ModelTrack *tr=0;
636 if(k==-1)
3e87ef69 637 tr = track;
6b8f6f6e 638 else
3e87ef69 639 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
6b8f6f6e 640 k++;
641 if(!tr) continue;
3e87ef69 642 if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
6b8f6f6e 643 {
644 all_fitted = kFALSE;
645 break;
646 }
647 }
648 if(all_fitted)
649 {
3e87ef69 650 if(fDebug)
651 cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
6b8f6f6e 652 return;
653 }
654
3e87ef69 655 //Allocate fit parameters array; this is interface to the C code
6b8f6f6e 656 plane = new DPOINT[FIT_PTS];
657 memset(plane,0,FIT_PTS*sizeof(DPOINT));
658
659 Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
660
661 //Fill the fit parameters:
662 Double_t a[FIT_MAXPAR];
663 Int_t lista[FIT_MAXPAR];
664 Double_t dev[FIT_MAXPAR],chisq_f;
665
6b8f6f6e 666 Int_t fit_pars=0;
667
6b8f6f6e 668 Int_t n_overlaps=0;
669 k=-1;
3e87ef69 670
6b8f6f6e 671 //Fill the overlapping tracks:
672 while(k < track->GetNOverlaps(fCurrentPadRow))
673 {
674 AliL3ModelTrack *tr=0;
675 if(k==-1)
676 tr = track;
677 else
3e87ef69 678 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
6b8f6f6e 679 k++;
680 if(!tr) continue;
681
3e87ef69 682 if(tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow)) continue;//Cluster fit failed before
683
684 //Use the local maxima as the input to the fitting routine.
685 //The local maxima is temporary stored in the cluster model:
686 Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad); //rint(tr->GetPadHit(fCurrentPadRow));
687 Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime); //rint(tr->GetTimeHit(fCurrentPadRow));
6b8f6f6e 688 Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
6b8f6f6e 689
3e87ef69 690 if(fDebug)
691 cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
692 <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" at local maxima in pad "<<hitpad
693 <<" time "<<hittime<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
694 <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
695
696 if(charge==0)
6b8f6f6e 697 {
3e87ef69 698 cerr<<"Charge still zero!"<<endl;
699 exit(5);
6b8f6f6e 700 }
3e87ef69 701
702 a[n_overlaps*NUM_PARS+2] = hitpad;
703 a[n_overlaps*NUM_PARS+4] = hittime;
6b8f6f6e 704
3e87ef69 705 if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
6b8f6f6e 706 {
707 a[n_overlaps*NUM_PARS+1] = charge;
3e87ef69 708 a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
709 a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
710 a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
6b8f6f6e 711 lista[n_overlaps*NUM_PARS + 1] = 1;
712 lista[n_overlaps*NUM_PARS + 2] = 1;
713 lista[n_overlaps*NUM_PARS + 3] = 0;
714 lista[n_overlaps*NUM_PARS + 4] = 1;
715 lista[n_overlaps*NUM_PARS + 5] = 0;
716 lista[n_overlaps*NUM_PARS + 6] = 0;
717 fit_pars += 3;
718 }
3e87ef69 719 else //Cluster was fitted before
6b8f6f6e 720 {
3e87ef69 721 if(!tr->IsPresent(fCurrentPadRow))
722 {
723 cerr<<"AliL3ClusterFitter::FindClusters : Cluster not present; there is a bug here"<<endl;
724 exit(5);
725 }
6b8f6f6e 726 Int_t charge;
727 Float_t xywidth,zwidth,pad,time;
728 tr->GetPad(fCurrentPadRow,pad);
729 tr->GetTime(fCurrentPadRow,time);
730 tr->GetClusterCharge(fCurrentPadRow,charge);
731 xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
732 zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
3e87ef69 733 if(fDebug)
734 cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
6b8f6f6e 735
736 a[n_overlaps*NUM_PARS+2] = pad;
737 a[n_overlaps*NUM_PARS+4] = time;
738 a[n_overlaps*NUM_PARS+1] = charge;
3e87ef69 739 a[n_overlaps*NUM_PARS+3] = sqrt(xywidth) * GetYWidthFactor();
740 a[n_overlaps*NUM_PARS+5] = sqrt(zwidth) * GetZWidthFactor();
741 a[n_overlaps*NUM_PARS+6] = sqrt(zwidth) * GetZWidthFactor();
6b8f6f6e 742
743 lista[n_overlaps*NUM_PARS + 1] = 1;
744 lista[n_overlaps*NUM_PARS + 2] = 0;
745 lista[n_overlaps*NUM_PARS + 3] = 0;
746 lista[n_overlaps*NUM_PARS + 4] = 0;
747 lista[n_overlaps*NUM_PARS + 5] = 0;
748 lista[n_overlaps*NUM_PARS + 6] = 0;
749 fit_pars += 1;
750 }
751 n_overlaps++;
752 }
753
754 if(n_overlaps==0) //No clusters here
3e87ef69 755 {
756 delete [] plane;
757 return;
758 }
6b8f6f6e 759
6b8f6f6e 760 Int_t pad_num=0;
761 Int_t time_num_max=0;
762 Int_t ndata=0;
763 Int_t tot_charge=0;
3e87ef69 764 if(fDebug)
765 cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
766 for(Int_t i=padrange[0]; i<=padrange[1]; i++)
6b8f6f6e 767 {
768 Int_t max_charge = 0;
769 Int_t time_num=0;
3e87ef69 770 for(Int_t j=timerange[0]; j<=timerange[1]; j++)
6b8f6f6e 771 {
772 Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
773
774 if(charge <= 0) continue;
775
3e87ef69 776 time_num++;
6b8f6f6e 777 if(charge > max_charge)
778 {
779 max_charge = charge;
3e87ef69 780 //time_num++;
6b8f6f6e 781 }
3e87ef69 782 if(fDebug)
783 cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
6b8f6f6e 784 tot_charge += charge;
785 ndata++;
786 if(ndata >= size)
3e87ef69 787 {
788 cerr<<"Too many points; row "<<fCurrentPadRow<<" padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "
789 <<timerange[0]<<" "<<timerange[1]<<endl;
790 exit(5);
791 }
792
6b8f6f6e 793 plane[ndata].u = (Double_t)i;
794 plane[ndata].v = (Double_t)j;
795 x[ndata]=ndata;
796 y[ndata]=charge;
797 s[ndata]= 1 + sqrt((Double_t)charge);
798 }
799 if(max_charge) //there was charge on this pad
800 pad_num++;
801 if(time_num_max < time_num)
802 time_num_max = time_num;
803 }
804
3e87ef69 805 if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps || ndata <= fit_pars) //too few to do fit
6b8f6f6e 806 {
3e87ef69 807 SetClusterfitFalse(track);
808 if(fDebug)
809 cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<" ndata "<<ndata<<" fit_pars "<<fit_pars<<endl;
810 delete [] plane;
6b8f6f6e 811 return;
812 }
3e87ef69 813
6b8f6f6e 814
815 Int_t npars = n_overlaps * NUM_PARS;
3e87ef69 816 if(fDebug)
817 cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
6b8f6f6e 818 Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
3e87ef69 819
820 if(ret<0)
6b8f6f6e 821 {
3e87ef69 822 SetClusterfitFalse(track);
823 fFailed++;
824 fFitError++;
825 delete [] plane;
826 return;
827 //exit(5);
6b8f6f6e 828 }
829
830 chisq_f /= (ndata-fit_pars);
3e87ef69 831 if(fDebug)
832 cout<<"Chisq "<<chisq_f<<endl;
833
6b8f6f6e 834 k=-1;
835 n_overlaps=0;
836 while(k < track->GetNOverlaps(fCurrentPadRow))
837 {
838 AliL3ModelTrack *tr=0;
839 if(k==-1)
840 tr = track;
841 else
3e87ef69 842 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
6b8f6f6e 843 k++;
844 if(!tr) continue;
6b8f6f6e 845 if(!tr->IsPresent(fCurrentPadRow))
846 {
3e87ef69 847 if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
848
6b8f6f6e 849 if(chisq_f < fChiSqMax)//cluster fit is good enough
850 {
851 tot_charge = (Int_t)(a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
3e87ef69 852 Float_t fpad = a[n_overlaps*NUM_PARS+2];
853 Float_t ftime = a[n_overlaps*NUM_PARS+4];
854 if(tot_charge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) ||
855 ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
856 {
857 if(fDebug)
858 cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
859 <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<tot_charge<<" xywidth "<<a[n_overlaps*NUM_PARS+3]
860 <<" zwidth "<<a[n_overlaps*NUM_PARS+5]<<" peakcharge "<<a[n_overlaps*NUM_PARS+1]<<endl;
861 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
862 fFailed++;
863 fResultError++;
864 continue;
865 }
866
867 tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,0,0,pad_num);
868 if(fDebug)
869 cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]<<" charge "<<tot_charge<<endl;
870 /*
871 //Set the digits to used:
872 for(Int_t i=padrange[0]; i<=padrange[1]; i++)
873 for(Int_t j=timerange[0]; j<=timerange[1]; j++)
874 fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
875 */
876 fFitted++;
6b8f6f6e 877 }
878 else //fit was too bad
3e87ef69 879 {
880 if(fDebug)
881 cout<<"Cluster fit was too bad"<<endl;
882 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
883 fBadFitError++;
884 fFailed++;
885 }
6b8f6f6e 886 }
887 n_overlaps++;
888 }
889
890 delete [] plane;
6b8f6f6e 891}
892
3e87ef69 893void AliL3ClusterFitter::SetClusterfitFalse(AliL3ModelTrack *track)
894{
895 //Cluster fit failed, so set the clusters to all the participating
896 //tracks to zero.
897
898 Int_t i=-1;
899 Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
900 while(i < track->GetNOverlaps(fCurrentPadRow))
901 {
902 AliL3ModelTrack *tr=0;
903 if(i==-1)
904 tr = track;
905 else
906 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[i]);
907 i++;
908 if(!tr) continue;
909
910 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
911 }
912}
913
914
915void AliL3ClusterFitter::AddClusters()
916{
917 if(!fClusters)
918 {
919 fClusters = new AliL3SpacePointData[fNMaxClusters];
920 fNClusters=0;
921 }
922
923 if(fDebug)
924 cout<<"Writing cluster in slice "<<fSlice<<" patch "<<fPatch<<endl;
925
926 AliL3TrackArray *tracks=0;
927 if(fSeeding==kTRUE)
928 tracks = fSeeds;
929 else
930 tracks = fTracks;
931
932 if(!tracks)
933 return;
934
935 for(Int_t i=0; i<tracks->GetNTracks(); i++)
936 {
937 AliL3ModelTrack *tr = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
938 if(!tr) continue;
939
940 UInt_t *hitids = tr->GetHitNumbers();
941 Int_t nhits = tr->GetNHits();
942 for(Int_t i=fRowMax; i>=fRowMin; i--)
943 {
944 if(fSeeding)
945 if(tr->GetClusterModel(i)->fSlice != fSlice) continue;
946 if(!tr->IsPresent(i)) continue;
947 fCurrentPadRow = i;
948 Float_t pad,time,xywidth,zwidth;
949 Int_t charge;
950 tr->GetPad(i,pad);
951 tr->GetTime(i,time);
952 tr->GetClusterCharge(i,charge);
953
954 if(pad < -1 || pad >= AliL3Transform::GetNPads(i) ||
955 time < -1 || time >= AliL3Transform::GetNTimeBins())
956 {
957 continue;
958 cout<<"slice "<<fSlice<<" row "<<i<<" pad "<<pad<<" time "<<time<<endl;
959 tr->Print();
960 exit(5);
961 }
962
963 tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
964
965 tr->GetXYWidth(i,xywidth);
966 tr->GetZWidth(i,zwidth);
967 Float_t xyz[3];
968 Int_t sector,row;
969 AliL3Transform::Slice2Sector(fSlice,i,sector,row);
970
971 AliL3Transform::Raw2Global(xyz,sector,row,pad,time);
972
973 if(fNClusters >= fNMaxClusters)
974 {
975 cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
976 exit(5);
977 }
978 fClusters[fNClusters].fX = xyz[0];
979 fClusters[fNClusters].fY = xyz[1];
980 fClusters[fNClusters].fZ = xyz[2];
981 fClusters[fNClusters].fCharge = charge;
982 fClusters[fNClusters].fPadRow = i;
983 Int_t pa = AliL3Transform::GetPatch(i);
984 if(xywidth==0 || zwidth==0)
985 cerr<<"AliL3ClusterFitter::AddClusters : Cluster with zero width"<<endl;
986 if(xywidth>0)
987 fClusters[fNClusters].fSigmaY2 = xywidth*pow(AliL3Transform::GetPadPitchWidth(pa),2);
988 else
989 fClusters[fNClusters].fSigmaY2 = 1;
990 if(zwidth>0)
991 fClusters[fNClusters].fSigmaZ2 = zwidth*pow(AliL3Transform::GetZWidth(),2);
992 else
993 fClusters[fNClusters].fSigmaZ2 = 1;
994 Int_t pat=fPatch;
995 if(fPatch==-1)
996 pat=0;
997 fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
998
999 if(nhits >= 159)
1000 {
1001 cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
1002 exit(5);
1003 }
1004 hitids[nhits++] = fClusters[fNClusters].fID;
1005
1006#ifdef do_mc
1007 Int_t trackID[3];
1008 Int_t fpad = (Int_t)rint(pad);
1009 Int_t ftime = (Int_t)rint(time);
1010 if(fpad < 0)
1011 fpad=0;
1012 if(fpad >= AliL3Transform::GetNPads(i))
1013 fpad = AliL3Transform::GetNPads(i)-1;
1014 if(ftime<0)
1015 ftime=0;
1016 if(ftime >= AliL3Transform::GetNTimeBins())
1017 ftime = AliL3Transform::GetNTimeBins()-1;
1018 GetTrackID(fpad,ftime,trackID);
1019 fClusters[fNClusters].fTrackID[0] = trackID[0];
1020 fClusters[fNClusters].fTrackID[1] = trackID[1];
1021 fClusters[fNClusters].fTrackID[2] = trackID[2];
1022#endif
1023 //cout<<"Setting id "<<trackID[0]<<" on pad "<<pad<<" time "<<time<<" row "<<i<<endl;
1024 fNClusters++;
1025 }
1026
1027 //Copy back the number of assigned clusters
1028 tr->SetNHits(nhits);
1029 }
1030}
1031
1032void AliL3ClusterFitter::WriteTracks()
1033{
1034 if(!fSeeds)
1035 return;
1036
1037 AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
1038 c->WriteFile(fSeeds,"tracks_after.raw");
1039 delete c;
1040
1041 Int_t clustercount=0;
1042 for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
1043 {
1044 AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
1045 if(!tr) continue;
1046 if(tr->GetNHits()==0)
1047 fSeeds->Remove(i);
1048 clustercount += tr->GetNHits();
1049 /*
1050 if(tr->GetPt() > 1 && tr->GetNPresentClusters() < 150)
1051 tr->Print();
1052 */
1053 }
1054 cout<<"Writing "<<clustercount<<" clusters"<<endl;
1055 fSeeds->Compress();
1056 AliL3MemHandler mem;
1057 Char_t filename[1024];
1058 sprintf(filename,"%s/fitter/tracks_0.raw",fPath);
1059 mem.SetBinaryOutput(filename);
1060 mem.TrackArray2Binary(fSeeds);
1061 mem.CloseBinaryOutput();
1062
1063}
1064
1065void AliL3ClusterFitter::WriteClusters()
1066{
1067 AliL3MemHandler mem;
1068 if(fDebug)
1069 cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
1070 Char_t filename[1024];
1071 sprintf(filename,"%s/fitter/points_0_%d_%d.raw",fPath,fSlice,fPatch);
1072 mem.SetBinaryOutput(filename);
1073 mem.Memory2Binary(fNClusters,fClusters);
1074 mem.CloseBinaryOutput();
1075 mem.Free();
1076
1077 delete [] fClusters;
1078 fClusters=0;
1079 fNClusters=0;
1080}