3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliL3StandardIncludes.h"
8 #include "AliL3ClusterFitter.h"
9 #include "AliL3FitUtilities.h"
10 #include "AliL3DigitData.h"
11 #include "AliL3ModelTrack.h"
12 #include "AliL3TrackArray.h"
13 #include "AliL3MemHandler.h"
14 #include "AliL3HoughTrack.h"
15 #include "AliL3SpacePointData.h"
21 /** /class AliL3ClusterFitter
23 //_____________________________________________________________
30 ClassImp(AliL3ClusterFitter)
32 Int_t AliL3ClusterFitter::fgBadFitError=0;
33 Int_t AliL3ClusterFitter::fgFitError=0;
34 Int_t AliL3ClusterFitter::fgResultError=0;
35 Int_t AliL3ClusterFitter::fgFitRangeError=0;
37 AliL3ClusterFitter::AliL3ClusterFitter()
41 fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
58 AliL3ClusterFitter::AliL3ClusterFitter(Char_t *path)
63 fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
80 AliL3ClusterFitter::~AliL3ClusterFitter()
88 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch,Int_t *rowrange,AliL3TrackArray *tracks)
90 //Assuming tracklets found by the line transform
95 if(rowrange[0] > AliL3Transform::GetLastRow(patch) || rowrange[1] < AliL3Transform::GetFirstRow(patch))
96 cerr<<"AliL3ClusterFitter::Init : Wrong rows "<<rowrange[0]<<" "<<rowrange[1]<<endl;
102 if(fRowMax > AliL3Transform::GetLastRow(fPatch))
103 fRowMax = AliL3Transform::GetLastRow(fPatch);
107 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
108 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
109 Int_t bounds = ntimes*npads;
112 fRow = new Digit[bounds];
116 fTracks = new AliL3TrackArray("AliL3ModelTrack");
118 for(Int_t i=0; i<tracks->GetNTracks(); i++)
120 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
122 AliL3ModelTrack *mtrack = (AliL3ModelTrack*)fTracks->NextTrack();
123 mtrack->Init(slice,patch);
124 mtrack->SetTgl(track->GetTgl());
125 mtrack->SetRowRange(rowrange[0],rowrange[1]);
126 for(Int_t j=fRowMin; j<=fRowMax; j++)
129 track->GetLineCrossingPoint(j,hit);
130 hit[0] += AliL3Transform::Row2X(track->GetFirstRow());
131 Float_t r = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
132 hit[2] = r*track->GetTgl();
134 AliL3Transform::Slice2Sector(slice,j,se,ro);
135 AliL3Transform::Local2Raw(hit,se,ro);
136 if(hit[1]<0 || hit[1]>=AliL3Transform::GetNPads(j) || hit[2]<0 || hit[2]>=AliL3Transform::GetNTimeBins())
138 mtrack->SetPadHit(j,-1);
139 mtrack->SetTimeHit(j,-1);
142 mtrack->SetPadHit(j,hit[1]);
143 mtrack->SetTimeHit(j,hit[2]);
144 mtrack->SetCrossingAngleLUT(j,fabs(track->GetPsiLine() - AliL3Transform::Pi()/2));
145 //if(mtrack->GetCrossingAngleLUT(j) > AliL3Transform::Deg2Rad(20))
146 // cout<<"Angle "<<mtrack->GetCrossingAngleLUT(j)<<" psiline "<<track->GetPsiLine()*180/3.1415<<endl;
147 mtrack->CalculateClusterWidths(j);
150 // cout<<"Copied "<<fTracks->GetNTracks()<<" tracks "<<endl;
153 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch)
158 fRowMin=AliL3Transform::GetFirstRow(patch);
159 fRowMax=AliL3Transform::GetLastRow(patch);
163 Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
164 Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
165 Int_t bounds = ntimes*npads;
168 fRow = new Digit[bounds];
171 fTracks = new AliL3TrackArray("AliL3ModelTrack");
175 void AliL3ClusterFitter::LoadLocalSegments()
177 Char_t filename[1024];
178 sprintf(filename,"%s/hough/tracks_ho_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
180 mem.SetBinaryInput(filename);
181 mem.Binary2TrackArray(fTracks);
182 mem.CloseBinaryInput();
184 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
186 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
189 track->CalculateHelix();
191 track->Init(fSlice,fPatch);
193 for(Int_t j=fRowMin; j<=fRowMax; j++)
195 //Calculate the crossing point between track and padrow
198 if(!track->GetCrossingPoint(j,xyzCross))
202 AliL3Transform::Slice2Sector(fSlice,j,sector,row);
203 AliL3Transform::Local2Raw(xyzCross,sector,row);
205 if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j) ||
206 xyzCross[2] < 0 || xyzCross[2] >= AliL3Transform::GetNTimeBins()) //track goes out of range
209 track->SetPadHit(j,xyzCross[1]);
210 track->SetTimeHit(j,xyzCross[2]);
212 Float_t crossingangle = track->GetCrossingAngle(j);
213 track->SetCrossingAngleLUT(j,crossingangle);
214 track->CalculateClusterWidths(j);
215 track->GetClusterModel(j)->fSlice = fSlice;
221 void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline,Int_t eventnr,Float_t zvertex)
223 //Function assumes _global_ tracks written to a single file.
224 cout<<"Loading the seeds"<<endl;
229 sprintf(fname,"%s/offline/tracks_%d.raw",fPath,fEvent);
231 sprintf(fname,"%s/hough/tracks_%d.raw",fPath,fEvent);
233 cout<<"AliL3ClusterFitter::LoadSeeds : Loading input tracks from "<<fname<<endl;
235 AliL3MemHandler tfile;
236 tfile.SetBinaryInput(fname);
240 fSeeds = new AliL3TrackArray("AliL3ModelTrack");
242 tfile.Binary2TrackArray(fSeeds);
243 tfile.CloseBinaryInput();
248 Int_t clustercount=0;
249 for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
251 AliL3ModelTrack *track = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
256 if(i==0) cerr<<"AliL3ClusterFitter::LoadSeeds : Cutting on pt of 4 GeV!!"<<endl;
257 if(track->GetPt() > 4.)
263 clustercount += track->GetNHits();
264 track->CalculateHelix();
266 Int_t nhits = track->GetNHits();
267 UInt_t *hitids = track->GetHitNumbers();
269 Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
271 if(i==0) cerr<<"Cluster fitter only in HALF TPC!!!"<<endl;
278 track->Init(origslice,-1);
279 Int_t slice = origslice;
281 //for(Int_t j=rowrange[1]; j>=rowrange[0]; j--)
282 for(Int_t j=rowrange[0]; j<=rowrange[1]; j++)
285 //Calculate the crossing point between track and padrow
286 Float_t angle = 0; //Perpendicular to padrow in local coordinates
287 AliL3Transform::Local2GlobalAngle(&angle,slice);
288 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
290 //cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
295 Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
296 xyzCross[2] += zvertex;
299 AliL3Transform::Slice2Sector(slice,j,sector,row);
300 AliL3Transform::Global2Raw(xyzCross,sector,row);
301 //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyzCross[1]<<" time "<<xyzCross[2]<<endl;
302 if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
307 Float_t lastcross=xyzCross[1];
326 if(slice < 0 || slice>35)
328 cerr<<"Wrong slice "<<slice<<" on row "<<j<<endl;
331 //cout<<"Track leaving, trying slice "<<slice<<endl;
333 AliL3Transform::Local2GlobalAngle(&angle,slice);
334 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
336 cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
341 xyzCross[0] = track->GetPointX();
342 xyzCross[1] = track->GetPointY();
343 xyzCross[2] = track->GetPointZ();
344 xyzCross[2] += zvertex;
346 AliL3Transform::Slice2Sector(slice,j,sector,row);
347 AliL3Transform::Global2Raw(xyzCross,sector,row);
348 if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
350 if(xyzCross[1] > 0 && lastcross > 0 || xyzCross[1] < 0 && lastcross < 0)
354 slice = tslice;//Track is on the border of two slices
360 if(xyzCross[2] < 0 || xyzCross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
363 if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j))
365 cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyzCross[1]<<" time "<<xyzCross[2]<<endl;
370 track->SetPadHit(j,xyzCross[1]);
371 track->SetTimeHit(j,xyzCross[2]);
373 AliL3Transform::Local2GlobalAngle(&angle,slice);
374 Float_t crossingangle = track->GetCrossingAngle(j,slice);
375 track->SetCrossingAngleLUT(j,crossingangle);
377 track->CalculateClusterWidths(j);
379 track->GetClusterModel(j)->fSlice = slice;
382 memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));//Reset the hitnumbers
387 cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
390 void AliL3ClusterFitter::FindClusters()
394 cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
399 cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
403 AliL3DigitRowData *rowPt = fRowData;
404 AliL3DigitData *digPt=0;
411 fRowMin = AliL3Transform::GetFirstRow(fPatch);
412 fRowMax = AliL3Transform::GetLastRow(fPatch);
414 for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
416 if((Int_t)rowPt->fRow < fRowMin)
418 AliL3MemHandler::UpdateRowPointer(rowPt);
421 else if((Int_t)rowPt->fRow > fRowMax)
423 else if((Int_t)rowPt->fRow != i)
425 cerr<<"AliL3ClusterFitter::FindClusters : Mismatching row numbering "<<i<<" "<<rowPt->fRow<<endl;
429 memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
430 digPt = (AliL3DigitData*)rowPt->fDigitData;
432 for(UInt_t j=0; j<rowPt->fNDigit; j++)
435 time = digPt[j].fTime;
436 charge = digPt[j].fCharge;
437 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
438 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
439 //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
442 for(Int_t it=0; it<2; it++)
446 fProcessTracks = fSeeds;
451 fProcessTracks = fTracks;
457 for(Int_t k=0; k<fProcessTracks->GetNTracks(); k++)
459 AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(k);
463 if(track->GetClusterModel(i)->fSlice != fSlice) continue;
465 if(track->GetPadHit(i) < 0 || track->GetPadHit(i) > AliL3Transform::GetNPads(i)-1 ||
466 track->GetTimeHit(i) < 0 || track->GetTimeHit(i) > AliL3Transform::GetNTimeBins()-1)
468 track->SetCluster(i,0,0,0,0,0,0);
472 if(CheckCluster(k) == kFALSE)
476 AliL3MemHandler::UpdateRowPointer(rowPt);
484 cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
485 cout<<"Distribution:"<<endl;
486 cout<<"Bad fit "<<fgBadFitError<<endl;
487 cout<<"Fit error "<<fgFitError<<endl;
488 cout<<"Result error "<<fgResultError<<endl;
489 cout<<"Fit range error "<<fgFitRangeError<<endl;
493 Bool_t AliL3ClusterFitter::CheckCluster(Int_t trackindex)
495 //Check if this is a single or overlapping cluster
497 AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(trackindex);
499 Int_t row = fCurrentPadRow;
501 if(track->IsSet(row)) //A fit has already be performed on this one
504 //Define the cluster region of this hit:
505 Int_t padr[2]={999,-1};
506 Int_t timer[2]={999,-1};
508 if(!SetFitRange(track,padr,timer))
510 track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
513 cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
514 <<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
515 <<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
520 //Check if any other track contributes to this cluster:
522 for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
524 AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
527 if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
529 if(tr->GetPadHit(row) > padr[0] && tr->GetPadHit(row) < padr[1] &&
530 tr->GetTimeHit(row) > timer[0] && tr->GetTimeHit(row) < timer[1])
532 if(SetFitRange(tr,padr,timer))
533 track->SetOverlap(row,t);
536 Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))*GetYWidthFactor());
537 Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))*GetZWidthFactor());
539 (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
540 tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
542 (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
543 tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
545 (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
546 tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) ||
548 (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
549 tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1])
553 if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
555 track->SetOverlap(row,t); //Set overlap
562 cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
564 FitClusters(track,padr,timer);
565 //CalculateWeightedMean(track,padr,timer);
569 Bool_t AliL3ClusterFitter::SetFitRange(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
571 Int_t row = fCurrentPadRow;
572 Int_t nt = AliL3Transform::GetNTimeBins()+1;
575 if(row < AliL3Transform::GetNRowLow())
577 else if(row < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
583 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
584 ,0,1,2,3,3,3,3,3,3,3
586 ,-3,-3,-3,-3,-3,-3,-1,-1};
587 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
588 ,-3,-3,-3,-3,-2,-1,0,1,2,3
589 ,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
592 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
593 ,-3,3,-3,3,-3,3,-3,3,-3,3
594 ,0,0,-1,-1,1,1,-2,-2,2,2,-3,-3,3,3};
595 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
596 ,0,0,-1,-1,1,1,-2,-2,2,2
597 ,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3};
599 Int_t padhit = (Int_t)rint(track->GetPadHit(row));
600 Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
604 for(Int_t index=0; index<nsearchbins; index++)
606 if(IsMaximum(padhit + padloop[index],timehit + timeloop[index]))
608 padmax = padhit + padloop[index];
609 timemax = timehit + timeloop[index];
615 //Define the cluster region of this hit:
616 //The region we look for, is centered at the local maxima
617 //and expanded around using the parametrized cluster width
618 //according to track parameters.
623 if(padmax>=0 && timemax>=0)
625 //Set the hit to the local maxima of the cluster.
626 //Store the maxima in the cluster model structure,
627 //-only temporary, it will be overwritten when calling SetCluster.
629 track->GetClusterModel(row)->fDPad = padmax;
630 track->GetClusterModel(row)->fDTime = timemax;
632 Int_t i=padmax,j=timemax;
633 for(Int_t pdir=-1; pdir<=1; pdir+=2)
636 while(abs(padmax-i) < xyw)
638 Bool_t chargeonpad=kFALSE;
639 for(Int_t tdir=-1; tdir<=1; tdir+=2)
642 while(abs(timemax-j) < zw)
644 if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) break;
645 if(fRow[nt*i+j].fCharge)
647 if(i < padrange[0]) padrange[0]=i;
648 if(i > padrange[1]) padrange[1]=i;
649 if(j < timerange[0]) timerange[0]=j;
650 if(j > timerange[1]) timerange[1]=j;
664 for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
666 for(Int_t j=timemax-zw; j<=timemax+zw; j++)
668 if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
669 if(fRow[nt*i+j].fCharge)
671 if(i < padrange[0]) padrange[0]=i;
672 if(i > padrange[1]) padrange[1]=i;
673 if(j < timerange[0]) timerange[0]=j;
674 if(j > timerange[1]) timerange[1]=j;
681 cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
687 Bool_t AliL3ClusterFitter::IsMaximum(Int_t pad,Int_t time)
689 if(pad<0 || pad >= AliL3Transform::GetNPads(fCurrentPadRow) ||
690 time<0 || time >= AliL3Transform::GetNTimeBins())
692 Int_t nt = AliL3Transform::GetNTimeBins()+1;
693 if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
694 Int_t charge = fRow[nt*pad+time].fCharge;
695 if(charge == 1024 || charge==0) return kFALSE;
696 //if(charge == 0) return kFALSE;
698 //fRow[nt*pad+time].fUsed = kTRUE;
701 if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
702 if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
703 if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
704 if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
705 if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
706 if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
707 if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
708 if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
709 fRow[nt*pad+time].fUsed = kTRUE;
713 void AliL3ClusterFitter::CalculateWeightedMean(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
717 Float_t pad=0,time=0;
718 Int_t nt = AliL3Transform::GetNTimeBins()+1;
719 for(Int_t i=padrange[0]; i<=padrange[1]; i++)
722 for(Int_t j=timerange[0]; j<=timerange[1]; j++)
724 lsum += fRow[nt*(i-1)+(j-1)].fCharge;
725 time += j*fRow[nt*(i-1)+(j-1)].fCharge;
735 track->SetCluster(fCurrentPadRow,pad,time,sum,0,0,npads);
738 track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
741 void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
743 //Handle single and overlapping clusters
746 if( (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) <= 1 ||
747 (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) >= AliL3Transform::GetNPads(fCurrentPadRow)-2 ||
748 (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) <= 1 ||
749 (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) >= AliL3Transform::GetNTimeBins()-2)
751 CalculateWeightedMean(track,padrange,timerange);
755 Int_t size = FIT_PTS;
756 Int_t maxTracks = FIT_MAXPAR/NUM_PARS;
757 if(track->GetNOverlaps(fCurrentPadRow) > maxTracks)
759 cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
762 Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
764 //Check if at least one cluster is not already fitted
765 Bool_t allFitted=kTRUE;
768 while(k < track->GetNOverlaps(fCurrentPadRow))
770 AliL3ModelTrack *tr=0;
774 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
777 if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
786 cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
790 //Allocate fit parameters array; this is interface to the C code
791 plane = new DPOINT[FIT_PTS];
792 memset(plane,0,FIT_PTS*sizeof(DPOINT));
794 Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
796 //Fill the fit parameters:
797 Double_t a[FIT_MAXPAR];
798 Int_t lista[FIT_MAXPAR];
799 Double_t dev[FIT_MAXPAR],chisqF;
806 //Fill the overlapping tracks:
807 while(k < track->GetNOverlaps(fCurrentPadRow))
809 AliL3ModelTrack *tr=0;
813 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
817 if(tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow)) continue;//Cluster fit failed before
819 //Use the local maxima as the input to the fitting routine.
820 //The local maxima is temporary stored in the cluster model:
821 Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
822 Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
823 Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
826 cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
827 <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" at local maxima in pad "<<hitpad
828 <<" time "<<hittime<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
829 <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
833 cerr<<"Charge still zero!"<<endl;
837 a[nOverlaps*NUM_PARS+2] = hitpad;
838 a[nOverlaps*NUM_PARS+4] = hittime;
840 if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
842 a[nOverlaps*NUM_PARS+1] = charge;
843 a[nOverlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
844 a[nOverlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
845 //a[nOverlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
846 lista[nOverlaps*NUM_PARS + 1] = 1;
847 lista[nOverlaps*NUM_PARS + 2] = 1;
848 lista[nOverlaps*NUM_PARS + 3] = 0;
849 lista[nOverlaps*NUM_PARS + 4] = 1;
850 lista[nOverlaps*NUM_PARS + 5] = 0;
851 //lista[nOverlaps*NUM_PARS + 6] = 0;
852 fitPars += 3; //<-------------------
854 else //Cluster was fitted before
856 if(!tr->IsPresent(fCurrentPadRow))
858 cerr<<"AliL3ClusterFitter::FindClusters : Cluster not present; there is a bug here"<<endl;
862 Float_t xywidth,zwidth,pad,time;
863 tr->GetPad(fCurrentPadRow,pad);
864 tr->GetTime(fCurrentPadRow,time);
865 tr->GetClusterCharge(fCurrentPadRow,charge);
866 xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
867 zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
869 cout<<endl<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
871 a[nOverlaps*NUM_PARS+2] = pad;
872 a[nOverlaps*NUM_PARS+4] = time;
873 a[nOverlaps*NUM_PARS+1] = charge;
874 a[nOverlaps*NUM_PARS+3] = xywidth * GetYWidthFactor();
875 a[nOverlaps*NUM_PARS+5] = zwidth * GetZWidthFactor();
876 //a[nOverlaps*NUM_PARS+6] = zwidth * GetZWidthFactor();
878 lista[nOverlaps*NUM_PARS + 1] = 1;
879 lista[nOverlaps*NUM_PARS + 2] = 0;
880 lista[nOverlaps*NUM_PARS + 3] = 0;
881 lista[nOverlaps*NUM_PARS + 4] = 0;
882 lista[nOverlaps*NUM_PARS + 5] = 0;
883 //lista[nOverlaps*NUM_PARS + 6] = 0;
889 if(nOverlaps==0) //No clusters here
900 cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
901 for(Int_t i=padrange[0]; i<=padrange[1]; i++)
905 for(Int_t j=timerange[0]; j<=timerange[1]; j++)
907 Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
909 if(charge <= 0) continue;
912 if(charge > maxCharge)
918 cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
923 cerr<<"Too many points; row "<<fCurrentPadRow<<" padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "
924 <<timerange[0]<<" "<<timerange[1]<<endl;
928 plane[ndata].u = (Double_t)i;
929 plane[ndata].v = (Double_t)j;
932 s[ndata]= 1 + sqrt((Double_t)charge);
934 if(maxCharge) //there was charge on this pad
936 if(timeNumMax < timeNum)
937 timeNumMax = timeNum;
940 if(padNum <= 1 || timeNumMax <=1 || nOverlaps > fNmaxOverlaps || ndata <= fitPars) //too few to do fit
942 SetClusterfitFalse(track);
944 cout<<"Too few digits or too many overlaps: "<<padNum<<" "<<timeNumMax<<" "<<nOverlaps<<" ndata "<<ndata<<" fitPars "<<fitPars<<endl;
950 Int_t npars = nOverlaps * NUM_PARS;
952 cout<<"Number of overlapping clusters "<<nOverlaps<<endl;
953 Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisqF, f2gauss5 );
957 SetClusterfitFalse(track);
965 chisqF /= (ndata-fitPars);
967 cout<<"Chisq "<<chisqF<<endl;
969 Bool_t overlapping=kFALSE;
970 if(track->GetNOverlaps(fCurrentPadRow) > 0)//There is a overlap
975 while(k < track->GetNOverlaps(fCurrentPadRow))
977 AliL3ModelTrack *tr=0;
981 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
984 if(!tr->IsPresent(fCurrentPadRow))
986 if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
989 if(fCurrentPadRow < AliL3Transform::GetNRowLow())
991 else if(fCurrentPadRow < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
996 //if(chisqF < fChiSqMax[(Int_t)overlapping])//cluster fit is good enough
997 if(chisqF < fChiSqMax[lpatch])//cluster fit is good enough
999 totCharge = (Int_t)(2*AliL3Transform::Pi() * a[nOverlaps*NUM_PARS+1] * a[nOverlaps*NUM_PARS+3] * a[nOverlaps*NUM_PARS+5]);
1000 Float_t fpad = a[nOverlaps*NUM_PARS+2];
1001 Float_t ftime = a[nOverlaps*NUM_PARS+4];
1002 if(totCharge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) ||
1003 ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
1006 cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
1007 <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<totCharge<<" xywidth "<<a[nOverlaps*NUM_PARS+3]
1008 <<" zwidth "<<a[nOverlaps*NUM_PARS+5]<<" peakcharge "<<a[nOverlaps*NUM_PARS+1]<<endl;
1009 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1015 tr->SetCluster(fCurrentPadRow,fpad,ftime,totCharge,0,0,padNum);
1017 tr->SetCluster(fCurrentPadRow,fpad,ftime,totCharge,
1018 pow(a[nOverlaps*NUM_PARS+3],2),
1019 pow(a[nOverlaps*NUM_PARS+5],2),padNum);
1023 cout<<"Setting cluster in slice "<<fSlice<<" row "<<fCurrentPadRow<<" pad "<<a[nOverlaps*NUM_PARS+2]<<" time "<<a[nOverlaps*NUM_PARS+4]
1024 <<" padwidth "<<a[nOverlaps*NUM_PARS+3]<<" timewidth "<<a[nOverlaps*NUM_PARS+5]
1025 <<" charge "<<totCharge<<endl;
1028 //Set the digits to used:
1029 for(Int_t i=padrange[0]; i<=padrange[1]; i++)
1030 for(Int_t j=timerange[0]; j<=timerange[1]; j++)
1031 fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
1035 else //fit was too bad
1038 cout<<"Cluster fit was too bad"<<endl;
1039 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1050 void AliL3ClusterFitter::SetClusterfitFalse(AliL3ModelTrack *track)
1052 //Cluster fit failed, so set the clusters to all the participating
1056 Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
1057 while(i < track->GetNOverlaps(fCurrentPadRow))
1059 AliL3ModelTrack *tr=0;
1063 tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[i]);
1067 //Set the digit data to unused, so it can be fitted to another bastard:
1068 Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
1069 Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
1070 fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fUsed = kFALSE;
1072 tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1078 void AliL3ClusterFitter::AddClusters()
1082 fClusters = new AliL3SpacePointData[fNMaxClusters];
1087 cout<<"Writing cluster in slice "<<fSlice<<" patch "<<fPatch<<endl;
1089 AliL3TrackArray *tracks=0;
1098 for(Int_t i=0; i<tracks->GetNTracks(); i++)
1100 AliL3ModelTrack *tr = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
1103 UInt_t *hitids = tr->GetHitNumbers();
1104 Int_t nhits = tr->GetNHits();
1105 for(Int_t i=fRowMax; i>=fRowMin; i--)
1108 if(tr->GetClusterModel(i)->fSlice != fSlice) continue;
1109 if(!tr->IsPresent(i)) continue;
1111 Float_t pad,time,xywidth,zwidth;
1114 tr->GetTime(i,time);
1115 tr->GetClusterCharge(i,charge);
1117 if(pad < -1 || pad >= AliL3Transform::GetNPads(i) ||
1118 time < -1 || time >= AliL3Transform::GetNTimeBins())
1121 cout<<"slice "<<fSlice<<" row "<<i<<" pad "<<pad<<" time "<<time<<endl;
1126 tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
1128 tr->GetSigmaY2(i,xywidth);
1129 tr->GetSigmaZ2(i,zwidth);
1132 AliL3Transform::Slice2Sector(fSlice,i,sector,row);
1134 AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
1136 if(fNClusters >= fNMaxClusters)
1138 cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
1142 fClusters[fNClusters].fX = xyz[0];
1143 fClusters[fNClusters].fY = xyz[1];
1144 fClusters[fNClusters].fZ = xyz[2];
1145 fClusters[fNClusters].fCharge = charge;
1146 fClusters[fNClusters].fPadRow = i;
1147 Int_t pa = AliL3Transform::GetPatch(i);
1148 if(xywidth==0 || zwidth==0)
1149 cerr<<"AliL3ClusterFitter::AddClusters : Cluster with zero width"<<endl;
1151 fClusters[fNClusters].fSigmaY2 = xywidth*pow(AliL3Transform::GetPadPitchWidth(pa),2);
1153 fClusters[fNClusters].fSigmaY2 = 1;
1155 fClusters[fNClusters].fSigmaZ2 = zwidth*pow(AliL3Transform::GetZWidth(),2);
1157 fClusters[fNClusters].fSigmaZ2 = 1;
1161 fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
1163 if(nhits >= AliL3Transform::GetNRows())
1165 cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
1168 hitids[nhits++] = fClusters[fNClusters].fID;
1172 Int_t fpad = (Int_t)rint(pad);
1173 Int_t ftime = (Int_t)rint(time);
1176 if(fpad >= AliL3Transform::GetNPads(i))
1177 fpad = AliL3Transform::GetNPads(i)-1;
1180 if(ftime >= AliL3Transform::GetNTimeBins())
1181 ftime = AliL3Transform::GetNTimeBins()-1;
1182 GetTrackID(fpad,ftime,trackID);
1183 fClusters[fNClusters].fTrackID[0] = trackID[0];
1184 fClusters[fNClusters].fTrackID[1] = trackID[1];
1185 fClusters[fNClusters].fTrackID[2] = trackID[2];
1187 //cout<<"Setting id "<<trackID[0]<<" on pad "<<pad<<" time "<<time<<" row "<<i<<endl;
1191 //Copy back the number of assigned clusters
1192 tr->SetNHits(nhits);
1197 void AliL3ClusterFitter::WriteTracks(Int_t minHits)
1202 AliL3TrackArray *fakes = new AliL3TrackArray();
1204 Int_t clustercount=0;
1205 for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
1207 AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
1209 if(tr->GetNHits() < minHits)
1214 clustercount += tr->GetNHits();
1217 cout<<"Writing "<<clustercount<<" clusters"<<endl;
1219 AliL3MemHandler mem;
1220 Char_t filename[1024];
1221 sprintf(filename,"%s/fitter/tracks_%d.raw",fPath,fEvent);
1222 mem.SetBinaryOutput(filename);
1223 mem.TrackArray2Binary(fSeeds);
1224 mem.CloseBinaryOutput();
1226 //Write the fake tracks to its own file
1228 sprintf(filename,"%s/fitter/tracks_fakes_%d.raw",fPath,fEvent);
1229 mem.SetBinaryOutput(filename);
1230 mem.TrackArray2Binary(fakes);
1231 mem.CloseBinaryOutput();
1235 void AliL3ClusterFitter::WriteClusters(Bool_t global)
1237 AliL3MemHandler mem;
1239 cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
1240 Char_t filename[1024];
1241 sprintf(filename,"%s/fitter/points_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
1242 mem.SetBinaryOutput(filename);
1244 mem.Transform(fNClusters,fClusters,fSlice);
1245 mem.Memory2Binary(fNClusters,fClusters);
1246 mem.CloseBinaryOutput();
1249 delete [] fClusters;