[u/mrichter/AliRoot.git] / HLT / comp / AliHLTOfflineDataCompressor.cxx

// @(#) $Id$

//_____________________________________________________________
//
//  AliHLTOfflineDataCompression
//
// Class to compress data with offline tracks
// as seeds.
//
// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ALICE HLT Group

#include "AliHLTStandardIncludes.h"

#include <AliTPCParamSR.h>
#include <AliTPCClustersArray.h>
#include <AliTPCcluster.h>
#include <AliTPCClustersRow.h>
#include <AliTPC.h>
#include <AliTPCv2.h>
#include <AliTPCcluster.h>
#include <AliTPCtracker.h>
#include <AliTPCclusterMI.h>
#include <AliTPCtrackerMI.h>
#include <AliKalmanTrack.h>
#include <AliRun.h>

#include <TTree.h>
#include <TFile.h>
#include <TMath.h>
#include <TDirectory.h>
#include <TSystem.h>
#include <TH1F.h>

#include "AliHLTTransform.h"
#include "AliHLTModelTrack.h"
#include "AliHLTCompress.h"
#include "AliHLTTrackArray.h"
#include "bitio.h"
#include "AliHLTOfflineDataCompressor.h"

#if __GNUC__ == 3
using namespace std;
#endif

ClassImp(AliHLTOfflineDataCompressor)

AliHLTOfflineDataCompressor::AliHLTOfflineDataCompressor()
{ //constructor
  fMarian = kFALSE;
  fTracker=0;
}

AliHLTOfflineDataCompressor::AliHLTOfflineDataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape,Bool_t MI) 
  : AliHLTDataCompressor(path,keep,writeshape)
{ //constructor
  fMarian = MI;
  fTracker=0;
}

AliHLTOfflineDataCompressor::~AliHLTOfflineDataCompressor()
{ //deconstructor
  if(fTracker)
    {
      fTracker->UnloadClusters();
      delete fTracker;
    }
}


