HLT/comp/AliL3OfflineDataCompressor.cxx

   1 // @(#) $Id$
   2
   3 //_____________________________________________________________
   4 //
   5 //  AliL3OfflineDataCompression
   6 //
   7 // Class to compress data with offline tracks
   8 // as seeds.
   9 //
  10 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
  11 //*-- Copyright &copy ALICE HLT Group
  12
  13 #include "AliL3StandardIncludes.h"
  14
  15 #include <AliTPCParamSR.h>
  16 #include <AliTPCClustersArray.h>
  17 #include <AliTPCcluster.h>
  18 #include <AliTPCClustersRow.h>
  19 #include <AliTPC.h>
  20 #include <AliTPCv2.h>
  21 #include <AliTPCcluster.h>
  22 #include <AliTPCtracker.h>
  23 #include <AliTPCclusterMI.h>
  24 #include <AliTPCtrackerMI.h>
  25 #include <AliKalmanTrack.h>
  26 #include <AliRun.h>
  27
  28 #include <TTree.h>
  29 #include <TFile.h>
  30 #include <TMath.h>
  31 #include <TDirectory.h>
  32 #include <TSystem.h>
  33 #include <TH1F.h>
  34
  35 #include "AliL3Transform.h"
  36 #include "AliL3ModelTrack.h"
  37 #include "AliL3Compress.h"
  38 #include "AliL3TrackArray.h"
  39 #include "bitio.h"
  40 #include "AliL3OfflineDataCompressor.h"
  41
  42 #if __GNUC__ == 3
  43 using namespace std;
  44 #endif
  45
  46 ClassImp(AliL3OfflineDataCompressor)
  47
  48 AliL3OfflineDataCompressor::AliL3OfflineDataCompressor()
  49 { //constructor
  50   fMarian = kFALSE;
  51   fTracker=0;
  52 }
  53
  54 AliL3OfflineDataCompressor::AliL3OfflineDataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape,Bool_t MI)
  55   : AliL3DataCompressor(path,keep,writeshape)
  56 { //constructor
  57   fMarian = MI;
  58   fTracker=0;
  59 }
  60
  61 AliL3OfflineDataCompressor::~AliL3OfflineDataCompressor()
  62 { //deconstructor
  63   if(fTracker)
  64     {
  65       fTracker->UnloadClusters();
  66       delete fTracker;
  67     }
  68 }
  69
  70
  71 void AliL3OfflineDataCompressor::LoadData(Int_t event,Bool_t sp)
  72 {
  73   //Take offline reconstructed tracks as an input.
  74   //In this case, no remaining clusters are written.
  75
  76   if(fTracker)
  77     {
  78       fTracker->UnloadClusters();
  79       delete fTracker;
  80     }
  81
  82   fSinglePatch = sp;
  83   fEvent = event;
  84
  85   char filename[1024];
  86   AliKalmanTrack::SetConvConst(1000/0.299792458/AliL3Transform::GetSolenoidField());
  87   if(fMarian==kFALSE)
  88     sprintf(filename,"%s/offline/AliTPCclusters.root",fPath);
  89   else
  90     sprintf(filename,"%s/offline/AliTPCclustersMI.root",fPath);
  91
  92   bool compressed=0;
  93   if(compressed)
  94     {
  95       cout<<"AliL3OfflineDataCompressor::LoadData : Taking compressed offline files!!"<<endl;
  96       sprintf(filename,"%s/comp/offline/AliTPCclusters.root",fPath);
  97     }
  98   TFile *in = TFile::Open(filename);
  99   AliTPCParam *param=(AliTPCParam*)in->Get("75x40_100x60_150x60");
 100
 101   if(fMarian==kFALSE)
 102     fTracker = new AliTPCtracker(param);
 103   else
 104     fTracker = new AliTPCtrackerMI(param);
 105   fTracker->SetEventNumber(event);
 106 #ifdef asvversion
 107   fTracker->LoadClusters();
 108 #endif
 109   if(fMarian==kTRUE)
 110     {
 111 #ifdef asvversion
 112       AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
 113       mitracker->LoadInnerSectors();
 114       mitracker->LoadOuterSectors();
 115 #endif
 116     }
 117
 118   const Int_t kMAX=20000;
 119   Int_t nentr=0,i=0; TObjArray tarray(kMAX);
 120   if(fMarian==kFALSE)
 121     sprintf(filename,"%s/offline/AliTPCtracks.root",fPath);
 122   else
 123     sprintf(filename,"%s/offline/AliTPCtracksMI.root",fPath);
 124
 125   if(compressed)
 126     sprintf(filename,"%s/comp/offline/AliTPCtracks.root",fPath);
 127
 128   TFile *tf=TFile::Open(filename);
 129
 130   char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
 131   TTree *tracktree=(TTree*)tf->Get(tname);
 132
 133   TBranch *tbranch=tracktree->GetBranch("tracks");
 134   nentr=(Int_t)tracktree->GetEntries();
 135   AliTPCtrack *iotrack=0;
 136
 137   for (i=0; i<nentr; i++) {
 138     iotrack=new AliTPCtrack;
 139     tbranch->SetAddress(&iotrack);
 140     tracktree->GetEvent(i);
 141     tarray.AddLast(iotrack);
 142   }
 143   delete tracktree;
 144   tf->Close();
 145
 146   AliL3TrackArray *comptracks = new AliL3TrackArray("AliL3ModelTrack");
 147   cout<<"Loaded "<<nentr<<" offline tracks"<<endl;
 148   Int_t slice,padrow;
 149   Int_t totcounter=0;
 150
 151   for(i=0; i<nentr; i++)
 152     {
 153
 154       AliTPCtrack *track=(AliTPCtrack*)tarray.UncheckedAt(i);
 155       Int_t nhits = track->GetNumberOfClusters();
 156       Int_t idx = track->GetClusterIndex(nhits-1);
 157       Int_t sec=(idx&0xff000000)>>24, row=(idx&0x00ff0000)>>16;
 158
 159       /*
 160         TPC sector numbering within the AliTPCtracker class:
 161         There are in total 18 inner sectors and 18 outer sectors.
 162         This means that one sector includes _both_ sides of the TPC.
 163         Example: sec=0 -> sector 0 and 18.
 164                  sec=18 -> sector 36 and 54
 165       */
 166
 167       if(fMarian==kFALSE)
 168         if(sec >= 18)
 169           sec += 18;
 170
 171       AliL3Transform::Sector2Slice(slice,padrow,sec,row);
 172       Double_t par[5],xk=AliL3Transform::Row2X(padrow);
 173       track->PropagateTo(xk);
 174       track->GetExternalParameters(xk,par);
 175       Double_t psi = TMath::ASin(par[2]) + track->GetAlpha();
 176       if (psi<-TMath::Pi()) psi+=2*TMath::Pi();
 177       if (psi>=TMath::Pi()) psi-=2*TMath::Pi();
 178       Float_t pt1=TMath::Abs(par[4]);
 179       Int_t charge = 1;
 180       if(par[4] > 0)
 181         charge=-1;
 182
 183       Float_t first[3];
 184       AliCluster *fcl=0;
 185       if(fMarian==kFALSE)
 186         fcl= fTracker->GetCluster(idx);
 187       else
 188         {
 189           AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
 190           fcl = mitracker->GetClusterMI(idx);
 191         }
 192       first[0] = xk;
 193       first[1] = par[0];
 194       first[2] = par[1];
 195
 196       AliL3Transform::Local2Global(first,slice);
 197
 198       AliL3ModelTrack *outtrack = (AliL3ModelTrack*)comptracks->NextTrack();
 199       outtrack->SetNHits(nhits);
 200       outtrack->SetFirstPoint(first[0],first[1],first[2]);
 201       outtrack->SetPt(1./pt1);
 202       outtrack->SetPsi(psi);
 203       outtrack->SetTgl(par[3]);
 204       outtrack->SetCharge(charge);
 205       outtrack->CalculateHelix();
 206       outtrack->Init(0,-1);
 207
 208       //cout<<"Loading track with "<<nhits<<" hits"<<endl;
 209       for(int j=nhits-1; j>=0; j--)
 210         {
 211
 212           Int_t index = track->GetClusterIndex(j);
 213           if(index == 0)
 214             continue;
 215
 216           Float_t xyz[3];
 217           Int_t clustercharge =0;
 218
 219           //AliTPCcluster *cluster = (AliTPCcluster*)tracker->GetCluster(index);
 220           AliCluster *cluster=0;
 221
 222           if(fMarian==kFALSE)
 223             cluster = fTracker->GetCluster(index);
 224           else
 225             {
 226               AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
 227               cluster = mitracker->GetClusterMI(index);
 228             }
 229
 230           xyz[1] = cluster->GetY();
 231           xyz[2] = cluster->GetZ();
 232           if(fMarian==kFALSE)
 233             {
 234               AliTPCcluster *cl = (AliTPCcluster*)cluster;
 235               clustercharge = (Int_t)cl->GetQ();
 236             }
 237           else
 238             {
 239               AliTPCclusterMI *cl = (AliTPCclusterMI*)cluster;
 240               clustercharge = (Int_t)cl->GetQ();
 241             }
 242
 243           cluster->Use();
 244
 245           sec=(index&0xff000000)>>24; row=(index&0x00ff0000)>>16;
 246
 247           if(fMarian==kFALSE)
 248             {
 249               if(sec >= 18)
 250                 sec += 18;
 251
 252               if(xyz[2] < 0)
 253                 sec += 18;
 254             }
 255
 256           //cout<<"sector "<<sec<<" row "<<row<<endl;
 257           if(!AliL3Transform::Sector2Slice(slice,padrow,sec,row))
 258             exit(5);
 259           xyz[0] = AliL3Transform::Row2X(padrow);
 260
 261           //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" index "<<index<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
 262           AliL3Transform::Local2Raw(xyz,sec,row);
 263           //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
 264
 265           if(xyz[1] < -1 || xyz[1] > AliL3Transform::GetNPads(padrow) ||
 266              xyz[2] < -1 || xyz[2] > AliL3Transform::GetNTimeBins())
 267             {
 268               cerr<<"AliL3DataCompressor::FillOfflineData : Wrong time "<<xyz[2]<<" in slice "
 269                   <<slice<<" padrow "<<padrow<<endl;
 270               cout<<"sector "<<sec<<" row "<<row<<endl;
 271               //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
 272               cout<<"Track hit "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 273               exit(5);
 274             }
 275
 276           Float_t angle = 0;
 277           AliL3Transform::Local2GlobalAngle(&angle,slice);
 278           if(!outtrack->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
 279             {
 280               cerr<<"AliL3DataCompressor::FillOfflineData : Error in crossing point calc on slice "
 281                   <<slice<<" row "<<padrow<<endl;
 282               exit(5);
 283             }
 284
 285           Float_t xyzcross[3] = {outtrack->GetPointX(),outtrack->GetPointY(),outtrack->GetPointZ()};
 286           AliL3Transform::Global2Raw(xyzcross,sec,row);
 287           /*
 288             if(fabs(xyzcross[1] - xyz[1]) > 10 ||
 289             fabs(xyzcross[2] - xyz[2]) > 10)
 290             {
 291             cout<<"AliL3DataCompressor::FillOfflineData : Wrong crossing slice "<<slice<<" padrow "
 292             <<padrow<<" pad "<<xyz[1]<<" padhit "<<xyzcross[1]<<" time "<<xyz[2]<<" timehit "<<xyzcross[2]<<endl;
 293             outtrack->Print();
 294             exit(5);
 295             }
 296           */
 297           //cout<<" crossing "<<xyzcross[0]<<" "<<xyzcross[1]<<" "<<xyzcross[2]<<endl;
 298           outtrack->SetPadHit(padrow,xyzcross[1]);
 299           outtrack->SetTimeHit(padrow,xyzcross[2]);
 300
 301           if(fWriteClusterShape)
 302             {
 303               Float_t angle = outtrack->GetCrossingAngle(padrow,slice);
 304               outtrack->SetCrossingAngleLUT(padrow,angle);
 305               outtrack->CalculateClusterWidths(padrow,kTRUE);
 306               Int_t patch = AliL3Transform::GetPatch(padrow);
 307               Float_t sigmaY2 = cluster->GetSigmaY2() / pow(AliL3Transform::GetPadPitchWidth(patch),2);
 308               Float_t sigmaZ2 = cluster->GetSigmaZ2() / pow(AliL3Transform::GetZWidth(),2);
 309               outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,sigmaY2,sigmaZ2,3);
 310             }
 311           else
 312             outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,0,0,3);
 313           totcounter++;
 314           outtrack->GetClusterModel(padrow)->fSlice = slice;
 315           fNusedClusters++;
 316         }
 317     }
 318
 319   cout<<"AliL3DataCompressor::FillOfflineData : Wrote "<<totcounter<<" clusters"<<endl;
 320   //Write tracks to file
 321   AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 322   comp->WriteFile(comptracks);
 323   delete comp;
 324   delete comptracks;
 325 }
 326
 327 void AliL3OfflineDataCompressor::WriteRemaining(Bool_t select)
 328 {
 329   //Write remaining clusters (not assigned to any tracks) to file
 330
 331   if(!fKeepRemaining)
 332     return;
 333
 334   if(select)
 335     SelectRemainingClusters();
 336
 337   Char_t filename[1024];
 338
 339   if(!fSinglePatch)
 340     {
 341       cerr<<"AliL3OfflineDataCompressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
 342       return;
 343     }
 344
 345   ofstream idfile;
 346   if(fWriteIdsToFile)
 347     {
 348       sprintf(filename,"%s/comp/remains_ids.txt",fPath);
 349       idfile.open(filename);
 350     }
 351
 352   cout<<"Writing remaining clusters "<<endl;
 353   Int_t nrows = AliL3Transform::GetNRows(),sector,row,sec;
 354 #ifdef asvversion
 355   AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
 356 #endif
 357   for(Int_t slice=0; slice<=35; slice++)
 358     {
 359       sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
 360       BIT_FILE *output = OpenOutputBitFile(filename);
 361       if(!output)
 362         {
 363           cerr<<"AliL3OfflineDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
 364           exit(5);
 365         }
 366
 367       //Write number of padrows with clusters
 368       OutputBits(output,nrows,8);
 369
 370       for(Int_t padrow=0; padrow < nrows; padrow++)
 371         {
 372           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 373           sec=sector;
 374
 375           if(fMarian == kFALSE)
 376             {
 377               if(slice >= 18)
 378                 sec -= 18;
 379               if(sec >= 18)
 380                 sec -= 18;
 381             }
 382           //cout<<"Getting clusters in sector "<<sec<<" row "<<row<<endl;
 383           Int_t ncl = 0;
 384 #ifdef asvversion
 385           ncl = tracker->GetNClusters(sec,row);
 386 #endif
 387
 388           Int_t counter=0;
 389           Int_t j;
 390           for(j=0; j<ncl; j++)
 391             {
 392               AliTPCcluster *cluster = 0;
 393 #ifdef asvversion
 394               cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
 395 #endif
 396               if(cluster->GetZ() < 0 && slice < 18) continue;
 397               if(cluster->GetZ() > 0 && slice > 17) continue;
 398               if(cluster->IsUsed())
 399                 continue;
 400               counter++;
 401             }
 402
 403           OutputBits(output,padrow,8);//Write padrow #
 404           OutputBits(output,counter,10);//Write number of clusters on this padrow
 405
 406           //cout<<"Found "<<counter<<" unused out of "<<ncl<<" clusters on slice "<<slice<<" padrow "<<padrow<<endl;
 407           for(j=0; j<ncl; j++)
 408             {
 409               AliTPCcluster *cluster = 0;
 410 #ifdef asvversion
 411               cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
 412 #endif
 413               if(cluster->GetZ() < 0 && slice < 18) continue;
 414               if(cluster->GetZ() > 0 && slice > 17) continue;
 415               if(cluster->IsUsed())
 416                 continue;
 417
 418               if(fWriteIdsToFile)
 419                 idfile << cluster->GetLabel(0)<<' ';
 420
 421               Float_t xyz[3] = {AliL3Transform::Row2X(padrow),cluster->GetY(),cluster->GetZ()};
 422               AliL3Transform::Local2Raw(xyz,sector,row);
 423
 424               Int_t patch = AliL3Transform::GetPatch(padrow);
 425               Float_t padw = cluster->GetSigmaY2()/pow(AliL3Transform::GetPadPitchWidth(patch),2);
 426               Float_t timew = cluster->GetSigmaZ2()/pow(AliL3Transform::GetZWidth(),2);
 427
 428               Int_t buff;
 429               //Write pad
 430               buff = (Int_t)rint(xyz[1]*10);
 431               if(buff<0)
 432                 {
 433                   cerr<<"AliL3OfflineDataCompressor:WriteRemaining : Wrong pad value "<<buff<<endl;
 434                   exit(5);
 435                 }
 436               OutputBits(output,buff,11);
 437
 438               //Write time
 439               buff = (Int_t)rint(xyz[2]*10);
 440               if(buff<0)
 441                 {
 442                   cerr<<"AliL3OfflineDataCompressor:WriteRemaining : Wrong time value "<<buff<<endl;
 443                   exit(5);
 444                 }
 445               OutputBits(output,buff,13);
 446
 447               //Write widths
 448               buff = (Int_t)rint(padw*100);
 449               OutputBits(output,buff,8);
 450               buff = (Int_t)rint(timew*100);
 451               OutputBits(output,buff,8);
 452
 453               //Write charge
 454               buff = (Int_t)cluster->GetQ();
 455               if(buff >= 1<<14)
 456                 buff = (1<<14)-1;
 457               OutputBits(output,buff,14);
 458
 459               fNunusedClusters++;
 460             }
 461         }
 462       CloseOutputBitFile(output);
 463     }
 464   if(fWriteIdsToFile)
 465     {
 466       idfile << endl;
 467       idfile.close();
 468     }
 469 }
 470
 471 void AliL3OfflineDataCompressor::SelectRemainingClusters()
 472 {
 473   //select the remaining clusters
 474   //which were not compressed
 475   cout<<"Cleaning up clusters"<<endl;
 476   Int_t nrows = AliL3Transform::GetNRows();
 477   Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
 478
 479   Int_t sector,row,sec;
 480
 481 #ifdef asvversion
 482   AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
 483 #endif
 484   for(Int_t slice=0; slice<36; slice++)
 485     {
 486       for(Int_t padrow=0; padrow < nrows; padrow++)
 487         {
 488
 489           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 490           sec=sector;
 491
 492           if(fMarian == kFALSE)
 493             {
 494               if(slice >= 18)
 495                 sec -= 18;
 496               if(sec >= 18)
 497                 sec -= 18;
 498             }
 499           Int_t ncl = 0;
 500 #ifdef asvversion
 501           ncl=tracker->GetNClusters(sec,row);
 502 #endif
 503           for(Int_t j=0; j<ncl; j++)
 504             {
 505               AliTPCcluster *cluster = 0;
 506 #ifdef asvversion
 507               cluster = (AliTPCcluster*)tracker->GetCluster(sec,row,j);
 508 #endif
 509               if(cluster->IsUsed())
 510                 continue;
 511
 512               //Check the widths (errors) of the cluster, and remove big bastards:
 513               Float_t xyw = cluster->GetSigmaY2() / pow(AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(padrow)),2);
 514               Float_t zw  = cluster->GetSigmaZ2() / pow(AliL3Transform::GetZWidth(),2);
 515               if(xyw >= 2.55 || zw >= 2.55)//Because we use 1 byte to store
 516                 {
 517                   cluster->Use();
 518                   continue;
 519                 }
 520
 521               //if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
 522
 523               if(padrow == nrows - 1 || padrow == nrows - 1 - gap ||                 //First seeding
 524                  padrow == nrows - 1 - shift || padrow == nrows - 1 - gap - shift)   //Second seeding
 525                 continue;
 526
 527               if(cluster->GetZ() < 0 && slice < 18) continue;
 528               if(cluster->GetZ() > 0 && slice > 17) continue;
 529               if(cluster->IsUsed())
 530                 continue;
 531
 532               cluster->Use();
 533             }
 534         }
 535
 536     }
 537 }