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