void AliHLTOfflineDataCompressor::LoadData(Int_t event,Bool_t sp)
{
  //Take offline reconstructed tracks as an input.
  //In this case, no remaining clusters are written.
  
  if(fTracker)
    {
      fTracker->UnloadClusters();
      delete fTracker;
    }
  
  fSinglePatch = sp;
  fEvent = event;

  char filename[1024];
  //AliKalmanTrack::SetConvConst(1000/0.299792458/AliHLTTransform::GetSolenoidField());
  if(fMarian==kFALSE)
    sprintf(filename,"%s/offline/AliTPCclusters.root",fPath);
  else
    sprintf(filename,"%s/offline/AliTPCclustersMI.root",fPath);
  
  bool compressed=0;
  if(compressed)
    {
      cout<<"AliHLTOfflineDataCompressor::LoadData : Taking compressed offline files!!"<<endl;
      sprintf(filename,"%s/comp/offline/AliTPCclusters.root",fPath);
    }
  TFile *in = TFile::Open(filename);
  AliTPCParam *param=(AliTPCParam*)in->Get("75x40_100x60_150x60");

  if(fMarian==kFALSE)
    fTracker = new AliTPCtracker(param);
  else
    fTracker = new AliTPCtrackerMI(param);
  //fTracker->SetEventNumber(event);
#ifdef asvversion
  fTracker->LoadClusters();
#endif  
  if(fMarian==kTRUE)
    {
#ifdef asvversion
      AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
      mitracker->LoadInnerSectors();
      mitracker->LoadOuterSectors();
#endif
    }
  
  const Int_t kMAX=20000;
  Int_t nentr=0,i=0; TObjArray tarray(kMAX);
  if(fMarian==kFALSE)
    sprintf(filename,"%s/offline/AliTPCtracks.root",fPath);
  else
    sprintf(filename,"%s/offline/AliTPCtracksMI.root",fPath);
  
  if(compressed)
    sprintf(filename,"%s/comp/offline/AliTPCtracks.root",fPath);
  
  TFile *tf=TFile::Open(filename);
  
  char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
  TTree *tracktree=(TTree*)tf->Get(tname);
  
  TBranch *tbranch=tracktree->GetBranch("tracks");
  nentr=(Int_t)tracktree->GetEntries();
  AliTPCtrack *iotrack=0;

  for (i=0; i<nentr; i++) {
    iotrack=new AliTPCtrack;
    tbranch->SetAddress(&iotrack);
    tracktree->GetEvent(i);
    tarray.AddLast(iotrack);
  }   
  delete tracktree; 
  tf->Close();
  
  AliHLTTrackArray *comptracks = new AliHLTTrackArray("AliHLTModelTrack");
  cout<<"Loaded "<<nentr<<" offline tracks"<<endl;
  Int_t slice,padrow;
  Int_t totcounter=0;

  for(i=0; i<nentr; i++)
    {
      
      AliTPCtrack *track=(AliTPCtrack*)tarray.UncheckedAt(i);
      Int_t nhits = track->GetNumberOfClusters();
      Int_t idx = track->GetClusterIndex(nhits-1);
      Int_t sec=(idx&0xff000000)>>24, row=(idx&0x00ff0000)>>16;
      
      /*
	TPC sector numbering within the AliTPCtracker class:
	There are in total 18 inner sectors and 18 outer sectors.
	This means that one sector includes _both_ sides of the TPC.
	Example: sec=0 -> sector 0 and 18.
	         sec=18 -> sector 36 and 54
      */
      
      if(fMarian==kFALSE)
	if(sec >= 18)
	  sec += 18;
      
      AliHLTTransform::Sector2Slice(slice,padrow,sec,row);
      Double_t par[5],xk=AliHLTTransform::Row2X(padrow);
      track->PropagateTo(xk);
      track->GetExternalParameters(xk,par);
      Double_t psi = TMath::ASin(par[2]) + track->GetAlpha();
      if (psi<-TMath::Pi()) psi+=2*TMath::Pi();
      if (psi>=TMath::Pi()) psi-=2*TMath::Pi();
      Float_t pt1=TMath::Abs(par[4]);
      Int_t charge = 1;
      if(par[4] > 0)
	charge=-1;

      Float_t first[3];
      AliCluster *fcl=0;
      if(fMarian==kFALSE)
	fcl= fTracker->GetCluster(idx);
      else
	{
	  AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
	  fcl = mitracker->GetClusterMI(idx);
	}
      first[0] = xk;
      first[1] = par[0];
      first[2] = par[1];

      AliHLTTransform::Local2Global(first,slice);
      
      AliHLTModelTrack *outtrack = (AliHLTModelTrack*)comptracks->NextTrack();
      outtrack->SetNHits(nhits);
      outtrack->SetFirstPoint(first[0],first[1],first[2]);
      outtrack->SetPt(1./pt1);
      outtrack->SetPsi(psi);
      outtrack->SetTgl(par[3]);
      outtrack->SetCharge(charge);
      outtrack->CalculateHelix();
      outtrack->Init(0,-1);
      
      //cout<<"Loading track with "<<nhits<<" hits"<<endl;
      for(int j=nhits-1; j>=0; j--)
	{

	  Int_t index = track->GetClusterIndex(j);
	  if(index == 0)
	    continue;

	  Float_t xyz[3];
	  Int_t clustercharge =0;

	  //AliTPCcluster *cluster = (AliTPCcluster*)tracker->GetCluster(index);
	  AliCluster *cluster=0;
	  
	  if(fMarian==kFALSE)
	    cluster = fTracker->GetCluster(index);
	  else
	    {
	      AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
	      cluster = mitracker->GetClusterMI(index);
	    }

	  xyz[1] = cluster->GetY();
	  xyz[2] = cluster->GetZ();
	  if(fMarian==kFALSE)
	    {
	      AliTPCcluster *cl = (AliTPCcluster*)cluster;
	      clustercharge = (Int_t)cl->GetQ();
	    }
	  else
	    {
	      AliTPCclusterMI *cl = (AliTPCclusterMI*)cluster;
	      clustercharge = (Int_t)cl->GetQ();
	    }

	  cluster->Use();
	  
	  sec=(index&0xff000000)>>24; row=(index&0x00ff0000)>>16;
	  
	  if(fMarian==kFALSE)
	    {
	      if(sec >= 18)
		sec += 18;
	      
	      if(xyz[2] < 0)
		sec += 18;
	    }
	    
	  //cout<<"sector "<<sec<<" row "<<row<<endl;
	  if(!AliHLTTransform::Sector2Slice(slice,padrow,sec,row))
	    exit(5);
	  xyz[0] = AliHLTTransform::Row2X(padrow);
	  
	  //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" index "<<index<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
	  AliHLTTransform::Local2Raw(xyz,sec,row);
	  //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
	  
	  if(xyz[1] < -1 || xyz[1] > AliHLTTransform::GetNPads(padrow) ||
	     xyz[2] < -1 || xyz[2] > AliHLTTransform::GetNTimeBins())
	    {
	      cerr<<"AliHLTDataCompressor::FillOfflineData : Wrong time "<<xyz[2]<<" in slice "
		  <<slice<<" padrow "<<padrow<<endl;
	      cout<<"sector "<<sec<<" row "<<row<<endl;
	      //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
	      cout<<"Track hit "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
	      exit(5);
	    }

	  Float_t angle = 0;
	  AliHLTTransform::Local2GlobalAngle(&angle,slice);
	  if(!outtrack->CalculateReferencePoint(angle,AliHLTTransform::Row2X(padrow)))
	    {
	      cerr<<"AliHLTDataCompressor::FillOfflineData : Error in crossing point calc on slice "
		  <<slice<<" row "<<padrow<<endl;
	      exit(5);
	    }

	  Float_t xyzcross[3] = {outtrack->GetPointX(),outtrack->GetPointY(),outtrack->GetPointZ()};
	  AliHLTTransform::Global2Raw(xyzcross,sec,row);
	  /*
	    if(fabs(xyzcross[1] - xyz[1]) > 10 ||
	    fabs(xyzcross[2] - xyz[2]) > 10)
	    {
	    cout<<"AliHLTDataCompressor::FillOfflineData : Wrong crossing slice "<<slice<<" padrow "
	    <<padrow<<" pad "<<xyz[1]<<" padhit "<<xyzcross[1]<<" time "<<xyz[2]<<" timehit "<<xyzcross[2]<<endl;
	    outtrack->Print();
	    exit(5);
	    }
	  */
	  //cout<<" crossing "<<xyzcross[0]<<" "<<xyzcross[1]<<" "<<xyzcross[2]<<endl;
	  outtrack->SetPadHit(padrow,xyzcross[1]);
	  outtrack->SetTimeHit(padrow,xyzcross[2]);
	  
	  if(fWriteClusterShape)
	    {
	      Float_t angle = outtrack->GetCrossingAngle(padrow,slice);
	      outtrack->SetCrossingAngleLUT(padrow,angle);
	      outtrack->CalculateClusterWidths(padrow,kTRUE);
	      Int_t patch = AliHLTTransform::GetPatch(padrow);
	      Float_t sigmaY2 = cluster->GetSigmaY2() / pow(AliHLTTransform::GetPadPitchWidth(patch),2);
	      Float_t sigmaZ2 = cluster->GetSigmaZ2() / pow(AliHLTTransform::GetZWidth(),2);
	      outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,sigmaY2,sigmaZ2,3);
	    }
	  else
	    outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,0,0,3);
	  totcounter++;
	  outtrack->GetClusterModel(padrow)->fSlice = slice;
	  fNusedClusters++;
	}
    }
  
  cout<<"AliHLTDataCompressor::FillOfflineData : Wrote "<<totcounter<<" clusters"<<endl;
  //Write tracks to file
  AliHLTCompress *comp = new AliHLTCompress(-1,-1,fPath,fWriteClusterShape,fEvent);
  comp->WriteFile(comptracks);
  delete comp;
  delete comptracks;
}

