HLT/comp/AliL3Modeller.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <stream.h>
   7 #include <iostream.h>
   8 #include <math.h>
   9
  10 #include "AliL3Modeller.h"
  11 #include "AliL3MemHandler.h"
  12 #include "AliL3TrackArray.h"
  13 #include "AliL3ModelTrack.h"
  14 #include "AliL3DigitData.h"
  15 #include "AliL3Transform.h"
  16
  17 #include "AliL3Defs.h"
  18
  19 ClassImp(AliL3Modeller)
  20
  21 AliL3Modeller::AliL3Modeller()
  22 {
  23   fMemHandler=0;
  24   fTracks=0;
  25   fTransform=0;
  26 }
  27
  28
  29 AliL3Modeller::~AliL3Modeller()
  30 {
  31   if(fMemHandler)
  32     delete fMemHandler;
  33   if(fTracks)
  34     delete fTracks;
  35   if(fTransform)
  36     delete fTransform;
  37
  38 }
  39
  40 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path)
  41 {
  42   fSlice = slice;
  43   fPatch = patch;
  44   fPadOverlap=6;
  45   fTimeOverlap=8;
  46   fTrackThreshold=10;
  47   sprintf(fPath,"%s",path);
  48
  49   fTransform = new AliL3Transform();
  50   fTracks = new AliL3TrackArray("AliL3ModelTrack");
  51
  52   Char_t fname[100];
  53   AliL3MemHandler *file = new AliL3MemHandler();
  54   sprintf(fname,"%s/tracks.raw",trackdata);
  55   if(!file->SetBinaryInput(fname))
  56     {
  57       cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
  58       return;
  59     }
  60   file->Binary2TrackArray(fTracks);
  61   file->CloseBinaryInput();
  62   delete file;
  63
  64   fTracks->QSort();
  65   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
  66     {
  67       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
  68       if(!track) continue;
  69       track->Init(fSlice,fPatch);
  70       track->Rotate(fSlice,kTRUE);
  71       track->CalculateHelix();
  72     }
  73
  74   CalculateCrossingPoints();
  75   CheckForOverlaps();
  76
  77   fMemHandler = new AliL3MemHandler();
  78   sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
  79   if(!fMemHandler->SetBinaryInput(fname))
  80     {
  81       cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
  82       return;
  83     }
  84   UInt_t ndigits;
  85   AliL3DigitRowData *digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
  86
  87   SetInputData(digits);
  88 }
  89
  90 void AliL3Modeller::FindClusters()
  91 {
  92
  93   if(!fTracks)
  94     {
  95       cerr<<"AliL3Modeller::Process : No tracks"<<endl;
  96       return;
  97     }
  98   if(!fRowData)
  99     {
 100       cerr<<"AliL3Modeller::Process : No data "<<endl;
 101       return;
 102     }
 103
 104   AliL3DigitRowData *rowPt = fRowData;
 105   AliL3DigitData *digPt=0;
 106
 107   Int_t ntimes = fTransform->GetNTimeBins()+1;
 108   Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
 109   Int_t bounds = ntimes*npads;
 110   Digit *row = new Digit[bounds];
 111
 112   Int_t seq_charge;
 113   Int_t pad,time,index;
 114   Short_t charge;
 115   Cluster cluster;
 116
 117   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
 118     {
 119       memset((void*)row,0,ntimes*npads*sizeof(Digit));
 120       digPt = (AliL3DigitData*)rowPt->fDigitData;
 121       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 122         {
 123           pad = digPt[j].fPad;
 124           time = digPt[j].fTime;
 125           charge = digPt[j].fCharge;
 126           row[ntimes*pad+time].fCharge = charge;
 127           row[ntimes*pad+time].fUsed = kFALSE;
 128           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 129         }
 130
 131       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 132         {
 133           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
 134           if(!track) continue;
 135           //if(track->GetOverlap(i)>=0) continue;//Track is overlapping
 136
 137           if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
 138             {
 139               track->SetCluster(i,0,0,0,0,0); //The track has left the patch.
 140               continue;
 141             }
 142
 143           Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
 144           Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
 145           //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
 146           pad = hitpad;
 147           time = hittime;
 148           Int_t padsign=-1;
 149           Int_t timesign=-1;
 150
 151           memset(&cluster,0,sizeof(Cluster));
 152
 153           while(1)//Process this padrow
 154             {
 155               if(pad < 0 || pad >= fTransform->GetNPads(i))
 156                 {
 157                   //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
 158                   FillCluster(track,&cluster,i);
 159                   break;
 160                 }
 161               seq_charge=0;
 162               timesign=-1;
 163               time = hittime;
 164               while(1) //Process sequence on this pad:
 165                 {
 166                   if(time < 0) break;
 167                   index = ntimes*pad + time;
 168                   if(index < 0 || index >= bounds)
 169                     {
 170                       cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
 171                         <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
 172                       break;
 173                     }
 174                   charge = row[index].fCharge;
 175                   if(charge==0 && timesign==-1)
 176                     {
 177                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
 178                         {
 179                           time--;
 180                           continue;
 181                         }
 182                       else
 183                         {
 184                           time = hittime+1;
 185                           timesign=1;
 186                           continue;
 187                         }
 188                     }
 189                   else if(charge==0 && timesign==1)
 190                     {
 191                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
 192                         {
 193                           time++;
 194                           continue;
 195                         }
 196                       else
 197                         {
 198                           //cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
 199                           break;
 200                         }
 201                     }
 202
 203                   //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 204
 205                   seq_charge += charge;
 206
 207                   cluster.fTime += time*charge;
 208                   cluster.fPad += pad*charge;
 209                   cluster.fCharge += charge;
 210                   cluster.fSigmaY2 += pad*pad*charge;
 211                   cluster.fSigmaZ2 += time*time*charge;
 212
 213                   row[ntimes*pad+time].fUsed = kTRUE;
 214                   time += timesign;
 215                 }
 216               //cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
 217               if(seq_charge)
 218                 pad += padsign;
 219               else //Nothing more on this pad, goto next pad
 220                 {
 221                   if(padsign==-1)
 222                     {
 223                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad > 0)
 224                         {
 225                           pad--; //In this case, we haven't found anything yet,
 226                         }        //so we will try to expand our search within the natural boundaries.
 227                       else
 228                         {
 229                           pad=hitpad+1;
 230                           padsign=1;
 231                         }
 232                       continue;
 233                     }
 234
 235                   else if(padsign==1)
 236                     {
 237                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad < fTransform->GetNPads(i)-2)
 238                         {
 239                           pad++;     //In this case, we haven't found anything yet,
 240                           continue;  //so we will try to expand our search within the natural boundaries.
 241                         }
 242                       else //We are out of range, or cluster if finished.
 243                         {
 244                           //cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
 245                           FillCluster(track,&cluster,i);
 246                           break;
 247                         }
 248                     }
 249                   else //Nothing more in this cluster
 250                     {
 251                       //cout<<"Filling final cluster"<<endl;
 252                       FillCluster(track,&cluster,i);
 253                       break;
 254                     }
 255                 }
 256             }
 257           //cout<<"done"<<endl;
 258         }
 259
 260       FillZeros(rowPt,row);
 261       fMemHandler->UpdateRowPointer(rowPt);
 262     }
 263   delete [] row;
 264   cout<<"done processing"<<endl;
 265
 266
 267   //Debug:
 268   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 269     {
 270       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 271       if(!track) continue;
 272       if(track->GetNClusters() != NumRows[fPatch])
 273         cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<NumRows[fPatch]<<endl<<endl;
 274     }
 275
 276 }
 277
 278 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row)
 279 {
 280   if(cluster->fCharge==0)
 281     {
 282       track->SetCluster(row,0,0,0,0,0);
 283       return;
 284     }
 285   Float_t fcharge = (Float_t)cluster->fCharge;
 286   Float_t fpad = ((Float_t)cluster->fPad/fcharge);
 287   Float_t ftime = ((Float_t)cluster->fTime/fcharge);
 288   Float_t sigmaY2,sigmaZ2;
 289   CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
 290   track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2);
 291
 292 }
 293
 294 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
 295 {
 296   //Fill zero where data has been used.
 297
 298   Int_t ntimes = fTransform->GetNTimeBins()+1;
 299   AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 300   for(UInt_t j=0; j<rowPt->fNDigit; j++)
 301     {
 302       Int_t pad = digPt[j].fPad;
 303       Int_t time = digPt[j].fTime;
 304       if(row[ntimes*pad+time].fUsed==kTRUE)
 305         digPt[j].fCharge = 0;
 306     }
 307 }
 308
 309 void AliL3Modeller::WriteRemaining()
 310 {
 311   //Write remaining (nonzero) digits to file.
 312
 313   AliL3DigitRowData *rowPt;
 314   rowPt = (AliL3DigitRowData*)fRowData;
 315   Int_t digitcount=0;
 316   Int_t ndigits[(NumRows[fPatch])];
 317   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
 318     {
 319       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 320       ndigits[(i-NRows[fPatch][0])]=0;
 321       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 322         {
 323           if(digPt[j].fCharge==0) continue;
 324
 325           digitcount++;
 326           ndigits[(i-NRows[fPatch][0])]++;
 327         }
 328       //cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
 329       fMemHandler->UpdateRowPointer(rowPt);
 330     }
 331
 332   Int_t size = digitcount*sizeof(AliL3DigitData) + NumRows[fPatch]*sizeof(AliL3DigitRowData);
 333   Byte_t *data = new Byte_t[size];
 334   memset(data,0,size);
 335   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
 336   rowPt = (AliL3DigitRowData*)fRowData;
 337
 338   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
 339     {
 340       Int_t localcount=0;
 341       tempPt->fRow = i;
 342       tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
 343       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 344       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 345         {
 346           if(digPt[j].fCharge==0) continue;
 347           if(localcount >= ndigits[(i-NRows[fPatch][0])])
 348             {
 349               cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
 350               return;
 351             }
 352           tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
 353           tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
 354           tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
 355           localcount++;
 356         }
 357       if(ndigits[(i-NRows[fPatch][0])] != localcount)
 358         {
 359           cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
 360           return;
 361         }
 362       fMemHandler->UpdateRowPointer(rowPt);
 363       Byte_t *tmp = (Byte_t*)tempPt;
 364       Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-NRows[fPatch][0])]*sizeof(AliL3DigitData);
 365       tmp += size;
 366       tempPt = (AliL3DigitRowData*)tmp;
 367     }
 368
 369   Char_t fname[100];
 370   AliL3MemHandler *mem = new AliL3MemHandler();
 371   sprintf(fname,"%s/remains_%d_%d.raw",fPath,fSlice,fPatch);
 372   mem->SetBinaryOutput(fname);
 373   mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
 374   mem->CloseBinaryOutput();
 375   delete mem;
 376   delete [] data;
 377 }
 378
 379
 380 void AliL3Modeller::CalculateCrossingPoints()
 381 {
 382   cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
 383   if(!fTracks)
 384     {
 385       cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
 386       return;
 387     }
 388   Float_t hit[3];
 389
 390
 391   //Remove tracks which are no good:
 392   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 393     {
 394       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 395       if(!track) continue;
 396       if(track->GetFirstPointX() > fTransform->Row2X(NRows[fPatch][1]) || track->GetPt()<0.1)
 397         fTracks->Remove(i);
 398     }
 399
 400   Int_t sector,row;
 401   for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
 402     {
 403       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
 404         {
 405           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 406           if(!track) continue;
 407
 408           if(track->GetNHits() < fTrackThreshold)
 409             {
 410               fTracks->Remove(j);
 411               continue;
 412             }
 413
 414           if(!track->GetCrossingPoint(i,hit))
 415             {
 416               cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
 417                 "First point "<<track->GetFirstPointX()<<
 418                 " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
 419                 //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
 420                 //"--------"<<endl;
 421               fTracks->Remove(j);
 422               continue;
 423             }
 424           //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
 425
 426           fTransform->Slice2Sector(fSlice,i,sector,row);
 427           fTransform->Local2Raw(hit,sector,row);
 428           if(hit[1]<0 || hit[1]>fTransform->GetNPads(i) ||
 429              hit[2]<0 || hit[2]>fTransform->GetNTimeBins())
 430             {//Track is leaving the patch, so flag the track hits (<0)
 431               track->SetPadHit(i,-1);
 432               track->SetTimeHit(i,-1);
 433               continue;
 434             }
 435
 436           track->SetPadHit(i,hit[1]);
 437           track->SetTimeHit(i,hit[2]);
 438
 439           //if(hit[1]<0 || hit[2]>445)
 440           //if(hit[2]<0 || hit[2]>445)
 441           //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
 442           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
 443         }
 444     }
 445   fTracks->Compress();
 446   cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 447 }
 448
 449 void AliL3Modeller::CheckForOverlaps()
 450 {
 451   //Flag the tracks that overlap
 452
 453   cout<<"Checking for overlaps...";
 454
 455   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 456     {
 457       AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 458       if(!track1) continue;
 459       for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
 460         {
 461           AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 462           if(!track2) continue;
 463           for(Int_t k=NRows[fPatch][0]; k<=NRows[fPatch][1]; k++)
 464             {
 465               if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
 466                  track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
 467                 continue;
 468               if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
 469                  fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
 470                 {
 471                   track2->SetOverlap(k,i);
 472                   track1->SetOverlap(k,j);
 473                 }
 474             }
 475         }
 476     }
 477   cout<<"done"<<endl;
 478 }
 479
 480
 481 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
 482 {
 483
 484   Float_t padw,timew;
 485   if(fPatch < 3)
 486     padw = fTransform->GetPadPitchWidthLow();
 487   else
 488     padw = fTransform->GetPadPitchWidthUp();
 489   Float_t charge = (Float_t)cl->fCharge;
 490   Float_t pad = (Float_t)cl->fPad/charge;
 491   Float_t time = (Float_t)cl->fTime/charge;
 492   Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
 493   sigmaY2 = (s2 + 1./12)*padw*padw;
 494
 495   if(s2 != 0)
 496     {
 497       sigmaY2 = sigmaY2*0.108;
 498       if(fPatch<3)
 499         sigmaY2 = sigmaY2*2.07;
 500     }
 501
 502   s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
 503   timew = fTransform->GetZWidth();
 504   sigmaZ2 = (s2 +1./12)*timew*timew;
 505   if(s2 != 0)
 506     {
 507       sigmaZ2 = sigmaZ2*0.169;
 508       if(fPatch < 3)
 509         sigmaZ2 = sigmaZ2*1.77;
 510     }
 511
 512 }