HLT/comp/AliL3OfflineDataCompressor.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Transform.h"
   9 #include "AliL3ModelTrack.h"
  10 #include "AliL3Compress.h"
  11 #include "AliL3TrackArray.h"
  12
  13 #include <AliKalmanTrack.h>
  14 #include <AliTPCtracker.h>
  15 #include <AliTPCtrackerMI.h>
  16 #include <AliTPCcluster.h>
  17 #include <AliTPCclusterMI.h>
  18 #include <AliTPCParamSR.h>
  19 #include <AliTPCClustersArray.h>
  20 #include <AliTPCcluster.h>
  21 #include <AliTPCClustersRow.h>
  22 #include <AliTPC.h>
  23 #include <AliTPCv2.h>
  24 #include <AliRun.h>
  25
  26 #include <TTree.h>
  27 #include <TFile.h>
  28 #include <TMath.h>
  29 #include <TDirectory.h>
  30 #include <TSystem.h>
  31 #include <TH2F.h>
  32
  33
  34 #include "AliL3OfflineDataCompressor.h"
  35
  36 #if GCCVERSION == 3
  37 using namespace std;
  38 #endif
  39
  40 //_____________________________________________________________
  41 //
  42 //  AliL3OfflineDataCompression
  43 //
  44
  45
  46 ClassImp(AliL3OfflineDataCompressor)
  47
  48 AliL3OfflineDataCompressor::AliL3OfflineDataCompressor()
  49 {
  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 {
  57   fMarian = MI;
  58   fTracker=0;
  59 }
  60
  61 AliL3OfflineDataCompressor::~AliL3OfflineDataCompressor()
  62 {
  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   fSinglePatch = sp;
  77
  78   char filename[1024];
  79   AliKalmanTrack::SetConvConst(1000/0.299792458/AliL3Transform::GetSolenoidField());
  80   if(fMarian==kFALSE)
  81     sprintf(filename,"%s/offline/AliTPCclusters.root",fPath);
  82   else
  83     sprintf(filename,"%s/offline/AliTPCclustersMI.root",fPath);
  84
  85   TFile *in = TFile::Open(filename);
  86   AliTPCParam *param=(AliTPCParam*)in->Get("75x40_100x60_150x60");
  87
  88   if(fMarian==kFALSE)
  89     fTracker = new AliTPCtracker(param);
  90   else
  91     fTracker = new AliTPCtrackerMI(param);
  92   fTracker->SetEventNumber(event);
  93   //Problems:
  94   //fTracker->LoadClusters();
  95
  96   if(fMarian==kTRUE)
  97     {
  98       AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
  99       mitracker->LoadInnerSectors();
 100       mitracker->LoadOuterSectors();
 101     }
 102
 103   const Int_t MAX=20000;
 104   Int_t nentr=0,i=0; TObjArray tarray(MAX);
 105   if(fMarian==kFALSE)
 106     sprintf(filename,"%s/offline/AliTPCtracks.root",fPath);
 107   else
 108     sprintf(filename,"%s/offline/AliTPCtracksMI.root",fPath);
 109   TFile *tf=TFile::Open(filename);
 110
 111   char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
 112   TTree *tracktree=(TTree*)tf->Get(tname);
 113
 114   TBranch *tbranch=tracktree->GetBranch("tracks");
 115   nentr=(Int_t)tracktree->GetEntries();
 116   AliTPCtrack *iotrack=0;
 117
 118   for (i=0; i<nentr; i++) {
 119     iotrack=new AliTPCtrack;
 120     tbranch->SetAddress(&iotrack);
 121     tracktree->GetEvent(i);
 122     tarray.AddLast(iotrack);
 123   }
 124   delete tracktree;
 125   tf->Close();
 126
 127   AliL3TrackArray *comptracks = new AliL3TrackArray("AliL3ModelTrack");
 128   cout<<"Loaded "<<nentr<<" offline tracks"<<endl;
 129   Int_t slice,padrow;
 130   Int_t totcounter=0;
 131   for(i=0; i<nentr; i++)
 132     {
 133
 134       AliTPCtrack *track=(AliTPCtrack*)tarray.UncheckedAt(i);
 135       Int_t nhits = track->GetNumberOfClusters();
 136       Int_t idx = track->GetClusterIndex(nhits-1);
 137       Int_t sec=(idx&0xff000000)>>24, row=(idx&0x00ff0000)>>16;
 138
 139       /*
 140         TPC sector numbering within the AliTPCtracker class:
 141         There are in total 18 inner sectors and 18 outer sectors.
 142         This means that one sector includes _both_ sides of the TPC.
 143         Example: sec=0 -> sector 0 and 18.
 144                  sec=18 -> sector 36 and 54
 145       */
 146
 147       if(fMarian==kFALSE)
 148         if(sec >= 18)
 149           sec += 18;
 150
 151       AliL3Transform::Sector2Slice(slice,padrow,sec,row);
 152       Double_t par[5],xk=AliL3Transform::Row2X(padrow);
 153       track->PropagateTo(xk);
 154       track->GetExternalParameters(xk,par);
 155       Double_t psi = TMath::ASin(par[2]) + track->GetAlpha();
 156       if (psi<-TMath::Pi()) psi+=2*TMath::Pi();
 157       if (psi>=TMath::Pi()) psi-=2*TMath::Pi();
 158       Float_t pt_1=TMath::Abs(par[4]);
 159       Int_t charge = 1;
 160       if(par[4] > 0)
 161         charge=-1;
 162
 163       Float_t first[3];
 164       AliCluster *fcl=0;
 165       if(fMarian==kFALSE)
 166         fcl= fTracker->GetCluster(idx);
 167       else
 168         {
 169           AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
 170           fcl = mitracker->GetClusterMI(idx);
 171         }
 172       first[0] = xk;
 173       first[1] = fcl->GetY();
 174       first[2] = fcl->GetZ();
 175
 176       AliL3Transform::Local2Global(first,slice);
 177
 178       AliL3ModelTrack *outtrack = (AliL3ModelTrack*)comptracks->NextTrack();
 179       outtrack->SetNHits(nhits);
 180       outtrack->SetFirstPoint(first[0],first[1],first[2]);
 181       outtrack->SetPt(1/pt_1);
 182       outtrack->SetPsi(psi);
 183       outtrack->SetTgl(par[3]);
 184       outtrack->SetCharge(charge);
 185       outtrack->CalculateHelix();
 186       outtrack->Init(0,-1);
 187
 188       //cout<<"Loading track with "<<nhits<<" hits"<<endl;
 189       for(int j=nhits-1; j>=0; j--)
 190         {
 191
 192           Int_t index = track->GetClusterIndex(j);
 193           if(index == 0)
 194             continue;
 195
 196           Float_t xyz[3];
 197           Int_t clustercharge =0;
 198
 199           //AliTPCcluster *cluster = (AliTPCcluster*)tracker->GetCluster(index);
 200           AliCluster *cluster=0;
 201
 202           if(fMarian==kFALSE)
 203             cluster = fTracker->GetCluster(index);
 204           else
 205             {
 206               AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
 207               cluster = mitracker->GetClusterMI(index);
 208             }
 209
 210           xyz[1] = cluster->GetY();
 211           xyz[2] = cluster->GetZ();
 212           if(fMarian==kFALSE)
 213             {
 214               AliTPCcluster *cl = (AliTPCcluster*)cluster;
 215               clustercharge = (Int_t)cl->GetQ();
 216             }
 217           else
 218             {
 219               AliTPCclusterMI *cl = (AliTPCclusterMI*)cluster;
 220               clustercharge = (Int_t)cl->GetQ();
 221             }
 222
 223           cluster->Use();
 224
 225           sec=(index&0xff000000)>>24; row=(index&0x00ff0000)>>16;
 226
 227           if(fMarian==kFALSE)
 228             {
 229               if(sec >= 18)
 230                 sec += 18;
 231
 232               if(xyz[2] < 0)
 233                 sec += 18;
 234             }
 235
 236           //cout<<"sector "<<sec<<" row "<<row<<endl;
 237           if(!AliL3Transform::Sector2Slice(slice,padrow,sec,row))
 238             exit(5);
 239           xyz[0] = AliL3Transform::Row2X(padrow);
 240
 241           //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" index "<<index<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
 242           AliL3Transform::Local2Raw(xyz,sec,row);
 243           //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
 244
 245           if(xyz[1] < -1 || xyz[1] > AliL3Transform::GetNPads(padrow) ||
 246              xyz[2] < -1 || xyz[2] > AliL3Transform::GetNTimeBins())
 247             {
 248               cerr<<"AliL3DataCompressor::FillOfflineData : Wrong time "<<xyz[2]<<" in slice "
 249                   <<slice<<" padrow "<<padrow<<endl;
 250               cout<<"sector "<<sec<<" row "<<row<<endl;
 251               //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
 252               cout<<"Track hit "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 253               exit(5);
 254             }
 255
 256           Float_t angle = 0;
 257           AliL3Transform::Local2GlobalAngle(&angle,slice);
 258           if(!outtrack->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
 259             {
 260               cerr<<"AliL3DataCompressor::FillOfflineData : Error in crossing point calc on slice "
 261                   <<slice<<" row "<<padrow<<endl;
 262               exit(5);
 263             }
 264           Float_t xyz_cross[3] = {outtrack->GetPointX(),outtrack->GetPointY(),outtrack->GetPointZ()};
 265           AliL3Transform::Global2Raw(xyz_cross,sec,row);
 266           /*
 267             if(fabs(xyz_cross[1] - xyz[1]) > 10 ||
 268             fabs(xyz_cross[2] - xyz[2]) > 10)
 269             {
 270             cout<<"AliL3DataCompressor::FillOfflineData : Wrong crossing slice "<<slice<<" padrow "
 271             <<padrow<<" pad "<<xyz[1]<<" padhit "<<xyz_cross[1]<<" time "<<xyz[2]<<" timehit "<<xyz_cross[2]<<endl;
 272             outtrack->Print();
 273             exit(5);
 274             }
 275           */
 276           //cout<<" crossing "<<xyz_cross[0]<<" "<<xyz_cross[1]<<" "<<xyz_cross[2]<<endl;
 277           outtrack->SetPadHit(padrow,xyz_cross[1]);
 278           outtrack->SetTimeHit(padrow,xyz_cross[2]);
 279
 280           if(fWriteClusterShape)
 281             {
 282               Float_t angle = outtrack->GetCrossingAngle(padrow,slice);
 283               outtrack->SetCrossingAngleLUT(padrow,angle);
 284               outtrack->CalculateClusterWidths(padrow,kTRUE);
 285               Int_t patch = AliL3Transform::GetPatch(padrow);
 286               Float_t sigmaY2 = cluster->GetSigmaY2() / pow(AliL3Transform::GetPadPitchWidth(patch),2);
 287               Float_t sigmaZ2 = cluster->GetSigmaZ2() / pow(AliL3Transform::GetZWidth(),2);
 288               outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,sigmaY2,sigmaZ2,3);
 289             }
 290           else
 291             outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,0,0,3);
 292           totcounter++;
 293           outtrack->GetClusterModel(padrow)->fSlice = slice;
 294         }
 295     }
 296
 297   cout<<"AliL3DataCompressor::FillOfflineData : Wrote "<<totcounter<<" clusters"<<endl;
 298   //Write tracks to file
 299   AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 300   comp->WriteFile(comptracks);
 301   delete comp;
 302   delete comptracks;
 303
 304 }
 305
 306 void AliL3OfflineDataCompressor::WriteRemaining(Bool_t select)
 307 {
 308   //Write remaining clusters (not assigned to any tracks) to file
 309
 310   if(!fKeepRemaining)
 311     return;
 312
 313   if(select)
 314     SelectRemainingClusters();
 315
 316   Char_t filename[1024];
 317
 318   if(!fSinglePatch)
 319     {
 320       cerr<<"AliL3OfflineDataCompressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
 321       return;
 322     }
 323
 324   cout<<"Writing remaining clusters "<<endl;
 325   Int_t nrows = AliL3Transform::GetNRows(),sector,row,sec;
 326   AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
 327   for(Int_t slice=0; slice<=35; slice++)
 328     {
 329       sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
 330       FILE *outfile = fopen(filename,"w");
 331       if(!outfile)
 332         {
 333           cerr<<"AliL3OfflineDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
 334           exit(5);
 335         }
 336
 337       for(Int_t padrow=0; padrow < nrows; padrow++)
 338         {
 339           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 340           sec=sector;
 341
 342           if(fMarian == kFALSE)
 343             {
 344               if(slice >= 18)
 345                 sec -= 18;
 346               if(sec >= 18)
 347                 sec -= 18;
 348             }
 349           //cout<<"Getting clusters in sector "<<sec<<" row "<<row<<endl;
 350           // Problems:
 351           Int_t ncl = 0;//tracker->GetNClusters(sec,row);
 352
 353           Int_t counter=0;
 354           Int_t j;
 355           for(j=0; j<ncl; j++)
 356             {
 357               AliTPCcluster *cluster;// = (AliTPCcluster*)tracker->GetCluster(sec,row,j);
 358               if(cluster->GetZ() < 0 && slice < 18) continue;
 359               if(cluster->GetZ() > 0 && slice > 17) continue;
 360               if(cluster->IsUsed())
 361                 continue;
 362               counter++;
 363             }
 364
 365           Int_t size = sizeof(AliL3RemainingRow) + counter*sizeof(AliL3RemainingCluster);
 366           Byte_t *data = new Byte_t[size];
 367           AliL3RemainingRow *tempPt = (AliL3RemainingRow*)data;
 368
 369           tempPt->fPadRow = padrow;
 370           tempPt->fNClusters = counter;
 371           //cout<<"Found "<<counter<<" unused out of "<<ncl<<" clusters on slice "<<slice<<" padrow "<<padrow<<endl;
 372           Int_t local_counter=0;
 373           for(j=0; j<ncl; j++)
 374             {
 375               //Problems:
 376               AliTPCcluster *cluster;// = (AliTPCcluster*)tracker->GetCluster(sec,row,j);
 377               if(cluster->GetZ() < 0 && slice < 18) continue;
 378               if(cluster->GetZ() > 0 && slice > 17) continue;
 379               if(cluster->IsUsed())
 380                 continue;
 381
 382               if(local_counter > counter)
 383                 {
 384                   cerr<<"AliL3OfflineDataCompressor::WriterRemaining : array out of range "<<local_counter<<" "<<counter<<endl;
 385                   return;
 386                 }
 387               Float_t xyz[3] = {AliL3Transform::Row2X(padrow),cluster->GetY(),cluster->GetZ()};
 388               AliL3Transform::Local2Raw(xyz,sector,row);
 389               //cout<<"y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
 390               tempPt->fClusters[local_counter].fY = cluster->GetY();
 391               tempPt->fClusters[local_counter].fZ = cluster->GetZ();
 392               tempPt->fClusters[local_counter].fCharge = (UShort_t)cluster->GetQ();
 393               tempPt->fClusters[local_counter].fSigmaY2 = cluster->GetSigmaY2();
 394               tempPt->fClusters[local_counter].fSigmaZ2 = cluster->GetSigmaZ2();
 395               local_counter++;
 396             }
 397
 398           fwrite(tempPt,size,1,outfile);
 399           delete [] data;
 400         }
 401       fclose(outfile);
 402     }
 403 }
 404
 405 void AliL3OfflineDataCompressor::SelectRemainingClusters()
 406 {
 407
 408   cout<<"Cleaning up clusters"<<endl;
 409   Int_t nrows = AliL3Transform::GetNRows();
 410   Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
 411
 412   Int_t sector,row,sec;
 413   AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
 414   for(Int_t slice=0; slice<36; slice++)
 415     {
 416       for(Int_t padrow=0; padrow < nrows; padrow++)
 417         {
 418           if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
 419
 420           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 421           sec=sector;
 422
 423           if(fMarian == kFALSE)
 424             {
 425               if(slice >= 18)
 426                 sec -= 18;
 427               if(sec >= 18)
 428                 sec -= 18;
 429             }
 430           //Problems:
 431           Int_t ncl = 0;//tracker->GetNClusters(sec,row);
 432           for(Int_t j=0; j<ncl; j++)
 433             {
 434               AliTPCcluster *cluster;// = (AliTPCcluster*)tracker->GetCluster(sec,row,j);
 435               if(cluster->GetZ() < 0 && slice < 18) continue;
 436               if(cluster->GetZ() > 0 && slice > 17) continue;
 437               if(cluster->IsUsed())
 438                 continue;
 439
 440               cluster->Use();
 441             }
 442         }
 443
 444     }
 445 }