void AliHLTOfflineDataCompressor::WriteRemaining(Bool_t select)
{
  //Write remaining clusters (not assigned to any tracks) to file
  
  if(!fKeepRemaining)
    return;
  
  if(select)
    SelectRemainingClusters();
  
  Char_t filename[1024];
  
  if(!fSinglePatch)
    {
      cerr<<"AliHLTOfflineDataCompressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
      return;
    }
  
  ofstream idfile;
  if(fWriteIdsToFile)
    {
      sprintf(filename,"%s/comp/remains_ids.txt",fPath);
      idfile.open(filename);
    }
  
  cout<<"Writing remaining clusters "<<endl;
  Int_t nrows = AliHLTTransform::GetNRows(),sector,row,sec;
#ifdef asvversion
  AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
#endif
  for(Int_t slice=0; slice<=35; slice++)
    {
      sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
      BIT_FILE *output = OpenOutputBitFile(filename);
      if(!output)
	{
	  cerr<<"AliHLTOfflineDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
	  exit(5);
	}
      
      //Write number of padrows with clusters
      OutputBits(output,nrows,8);
      
      for(Int_t padrow=0; padrow < nrows; padrow++)
	{
	  AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
	  sec=sector;
	  
	  if(fMarian == kFALSE)
	    {
	      if(slice >= 18)
		sec -= 18;
	      if(sec >= 18)
		sec -= 18;
	    }
	  //cout<<"Getting clusters in sector "<<sec<<" row "<<row<<endl;
	  Int_t ncl = 0;
#ifdef asvversion
	  ncl = tracker->GetNClusters(sec,row);
#endif
	  	  
	  Int_t counter=0;
	  Int_t j;
	  for(j=0; j<ncl; j++)
	    {
	      AliTPCcluster *cluster = 0;
#ifdef asvversion
	      cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
#endif
	      if(cluster->GetZ() < 0 && slice < 18) continue;
	      if(cluster->GetZ() > 0 && slice > 17) continue;
	      if(cluster->IsUsed())
		continue;
	      counter++;
	    }
	  
	  OutputBits(output,padrow,8);//Write padrow #
	  OutputBits(output,counter,10);//Write number of clusters on this padrow

	  //cout<<"Found "<<counter<<" unused out of "<<ncl<<" clusters on slice "<<slice<<" padrow "<<padrow<<endl;
	  for(j=0; j<ncl; j++)
	    {
	      AliTPCcluster *cluster = 0;
#ifdef asvversion
	      cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
#endif
	      if(cluster->GetZ() < 0 && slice < 18) continue;
	      if(cluster->GetZ() > 0 && slice > 17) continue;
	      if(cluster->IsUsed())
		continue;
	      
	      if(fWriteIdsToFile)
		idfile << cluster->GetLabel(0)<<' ';
	      
	      Float_t xyz[3] = {AliHLTTransform::Row2X(padrow),cluster->GetY(),cluster->GetZ()};
	      AliHLTTransform::Local2Raw(xyz,sector,row);

	      Int_t patch = AliHLTTransform::GetPatch(padrow);
	      Float_t padw = cluster->GetSigmaY2()/pow(AliHLTTransform::GetPadPitchWidth(patch),2);
	      Float_t timew = cluster->GetSigmaZ2()/pow(AliHLTTransform::GetZWidth(),2);
	      
	      Int_t buff;
	      //Write pad
	      buff = (Int_t)rint(xyz[1]*10);
	      if(buff<0)
		{
		  cerr<<"AliHLTOfflineDataCompressor:WriteRemaining : Wrong pad value "<<buff<<endl;
		  exit(5);
		}
	      OutputBits(output,buff,11);

	      //Write time
	      buff = (Int_t)rint(xyz[2]*10);
	      if(buff<0)
		{
		  cerr<<"AliHLTOfflineDataCompressor:WriteRemaining : Wrong time value "<<buff<<endl;
		  exit(5);
		}
	      OutputBits(output,buff,13);

	      //Write widths
	      buff = (Int_t)rint(padw*100);
	      OutputBits(output,buff,8);
	      buff = (Int_t)rint(timew*100);
	      OutputBits(output,buff,8);
	      
	      //Write charge 
	      buff = (Int_t)cluster->GetQ();
	      if(buff >= 1<<14)
		buff = (1<<14)-1;
	      OutputBits(output,buff,14);
	      
	      fNunusedClusters++;
	    }
	}
      CloseOutputBitFile(output);
    }
  if(fWriteIdsToFile)
    {
      idfile << endl;
      idfile.close();
    }
}

