HLT/comp/AliHLTClusterFitter.cxx

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