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