void AliHLTOfflineDataCompressor::SelectRemainingClusters()
{
  //select the remaining clusters 
  //which were not compressed  
  cout<<"Cleaning up clusters"<<endl;
  Int_t nrows = AliHLTTransform::GetNRows();
  Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
  
  Int_t sector,row,sec;

#ifdef asvversion
  AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
#endif
  for(Int_t slice=0; slice<36; slice++)
    {
      for(Int_t padrow=0; padrow < nrows; padrow++)
	{
	  	  
	  AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
	  sec=sector;
	  
	  if(fMarian == kFALSE)
	    {
	      if(slice >= 18)
		sec -= 18;
	      if(sec >= 18)
		sec -= 18;
	    }
	  Int_t ncl = 0;
#ifdef asvversion
	  ncl=tracker->GetNClusters(sec,row);
#endif
	  for(Int_t j=0; j<ncl; j++)
	    {
	      AliTPCcluster *cluster = 0;
#ifdef asvversion
	      cluster = (AliTPCcluster*)tracker->GetCluster(sec,row,j);
#endif
	      if(cluster->IsUsed())
		continue;

	      //Check the widths (errors) of the cluster, and remove big bastards:
	      Float_t xyw = cluster->GetSigmaY2() / pow(AliHLTTransform::GetPadPitchWidth(AliHLTTransform::GetPatch(padrow)),2);
	      Float_t zw  = cluster->GetSigmaZ2() / pow(AliHLTTransform::GetZWidth(),2);
	      if(xyw >= 2.55 || zw >= 2.55)//Because we use 1 byte to store
		{
		  cluster->Use();
		  continue;
		}
	      
	      //if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
	      
	      if(padrow == nrows - 1 || padrow == nrows - 1 - gap ||                 //First seeding
		 padrow == nrows - 1 - shift || padrow == nrows - 1 - gap - shift)   //Second seeding
		continue;
	      
	      if(cluster->GetZ() < 0 && slice < 18) continue;
	      if(cluster->GetZ() > 0 && slice > 17) continue;
	      if(cluster->IsUsed())
		continue;
	      
	      cluster->Use();
	    }
	}
      
    }
}
Commit	Line	Data
0a86fbb7	1	// @(#) $Id$
0a86fbb7	2
c3d96f5a	3	//_____________________________________________________________
c3d96f5a	4	//
4aa41877	5	// AliHLTOfflineDataCompression
c3d96f5a	6	//
	7	// Class to compress data with offline tracks
	8	// as seeds.
	9	//
0a86fbb7	10	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	11	//*-- Copyright &copy ALICE HLT Group
	12
4aa41877	13	#include "AliHLTStandardIncludes.h"
0a86fbb7	14
0a86fbb7	15	#include <AliTPCParamSR.h>
	16	#include <AliTPCClustersArray.h>
	17	#include <AliTPCcluster.h>
	18	#include <AliTPCClustersRow.h>
	19	#include <AliTPC.h>
	20	#include <AliTPCv2.h>
b2a02bce	21	#include <AliTPCcluster.h>
	22	#include <AliTPCtracker.h>
	23	#include <AliTPCclusterMI.h>
	24	#include <AliTPCtrackerMI.h>
	25	#include <AliKalmanTrack.h>
0a86fbb7	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>
1f1942b8	33	#include <TH1F.h>
0a86fbb7	34
4aa41877	35	#include "AliHLTTransform.h"
	36	#include "AliHLTModelTrack.h"
	37	#include "AliHLTCompress.h"
	38	#include "AliHLTTrackArray.h"
6f388e0d	39	#include "bitio.h"
4aa41877	40	#include "AliHLTOfflineDataCompressor.h"
0a86fbb7	41
0bd0c1ef	42	#if __GNUC__ == 3
0a86fbb7	43	using namespace std;
	44	#endif
	45
4aa41877	46	ClassImp(AliHLTOfflineDataCompressor)
0a86fbb7	47
4aa41877	48	AliHLTOfflineDataCompressor::AliHLTOfflineDataCompressor()
c3d96f5a	49	{ //constructor
0a86fbb7	50	fMarian = kFALSE;
	51	fTracker=0;
	52	}
	53
4aa41877	54	AliHLTOfflineDataCompressor::AliHLTOfflineDataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape,Bool_t MI)
4aa41877	55	: AliHLTDataCompressor(path,keep,writeshape)
c3d96f5a	56	{ //constructor
0a86fbb7	57	fMarian = MI;
	58	fTracker=0;
	59	}
	60
4aa41877	61	AliHLTOfflineDataCompressor::~AliHLTOfflineDataCompressor()
c3d96f5a	62	{ //deconstructor
0a86fbb7	63	if(fTracker)
	64	{
	65	fTracker->UnloadClusters();
	66	delete fTracker;
	67	}
	68	}
	69
	70
4aa41877	71	void AliHLTOfflineDataCompressor::LoadData(Int_t event,Bool_t sp)
0a86fbb7	72	{
	73	//Take offline reconstructed tracks as an input.
	74	//In this case, no remaining clusters are written.
	75
6f388e0d	76	if(fTracker)
	77	{
	78	fTracker->UnloadClusters();
	79	delete fTracker;
	80	}
0a86fbb7	81
6f388e0d	82	fSinglePatch = sp;
	83	fEvent = event;
	84
0a86fbb7	85	char filename[1024];
4aa41877	86	//AliKalmanTrack::SetConvConst(1000/0.299792458/AliHLTTransform::GetSolenoidField());
0a86fbb7	87	if(fMarian==kFALSE)
	88	sprintf(filename,"%s/offline/AliTPCclusters.root",fPath);
	89	else
	90	sprintf(filename,"%s/offline/AliTPCclustersMI.root",fPath);
	91
1f1942b8	92	bool compressed=0;
	93	if(compressed)
	94	{
4aa41877	95	cout<<"AliHLTOfflineDataCompressor::LoadData : Taking compressed offline files!!"<<endl;
1f1942b8	96	sprintf(filename,"%s/comp/offline/AliTPCclusters.root",fPath);
1f1942b8	97	}
0a86fbb7	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);
c84a5e9e	105	//fTracker->SetEventNumber(event);
b2a02bce	106	#ifdef asvversion
	107	fTracker->LoadClusters();
	108	#endif
0a86fbb7	109	if(fMarian==kTRUE)
0a86fbb7	110	{
b2a02bce	111	#ifdef asvversion
1f1942b8	112	AliTPCtrackerMI mitracker = (AliTPCtrackerMI)fTracker;
0a86fbb7	113	mitracker->LoadInnerSectors();
0a86fbb7	114	mitracker->LoadOuterSectors();
b2a02bce	115	#endif
0a86fbb7	116	}
0a86fbb7	117
c3d96f5a	118	const Int_t kMAX=20000;
c3d96f5a	119	Int_t nentr=0,i=0; TObjArray tarray(kMAX);
0a86fbb7	120	if(fMarian==kFALSE)
	121	sprintf(filename,"%s/offline/AliTPCtracks.root",fPath);
	122	else
	123	sprintf(filename,"%s/offline/AliTPCtracksMI.root",fPath);
1f1942b8	124
	125	if(compressed)
	126	sprintf(filename,"%s/comp/offline/AliTPCtracks.root",fPath);
	127
0a86fbb7	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
4aa41877	146	AliHLTTrackArray *comptracks = new AliHLTTrackArray("AliHLTModelTrack");
0a86fbb7	147	cout<<"Loaded "<<nentr<<" offline tracks"<<endl;
	148	Int_t slice,padrow;
	149	Int_t totcounter=0;
1f1942b8	150
0a86fbb7	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
4aa41877	171	AliHLTTransform::Sector2Slice(slice,padrow,sec,row);
4aa41877	172	Double_t par[5],xk=AliHLTTransform::Row2X(padrow);
0a86fbb7	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();
c3d96f5a	178	Float_t pt1=TMath::Abs(par[4]);
0a86fbb7	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;
1f1942b8	193	first[1] = par[0];
1f1942b8	194	first[2] = par[1];
0a86fbb7	195
4aa41877	196	AliHLTTransform::Local2Global(first,slice);
0a86fbb7	197
4aa41877	198	AliHLTModelTrack outtrack = (AliHLTModelTrack)comptracks->NextTrack();
0a86fbb7	199	outtrack->SetNHits(nhits);
0a86fbb7	200	outtrack->SetFirstPoint(first[0],first[1],first[2]);
c3d96f5a	201	outtrack->SetPt(1./pt1);
0a86fbb7	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;
4aa41877	257	if(!AliHLTTransform::Sector2Slice(slice,padrow,sec,row))
0a86fbb7	258	exit(5);
4aa41877	259	xyz[0] = AliHLTTransform::Row2X(padrow);
0a86fbb7	260
0a86fbb7	261	//cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" index "<<index<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
4aa41877	262	AliHLTTransform::Local2Raw(xyz,sec,row);
0a86fbb7	263	//cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
0a86fbb7	264
4aa41877	265	if(xyz[1] < -1 \|\| xyz[1] > AliHLTTransform::GetNPads(padrow) \|\|
4aa41877	266	xyz[2] < -1 \|\| xyz[2] > AliHLTTransform::GetNTimeBins())
0a86fbb7	267	{
4aa41877	268	cerr<<"AliHLTDataCompressor::FillOfflineData : Wrong time "<<xyz[2]<<" in slice "
0a86fbb7	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	}
1f1942b8	275
0a86fbb7	276	Float_t angle = 0;
4aa41877	277	AliHLTTransform::Local2GlobalAngle(&angle,slice);
4aa41877	278	if(!outtrack->CalculateReferencePoint(angle,AliHLTTransform::Row2X(padrow)))
0a86fbb7	279	{
4aa41877	280	cerr<<"AliHLTDataCompressor::FillOfflineData : Error in crossing point calc on slice "
0a86fbb7	281	<<slice<<" row "<<padrow<<endl;
	282	exit(5);
	283	}
1f1942b8	284
c3d96f5a	285	Float_t xyzcross[3] = {outtrack->GetPointX(),outtrack->GetPointY(),outtrack->GetPointZ()};
4aa41877	286	AliHLTTransform::Global2Raw(xyzcross,sec,row);
0a86fbb7	287	/*
c3d96f5a	288	if(fabs(xyzcross[1] - xyz[1]) > 10 \|\|
c3d96f5a	289	fabs(xyzcross[2] - xyz[2]) > 10)
0a86fbb7	290	{
4aa41877	291	cout<<"AliHLTDataCompressor::FillOfflineData : Wrong crossing slice "<<slice<<" padrow "
c3d96f5a	292	<<padrow<<" pad "<<xyz[1]<<" padhit "<<xyzcross[1]<<" time "<<xyz[2]<<" timehit "<<xyzcross[2]<<endl;
0a86fbb7	293	outtrack->Print();
	294	exit(5);
	295	}
	296	*/
c3d96f5a	297	//cout<<" crossing "<<xyzcross[0]<<" "<<xyzcross[1]<<" "<<xyzcross[2]<<endl;
	298	outtrack->SetPadHit(padrow,xyzcross[1]);
	299	outtrack->SetTimeHit(padrow,xyzcross[2]);
0a86fbb7	300
	301	if(fWriteClusterShape)
	302	{
	303	Float_t angle = outtrack->GetCrossingAngle(padrow,slice);
	304	outtrack->SetCrossingAngleLUT(padrow,angle);
	305	outtrack->CalculateClusterWidths(padrow,kTRUE);
4aa41877	306	Int_t patch = AliHLTTransform::GetPatch(padrow);
	307	Float_t sigmaY2 = cluster->GetSigmaY2() / pow(AliHLTTransform::GetPadPitchWidth(patch),2);
	308	Float_t sigmaZ2 = cluster->GetSigmaZ2() / pow(AliHLTTransform::GetZWidth(),2);
0a86fbb7	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;
b2a02bce	315	fNusedClusters++;
0a86fbb7	316	}
	317	}
	318
4aa41877	319	cout<<"AliHLTDataCompressor::FillOfflineData : Wrote "<<totcounter<<" clusters"<<endl;
0a86fbb7	320	//Write tracks to file
4aa41877	321	AliHLTCompress *comp = new AliHLTCompress(-1,-1,fPath,fWriteClusterShape,fEvent);
0a86fbb7	322	comp->WriteFile(comptracks);
	323	delete comp;
	324	delete comptracks;
0a86fbb7	325	}
0a86fbb7	326
4aa41877	327	void AliHLTOfflineDataCompressor::WriteRemaining(Bool_t select)
0a86fbb7	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	{
4aa41877	341	cerr<<"AliHLTOfflineDataCompressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
0a86fbb7	342	return;
0a86fbb7	343	}
b2a02bce	344
	345	ofstream idfile;
	346	if(fWriteIdsToFile)
	347	{
	348	sprintf(filename,"%s/comp/remains_ids.txt",fPath);
	349	idfile.open(filename);
	350	}
	351
0a86fbb7	352	cout<<"Writing remaining clusters "<<endl;
4aa41877	353	Int_t nrows = AliHLTTransform::GetNRows(),sector,row,sec;
1f1942b8	354	#ifdef asvversion
0a86fbb7	355	AliTPCtracker tracker = (AliTPCtracker)fTracker;
1f1942b8	356	#endif
0a86fbb7	357	for(Int_t slice=0; slice<=35; slice++)
	358	{
	359	sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
6f388e0d	360	BIT_FILE *output = OpenOutputBitFile(filename);
6f388e0d	361	if(!output)
0a86fbb7	362	{
4aa41877	363	cerr<<"AliHLTOfflineDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
0a86fbb7	364	exit(5);
	365	}
	366
6f388e0d	367	//Write number of padrows with clusters
	368	OutputBits(output,nrows,8);
	369
0a86fbb7	370	for(Int_t padrow=0; padrow < nrows; padrow++)
0a86fbb7	371	{
4aa41877	372	AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
0a86fbb7	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;
b2a02bce	383	Int_t ncl = 0;
b2a02bce	384	#ifdef asvversion
6f388e0d	385	ncl = tracker->GetNClusters(sec,row);
b2a02bce	386	#endif
0a86fbb7	387
	388	Int_t counter=0;
	389	Int_t j;
	390	for(j=0; j<ncl; j++)
	391	{
b2a02bce	392	AliTPCcluster *cluster = 0;
	393	#ifdef asvversion
	394	cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
	395	#endif
0a86fbb7	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	}
0a86fbb7	402
6f388e0d	403	OutputBits(output,padrow,8);//Write padrow #
	404	OutputBits(output,counter,10);//Write number of clusters on this padrow
	405
0a86fbb7	406	//cout<<"Found "<<counter<<" unused out of "<<ncl<<" clusters on slice "<<slice<<" padrow "<<padrow<<endl;
0a86fbb7	407	for(j=0; j<ncl; j++)
0a86fbb7	408	{
b2a02bce	409	AliTPCcluster *cluster = 0;
	410	#ifdef asvversion
	411	cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
	412	#endif
0a86fbb7	413	if(cluster->GetZ() < 0 && slice < 18) continue;
	414	if(cluster->GetZ() > 0 && slice > 17) continue;
	415	if(cluster->IsUsed())
1f1942b8	416	continue;
b2a02bce	417
	418	if(fWriteIdsToFile)
	419	idfile << cluster->GetLabel(0)<<' ';
	420
4aa41877	421	Float_t xyz[3] = {AliHLTTransform::Row2X(padrow),cluster->GetY(),cluster->GetZ()};
4aa41877	422	AliHLTTransform::Local2Raw(xyz,sector,row);
6f388e0d	423
4aa41877	424	Int_t patch = AliHLTTransform::GetPatch(padrow);
	425	Float_t padw = cluster->GetSigmaY2()/pow(AliHLTTransform::GetPadPitchWidth(patch),2);
	426	Float_t timew = cluster->GetSigmaZ2()/pow(AliHLTTransform::GetZWidth(),2);
6f388e0d	427
	428	Int_t buff;
	429	//Write pad
	430	buff = (Int_t)rint(xyz[1]*10);
	431	if(buff<0)
	432	{
4aa41877	433	cerr<<"AliHLTOfflineDataCompressor:WriteRemaining : Wrong pad value "<<buff<<endl;
6f388e0d	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	{
4aa41877	442	cerr<<"AliHLTOfflineDataCompressor:WriteRemaining : Wrong time value "<<buff<<endl;
6f388e0d	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
b2a02bce	459	fNunusedClusters++;
0a86fbb7	460	}
0a86fbb7	461	}
6f388e0d	462	CloseOutputBitFile(output);
0a86fbb7	463	}
b2a02bce	464	if(fWriteIdsToFile)
	465	{
	466	idfile << endl;
	467	idfile.close();
	468	}
0a86fbb7	469	}
0a86fbb7	470
4aa41877	471	void AliHLTOfflineDataCompressor::SelectRemainingClusters()
0a86fbb7	472	{
c3d96f5a	473	//select the remaining clusters
c3d96f5a	474	//which were not compressed
0a86fbb7	475	cout<<"Cleaning up clusters"<<endl;
4aa41877	476	Int_t nrows = AliHLTTransform::GetNRows();
0a86fbb7	477	Int_t gap=(Int_t)(0.125nrows), shift=(Int_t)(0.5gap);
	478
	479	Int_t sector,row,sec;
1f1942b8	480
1f1942b8	481	#ifdef asvversion
0a86fbb7	482	AliTPCtracker tracker = (AliTPCtracker)fTracker;
1f1942b8	483	#endif
0a86fbb7	484	for(Int_t slice=0; slice<36; slice++)
	485	{
	486	for(Int_t padrow=0; padrow < nrows; padrow++)
	487	{
6f388e0d	488
4aa41877	489	AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
0a86fbb7	490	sec=sector;
	491
	492	if(fMarian == kFALSE)
	493	{
	494	if(slice >= 18)
	495	sec -= 18;
	496	if(sec >= 18)
	497	sec -= 18;
	498	}
b2a02bce	499	Int_t ncl = 0;
	500	#ifdef asvversion
	501	ncl=tracker->GetNClusters(sec,row);
	502	#endif
0a86fbb7	503	for(Int_t j=0; j<ncl; j++)
0a86fbb7	504	{
6f388e0d	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:
4aa41877	513	Float_t xyw = cluster->GetSigmaY2() / pow(AliHLTTransform::GetPadPitchWidth(AliHLTTransform::GetPatch(padrow)),2);
4aa41877	514	Float_t zw = cluster->GetSigmaZ2() / pow(AliHLTTransform::GetZWidth(),2);
6f388e0d	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
0a86fbb7	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	}