[u/mrichter/AliRoot.git] / HLT / BASE / AliHLTSpacePointContainer.cxx

// $Id$

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
//*                  for The ALICE HLT Project.                            *
//*                                                                        *
//* Permission to use, copy, modify and distribute this software and its   *
//* documentation strictly for non-commercial purposes is hereby granted   *
//* without fee, provided that the above copyright notice appears in all   *
//* copies and that both the copyright notice and this permission notice   *
//* appear in the supporting documentation. The authors make no claims     *
//* about the suitability of this software for any purpose. It is          *
//* provided "as is" without express or implied warranty.                  *
//**************************************************************************

/// @file   AliHLTSpacePointContainer.cxx
/// @author Matthias Richter
/// @date   2011-04-29
/// @brief  Base helper class for handling of HLT space point data blocks
///

#include "AliHLTSpacePointContainer.h"
#include "AliHLTComponent.h"
#include "TFile.h"
#include "TString.h"
#include "TArrayC.h"
#include "TObjArray.h"
#include "TH2F.h"
#include "TMath.h"
#include "TMarker.h"
#include "TTree.h"
#include <memory>
#include <algorithm>
#include <iostream>
#include <iomanip>

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTSpacePointContainer)

AliHLTSpacePointContainer::AliHLTSpacePointContainer()
  : TObject(), AliHLTLogging()
  , fBuffers()
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
}

AliHLTSpacePointContainer::AliHLTSpacePointContainer(const AliHLTSpacePointContainer& c)
  : TObject(c), AliHLTLogging()
  , fBuffers()
{
  /// copy constructor

  // currently no need to copy or allocate, the fBuffers of the source object are expected
  // to be persistent throughout the lifetime of the new object. Only pointers are copied.
}

AliHLTSpacePointContainer& AliHLTSpacePointContainer::operator=(const AliHLTSpacePointContainer& c)
{
  /// assignment operator
  if (&c==this) return *this;

  // currently no need to copy or allocate
  return *this;
}

AliHLTSpacePointContainer::~AliHLTSpacePointContainer()
{
  // destructor
  Clear();
}

void AliHLTSpacePointContainer::Clear(Option_t * /*option*/)
{
  // internal cleanup
  vector<TArrayC*>::iterator element=fBuffers.begin();
  while (element!=fBuffers.end()) {
    TArrayC* pBuffer=*element;
    element=fBuffers.erase(element);
    delete pBuffer;
  }
}

void AliHLTSpacePointContainer::Print(Option_t *option) const
{
  // print info
  Print(cout, option);
}

void AliHLTSpacePointContainer::Print(ostream& out, Option_t */*option*/) const
{
  // print to stream
  out << "AliHLTSpacePointContainer::Print" << endl;
}

int AliHLTSpacePointContainer::AddInputBlock(const char* filename, AliHLTComponentDataType dt, unsigned specification)
{
  // open block from file and add to collection
  if (!filename) return -EINVAL;
  
  TString input=filename;
  input+="?filetype=raw";
  std::auto_ptr<TFile> pFile(new TFile(input));
  if (!pFile.get()) return -ENOMEM;
  if (pFile->IsZombie()) return -ENOENT;

  int iResult=0;
  pFile->Seek(0);
  std::auto_ptr<TArrayC> buffer(new TArrayC);
  if (!buffer.get()) return -ENOMEM;

  buffer->Set(pFile->GetSize());
  if (pFile->ReadBuffer(buffer->GetArray(), buffer->GetSize())==0) {
    AliHLTComponentBlockData bd;
    AliHLTComponent::FillBlockData(bd);
    bd.fPtr=buffer->GetArray();
    bd.fSize=buffer->GetSize();
    bd.fDataType=dt;
    bd.fSpecification=specification;
    HLTDebug("adding data block %d byte(s) from file %s", pFile->GetSize(), filename);
    iResult=AddInputBlock(&bd);
  } else {
    HLTError("failed reading %d byte(s) from file %s", pFile->GetSize(), filename);
    iResult=-ENODATA;
  }

  fBuffers.push_back(buffer.release());
  return iResult;
}

int AliHLTSpacePointContainer::AddInputBlocks(const char* listfile, AliHLTComponentDataType dt)
{
  // open blank separated list of files and add data
  ifstream list(listfile);
  if (!list.good()) return -ENOENT;

  int count=0;
  TString file;
  while (file.ReadLine(list)) {
    HLTInfo("adding data from file %s", file.Data());
    int iResult=AddInputBlock(file.Data(), dt, kAliHLTVoidDataSpec);
    if (iResult<0) {
      HLTInfo("failed to add data from file %s: error %d", file.Data(), iResult);
      return iResult;
    }
    count+=iResult;
  }

  // std::auto_ptr<TObjArray> tokens(files.Tokenize(" "));
  // if (!tokens.get()) return 0;
  // for (int i=0; i<tokens->GetEntriesFast(); i++) {
  //   if (!tokens->At(i)) continue;
  //   cout << "adding data from file " << tokens->At(i)->GetName() << endl;
  //   AddInputBlock(tokens->At(i)->GetName(), dt, kAliHLTVoidDataSpec);
  // }

  return count;
}

AliHLTUInt8_t* AliHLTSpacePointContainer::Alloc(int size)
{
  /// alloc memory for a space point data block
  TArrayC* buffer=new TArrayC;
  if (!buffer) return NULL;

  buffer->Set(size);
  fBuffers.push_back(buffer);
  return reinterpret_cast<AliHLTUInt8_t*>(buffer->GetArray());
}

AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByMask(AliHLTUInt32_t /*mask*/, bool /*bAlloc*/) const
{
  /// create a collection of clusters for a space point mask
  /// default implementation, nothing to do
  return NULL;
}

AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByTrack(int /*trackId*/, bool /*bAlloc*/) const
{
  /// create a collection of clusters for a specific track
  /// default implementation, nothing to do
  return NULL;
}

AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByMC(int /*mcId*/, bool /*bAlloc*/) const
{
  /// create a collection of clusters for a specific MC track
  /// default implementation, nothing to do
  return NULL;
}

AliHLTSpacePointContainer* AliHLTSpacePointContainer::UsedClusters(bool /*bAlloc*/) const
{
  /// create a collection of all used clusters
  /// default implementation, nothing to do
  return NULL;
}

AliHLTSpacePointContainer* AliHLTSpacePointContainer::UnusedClusters(bool /*bAlloc*/) const
{
  /// create a collection of all unused clusters
  /// default implementation, nothing to do
  return NULL;
}

int AliHLTSpacePointContainer::MarkUsed(const AliHLTUInt32_t* /*clusterIDs*/, int /*arraySize*/)
{
  /// default implementation, nothing to do
  return -ENOSYS;
}

int AliHLTSpacePointContainer::SetTrackID(int /*trackID*/, const AliHLTUInt32_t* /*clusterIDs*/, int /*arraySize*/)
{
  /// default implementation, nothing to do
  return -ENOSYS;
}

int AliHLTSpacePointContainer::SetMCID(int /*mcID*/, const AliHLTUInt32_t* /*clusterIDs*/, int /*arraySize*/)
{
  /// default implementation, nothing to do
  return -ENOSYS;
}

int AliHLTSpacePointContainer::GetNumberOfSpacePoints() const
{
  // get number of space points
  vector<AliHLTUInt32_t> clusterIDs;
  if (GetClusterIDs(clusterIDs)<0) return 0;
  return clusterIDs.size();
}

bool AliHLTSpacePointContainer::Check(AliHLTUInt32_t id) const
{
  // check if space point exists
  vector<AliHLTUInt32_t> clusterIDs;
  if (GetClusterIDs(clusterIDs)<0) return false;
  return find(clusterIDs.begin(), clusterIDs.end(), id)!=clusterIDs.end();
}

TH1* AliHLTSpacePointContainer::DrawProjection(const char* plane, const vector<AliHLTUInt32_t>& selection) const
{
  // draw the projection of space points in a specified plane
  // "xy", "xz", "yz"
  
  int mode=0;
  if (strcmp(plane, "xy")==0) mode=0;
  else if (strcmp(plane, "xz")==0) mode=1;
  else if (strcmp(plane, "yz")==0) mode=2;
  else {
    HLTError("invalid plane specification %s, allowed 'xy', 'xz', 'yz'", plane);
    return NULL;
  }
  
  float maxXAxis=100.0;
  float maxYAxis=100.0;
  vector<AliHLTUInt32_t> clusterIDs;
  GetClusterIDs(clusterIDs);
  vector<AliHLTUInt32_t>::iterator clusterID=clusterIDs.begin();
  while (clusterID!=clusterIDs.end()) {
    vector<AliHLTUInt32_t>::const_iterator element=selection.begin();
    for (; element!=selection.end(); element++) {
      if (((*element)&0xffff0000)==((*clusterID)&0xffff0000)) break;
    }
    if (element==selection.end() && selection.size()>0) {
      clusterID=clusterIDs.erase(clusterID);
      continue;
    }
    float XValue=0.0;
    float YValue=0.0;
    if (mode==0) {
      XValue=GetX(*clusterID);
      YValue=GetY(*clusterID);
    } else if (mode==1) {
      XValue=GetX(*clusterID);
      YValue=GetZ(*clusterID);
    } else if (mode==2) {
      XValue=GetY(*clusterID);
      YValue=GetZ(*clusterID);
    }
    if (maxXAxis<XValue) maxXAxis=XValue;
    if (maxYAxis<YValue) maxYAxis=YValue;
    clusterID++;
  }
  if (maxXAxis<maxYAxis) maxXAxis=maxYAxis;
  else maxYAxis=maxXAxis;

  TH2F* histogram=new TH2F(plane, plane, 1000, -maxXAxis, maxXAxis, 1000, -maxYAxis, maxYAxis);
  if (!histogram) return NULL;
    if (mode==0) {
      histogram->GetXaxis()->SetTitle("X [cm]");
      histogram->GetYaxis()->SetTitle("Y [cm]");
    } else if (mode==1) {
      histogram->GetXaxis()->SetTitle("Z [cm]");
      histogram->GetYaxis()->SetTitle("X [cm]");
    } else if (mode==2) {
      histogram->GetXaxis()->SetTitle("Z [cm]");
      histogram->GetYaxis()->SetTitle("Y [cm]");
    }

  for (clusterID=clusterIDs.begin(); clusterID!=clusterIDs.end(); clusterID++) {
    float XValue=GetX(*clusterID);
    float YValue=GetY(*clusterID);
    float ZValue=GetZ(*clusterID);
    Float_t phi     = GetPhi(*clusterID);
    Float_t cos     = TMath::Cos( phi );
    Float_t sin     = TMath::Sin( phi );
    if (mode==0) {
      histogram->SetTitle("XY projection");
      histogram->Fill(cos*XValue - sin*YValue, sin*XValue + cos*YValue);
    } else if (mode==1) {
      // should be maybe 'ZX' but somehow 'XZ' is more 'lingual convenient'
      histogram->SetTitle("XZ projection");
      histogram->Fill(ZValue, cos*XValue - sin*YValue);
    } else if (mode==2) {
      // same comment
      histogram->SetTitle("YZ projection");
      histogram->Fill(ZValue, sin*XValue + cos*YValue);
    }
  }

  return histogram;
}

void AliHLTSpacePointContainer::Draw(Option_t *option)
{
  /// Inherited from TObject, draw clusters
  float scale=250;
  float center[2]={0.5,0.5};
  int markersize=1;
  int markercolor=5;

  TString strOption(option);
  std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
  if (!tokens.get()) return;
  for (int i=0; i<tokens->GetEntriesFast(); i++) {
    if (!tokens->At(i)) continue;
    const char* key="";
    TString arg=tokens->At(i)->GetName();

    key="scale=";
    if (arg.BeginsWith(key)) {
      arg.ReplaceAll(key, "");
      scale=arg.Atof();
    }
    key="centerx=";
    if (arg.BeginsWith(key)) {
      arg.ReplaceAll(key, "");
      center[0]=arg.Atof();
    }
    key="centery=";
    if (arg.BeginsWith(key)) {
      arg.ReplaceAll(key, "");
      center[1]=arg.Atof();
    }
    key="markersize=";
    if (arg.BeginsWith(key)) {
      arg.ReplaceAll(key, "");
      markersize=arg.Atoi();
    }
    key="markercolor=";
    if (arg.BeginsWith(key)) {
      arg.ReplaceAll(key, "");
      markercolor=arg.Atoi();
    }
  }

  vector<AliHLTUInt32_t> clusterIDs;
  GetClusterIDs(clusterIDs);
  for (vector<AliHLTUInt32_t>::iterator clusterID=clusterIDs.begin();
       clusterID!=clusterIDs.end();
       clusterID++) {
    float clusterphi=GetPhi(*clusterID);
    float cosphi=TMath::Cos(clusterphi);
    float sinphi=TMath::Sin(clusterphi);
    float clusterx=GetX(*clusterID);
    float clustery=GetY(*clusterID);
    //cout << " " << *clusterID << " " << clusterx << " " << clustery << " " << clusterphi << endl;
    // rotate
    float pointx= clusterx*sinphi + clustery*cosphi;
    float pointy=-clusterx*cosphi + clustery*sinphi;

    TMarker* m=new TMarker(pointx/(2*scale)+center[0], pointy/(2*scale)+center[1], 3);
    m->SetMarkerSize(markersize);
    m->SetMarkerColor(markercolor);
    m->Draw("same");    
  }  
}

TTree* AliHLTSpacePointContainer::FillTree(const char* name, const char* title)
{
  // create tree and fill the space point coordinates
  TString treename=name;
  TString treetitle=title;
  if (treename.IsNull()) treename="spacepoints";
  if (treetitle.IsNull()) treetitle="HLT space point coordinates";
  std::auto_ptr<TTree> tree(new TTree(treename, treetitle));
  if (!tree.get()) return NULL;

  const unsigned dimension=8;
  float values[dimension];
  memset(values, 0, sizeof(values));
  const char* names[dimension]={
    "x", "y", "z", "sigmaY2", "sigmaZ2", "charge", "qmax", "alpha"
  };

  for (unsigned i=0; i<dimension; i++) {
    TString type=names[i]; type+="/F";
    tree->Branch(names[i], &values[i], type);
  }

  const vector<AliHLTUInt32_t>* clusterIDs=GetClusterIDs(0xffffffff);
  for (vector<AliHLTUInt32_t>::const_iterator clusterID=clusterIDs->begin();
       clusterID!=clusterIDs->end();
       clusterID++) {
    unsigned pos=0;
    values[pos++]=GetX(*clusterID);
    values[pos++]=GetY(*clusterID);
    values[pos++]=GetZ(*clusterID);
    values[pos++]=GetYWidth(*clusterID);
    values[pos++]=GetZWidth(*clusterID);
    values[pos++]=GetCharge(*clusterID);
    values[pos++]=GetMaxSignal(*clusterID);
    values[pos++]=GetPhi(*clusterID);

    tree->Fill();
  }

  return tree.release();
}

ClassImp(AliHLTSpacePointContainer::AliHLTSpacePointGrid)

AliHLTSpacePointContainer::AliHLTSpacePointGrid::AliHLTSpacePointGrid(float maxX, float stepX,
								      float maxY, float stepY,
								      float maxZ, float stepZ,
								      int initialDataSize)
  : fMaxX(maxX)
  , fStepX(stepX)
  , fMaxY(maxY)
  , fStepY(stepY)
  , fMaxZ(maxZ)
  , fStepZ(stepZ)
  , fDimX(0)
  , fDimY(0)
  , fDimZ(0)
  , fCells(NULL)
  , fCellDimension(0)
  , fData(NULL)
  , fDataDimension(initialDataSize)
  , fCount(0)
  , fIterator()
  , fIteratorEnd()
{
  // constructor
  if (fMaxX>0. && fMaxY>0. && fMaxZ>0 &&
      fStepX>0. && fStepY>0. && fStepZ>0) {
    fDimX=ceil(fMaxX/fStepX);
    fDimY=ceil(fMaxY/fStepY);
    fDimZ=ceil(fMaxZ/fStepZ);

    fCellDimension=fDimX*fDimY*fDimZ;
    fCells=new AliHLTSpacePointCell[fCellDimension];
    if (fDataDimension<0) fDataDimension=fgkDefaultDataSize;
    fData=new AliHLTSpacePointGrid::ValueType[fDataDimension];
    Clear();
  }
}

const int AliHLTSpacePointContainer::AliHLTSpacePointGrid::fgkDefaultDataSize=10000;

AliHLTSpacePointContainer::AliHLTSpacePointGrid::~AliHLTSpacePointGrid()
{
  // destructor
  if (fData) delete [] fData;
  if (fCells) delete [] fCells;
}

int AliHLTSpacePointContainer::AliHLTSpacePointGrid::CountSpacePoint(float x, float y, float z)
{
  // increment counter of the cell where the spacepoint is
  int cell=GetCellIndex(x, y, z);
  if (cell<0 || !fCells || cell>=fCellDimension) return -EFAULT;
  if (fCells[cell].fCount<0) fCells[cell].fCount=1;
  else fCells[cell].fCount++;
  return 0;
}

int AliHLTSpacePointContainer::AliHLTSpacePointGrid::IndexCells()
{
  // set the start index for data of every cell based on the counts
  if (!fCells || fCellDimension<=0) return -ENOBUFS;
  int offset=0;
  int cell=0;
  for (; cell<fCellDimension; cell++) {
    if (fCells[cell].fCount<0) continue;
    fCells[cell].fStartIndex=offset;
    offset+=fCells[cell].fCount;
    fCells[cell].fFilled=0;
  }

  if (offset>fDataDimension) {
    // grow the data array
    auto_ptr<AliHLTSpacePointGrid::ValueType> newArray(new AliHLTSpacePointGrid::ValueType[offset]);
    if (newArray.get()) {
      memcpy(newArray.get(), fData, fDataDimension);
      memset(newArray.get()+fDataDimension, 0, (offset-fDataDimension)*sizeof(AliHLTSpacePointGrid::ValueType));
      delete fData;
      fData=newArray.release();
      fDataDimension=offset;
    } else {
      for (cell=0; cell<fCellDimension; cell++) {
	fCells[cell].fStartIndex=-1;
      }
    }
  }
  return 0;
}

int AliHLTSpacePointContainer::AliHLTSpacePointGrid::AddSpacePoint(AliHLTSpacePointContainer::AliHLTSpacePointGrid::ValueType t,
								   float x, float y, float z)
{
  // add spacepoint, all spacepoints must have been counted before
  int cell=GetCellIndex(x, y, z);
  if (cell<0 || !fCells || cell>=fCellDimension) return -EFAULT;
  if (fCells[cell].fFilled==fCells[cell].fCount) return -ENOSPC;
  if (fCells[cell].fStartIndex<0 && IndexCells()<0) return -EACCES;
  int offset=fCells[cell].fStartIndex+fCells[cell].fFilled;
  fData[offset]=t;
  fCells[cell].fFilled++;
  fCount++;
  return 0;
}

void AliHLTSpacePointContainer::AliHLTSpacePointGrid::Clear(const char* /*option*/)
{
  // clear internal data
  if (fCells) memset(fCells, 0xff, fCellDimension*sizeof(AliHLTSpacePointCell));
  if (fData) memset(fData, 0, fDataDimension*sizeof(AliHLTSpacePointGrid::ValueType));
  fCount=0;
}

void AliHLTSpacePointContainer::AliHLTSpacePointGrid::Print(const char* /*option*/)
{
  // print info
  bool bPrintEmpty=false;
  cout << "AliHLTSpacePointGrid: " << (fCells?fCellDimension:0) << " cells" << endl;
  cout << "   x: " << fDimX << " [0," << fMaxX << "]" << endl;
  cout << "   y: " << fDimY << " [0," << fMaxY << "]" << endl;
  cout << "   z: " << fDimZ << " [0," << fMaxZ << "]" << endl;
  cout << "   " << GetNumberOfSpacePoints(0, fCellDimension) << " point(s)" << endl;
  if (fCells) {
    for (int i=0; i<fCellDimension; i++) {
      if (!bPrintEmpty && fCells[i].fCount<=0) continue;
      cout << "     " << setfill(' ') << setw(7) << setprecision(0) << i << " (" 
	   << " " << setw(3) << GetLowerBoundX(i)
	   << " " << setw(3) << GetLowerBoundY(i)
	   << " " << setw(4) << GetLowerBoundZ(i)
	   << "): ";
      cout << setw(3) << fCells[i].fCount << " entries, " << setw(3) << fCells[i].fFilled << " filled";
      cout << "  start index " << setw(5) << fCells[i].fStartIndex;
      cout << endl;
      if (fCells[i].fCount>0) {
	cout << "          ";
	for (iterator id=begin(GetLowerBoundX(i), GetLowerBoundY(i), GetLowerBoundZ(i));
	     id!=end(); id++) {
	  cout << " 0x" << hex << setw(8) << setfill('0') << id.Data();
	}
	cout  << dec << endl;
      }
    }
  }
}

ostream& operator<<(ostream &out, const AliHLTSpacePointContainer& c)
{
  c.Print(out);
  return out;
}
Commit	Line	Data
488caaa0	1	// $Id$
	2
	3	//**************************************************************************
	4	//* This file is property of and copyright by the ALICE HLT Project *
	5	//* ALICE Experiment at CERN, All rights reserved. *
	6	//* *
	7	//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
	8	//* for The ALICE HLT Project. *
	9	//* *
	10	//* Permission to use, copy, modify and distribute this software and its *
	11	//* documentation strictly for non-commercial purposes is hereby granted *
	12	//* without fee, provided that the above copyright notice appears in all *
	13	//* copies and that both the copyright notice and this permission notice *
	14	//* appear in the supporting documentation. The authors make no claims *
	15	//* about the suitability of this software for any purpose. It is *
	16	//* provided "as is" without express or implied warranty. *
	17	//**************************************************************************
	18
	19	/// @file AliHLTSpacePointContainer.cxx
	20	/// @author Matthias Richter
	21	/// @date 2011-04-29
	22	/// @brief Base helper class for handling of HLT space point data blocks
	23	///
	24
	25	#include "AliHLTSpacePointContainer.h"
	26	#include "AliHLTComponent.h"
	27	#include "TFile.h"
	28	#include "TString.h"
	29	#include "TArrayC.h"
	30	#include "TObjArray.h"
27d45f0c	31	#include "TH2F.h"
27d45f0c	32	#include "TMath.h"
5f195a83	33	#include "TMarker.h"
6fa220b4	34	#include "TTree.h"
488caaa0	35	#include <memory>
e1e03704	36	#include <algorithm>
488caaa0	37	#include <iostream>
86e69fe5	38	#include <iomanip>
488caaa0	39
	40	/** ROOT macro for the implementation of ROOT specific class methods */
	41	ClassImp(AliHLTSpacePointContainer)
	42
	43	AliHLTSpacePointContainer::AliHLTSpacePointContainer()
	44	: TObject(), AliHLTLogging()
	45	, fBuffers()
	46	{
	47	// see header file for class documentation
	48	// or
	49	// refer to README to build package
	50	// or
	51	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	52	}
	53
	54	AliHLTSpacePointContainer::AliHLTSpacePointContainer(const AliHLTSpacePointContainer& c)
	55	: TObject(c), AliHLTLogging()
	56	, fBuffers()
	57	{
	58	/// copy constructor
	59
	60	// currently no need to copy or allocate, the fBuffers of the source object are expected
	61	// to be persistent throughout the lifetime of the new object. Only pointers are copied.
	62	}
	63
	64	AliHLTSpacePointContainer& AliHLTSpacePointContainer::operator=(const AliHLTSpacePointContainer& c)
	65	{
	66	/// assignment operator
	67	if (&c==this) return *this;
	68
	69	// currently no need to copy or allocate
	70	return *this;
	71	}
	72
	73	AliHLTSpacePointContainer::~AliHLTSpacePointContainer()
	74	{
	75	// destructor
	76	Clear();
	77	}
	78
	79	void AliHLTSpacePointContainer::Clear(Option_t * /option/)
	80	{
	81	// internal cleanup
	82	vector<TArrayC*>::iterator element=fBuffers.begin();
	83	while (element!=fBuffers.end()) {
	84	TArrayC* pBuffer=*element;
	85	element=fBuffers.erase(element);
	86	delete pBuffer;
	87	}
	88	}
	89
	90	void AliHLTSpacePointContainer::Print(Option_t *option) const
	91	{
	92	// print info
	93	Print(cout, option);
	94	}
	95
	96	void AliHLTSpacePointContainer::Print(ostream& out, Option_t /option*/) const
	97	{
	98	// print to stream
	99	out << "AliHLTSpacePointContainer::Print" << endl;
	100	}
	101
	102	int AliHLTSpacePointContainer::AddInputBlock(const char* filename, AliHLTComponentDataType dt, unsigned specification)
103	{
104	// open block from file and add to collection
105	if (!filename) return -EINVAL;
106
107	TString input=filename;
108	input+="?filetype=raw";
109	std::auto_ptr<TFile> pFile(new TFile(input));
110	if (!pFile.get()) return -ENOMEM;
111	if (pFile->IsZombie()) return -ENOENT;
112
113	int iResult=0;
114	pFile->Seek(0);
115	std::auto_ptr<TArrayC> buffer(new TArrayC);
116	if (!buffer.get()) return -ENOMEM;
117
118	buffer->Set(pFile->GetSize());
119	if (pFile->ReadBuffer(buffer->GetArray(), buffer->GetSize())==0) {
120	AliHLTComponentBlockData bd;
121	AliHLTComponent::FillBlockData(bd);
122	bd.fPtr=buffer->GetArray();
123	bd.fSize=buffer->GetSize();
124	bd.fDataType=dt;
125	bd.fSpecification=specification;
126	HLTDebug("adding data block %d byte(s) from file %s", pFile->GetSize(), filename);
127	iResult=AddInputBlock(&bd);
128	} else {
129	HLTError("failed reading %d byte(s) from file %s", pFile->GetSize(), filename);
130	iResult=-ENODATA;
131	}
132
133	fBuffers.push_back(buffer.release());
134	return iResult;
135	}
136
137	int AliHLTSpacePointContainer::AddInputBlocks(const char* listfile, AliHLTComponentDataType dt)
138	{
139	// open blank separated list of files and add data
140	ifstream list(listfile);
141	if (!list.good()) return -ENOENT;
142
143	int count=0;
144	TString file;
145	while (file.ReadLine(list)) {
146	HLTInfo("adding data from file %s", file.Data());
147	int iResult=AddInputBlock(file.Data(), dt, kAliHLTVoidDataSpec);
148	if (iResult<0) {
149	HLTInfo("failed to add data from file %s: error %d", file.Data(), iResult);
150	return iResult;
151	}
152	count+=iResult;
153	}
154
155	// std::auto_ptr<TObjArray> tokens(files.Tokenize(" "));
156	// if (!tokens.get()) return 0;
157	// for (int i=0; i<tokens->GetEntriesFast(); i++) {
158	// if (!tokens->At(i)) continue;
159	// cout << "adding data from file " << tokens->At(i)->GetName() << endl;
160	// AddInputBlock(tokens->At(i)->GetName(), dt, kAliHLTVoidDataSpec);
161	// }
162
163	return count;
164	}
165
3488f70b	166	AliHLTUInt8_t* AliHLTSpacePointContainer::Alloc(int size)
	167	{
	168	/// alloc memory for a space point data block
	169	TArrayC* buffer=new TArrayC;
	170	if (!buffer) return NULL;
	171
	172	buffer->Set(size);
	173	fBuffers.push_back(buffer);
	174	return reinterpret_cast<AliHLTUInt8_t*>(buffer->GetArray());
	175	}
	176
54ff4c01	177	AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByMask(AliHLTUInt32_t /mask/, bool /bAlloc/) const
	178	{
	179	/// create a collection of clusters for a space point mask
	180	/// default implementation, nothing to do
	181	return NULL;
	182	}
	183
488caaa0	184	AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByTrack(int /trackId/, bool /bAlloc/) const
	185	{
	186	/// create a collection of clusters for a specific track
	187	/// default implementation, nothing to do
	188	return NULL;
	189	}
	190
	191	AliHLTSpacePointContainer* AliHLTSpacePointContainer::SelectByMC(int /mcId/, bool /bAlloc/) const
	192	{
	193	/// create a collection of clusters for a specific MC track
	194	/// default implementation, nothing to do
	195	return NULL;
	196	}
	197
	198	AliHLTSpacePointContainer* AliHLTSpacePointContainer::UsedClusters(bool /bAlloc/) const
	199	{
	200	/// create a collection of all used clusters
	201	/// default implementation, nothing to do
	202	return NULL;
	203	}
	204
	205	AliHLTSpacePointContainer* AliHLTSpacePointContainer::UnusedClusters(bool /bAlloc/) const
	206	{
	207	/// create a collection of all unused clusters
	208	/// default implementation, nothing to do
	209	return NULL;
	210	}
	211
	212	int AliHLTSpacePointContainer::MarkUsed(const AliHLTUInt32_t* /clusterIDs/, int /arraySize/)
	213	{
	214	/// default implementation, nothing to do
	215	return -ENOSYS;
	216	}
	217
	218	int AliHLTSpacePointContainer::SetTrackID(int /trackID/, const AliHLTUInt32_t* /clusterIDs/, int /arraySize/)
	219	{
	220	/// default implementation, nothing to do
	221	return -ENOSYS;
	222	}
	223
	224	int AliHLTSpacePointContainer::SetMCID(int /mcID/, const AliHLTUInt32_t* /clusterIDs/, int /arraySize/)
	225	{
	226	/// default implementation, nothing to do
	227	return -ENOSYS;
	228	}
	229
5f195a83	230	int AliHLTSpacePointContainer::GetNumberOfSpacePoints() const
	231	{
	232	// get number of space points
	233	vector<AliHLTUInt32_t> clusterIDs;
	234	if (GetClusterIDs(clusterIDs)<0) return 0;
	235	return clusterIDs.size();
	236	}
	237
e1e03704	238	bool AliHLTSpacePointContainer::Check(AliHLTUInt32_t id) const
	239	{
	240	// check if space point exists
	241	vector<AliHLTUInt32_t> clusterIDs;
	242	if (GetClusterIDs(clusterIDs)<0) return false;
	243	return find(clusterIDs.begin(), clusterIDs.end(), id)!=clusterIDs.end();
	244	}
	245
27d45f0c	246	TH1* AliHLTSpacePointContainer::DrawProjection(const char* plane, const vector<AliHLTUInt32_t>& selection) const
	247	{
	248	// draw the projection of space points in a specified plane
	249	// "xy", "xz", "yz"
	250
	251	int mode=0;
	252	if (strcmp(plane, "xy")==0) mode=0;
	253	else if (strcmp(plane, "xz")==0) mode=1;
	254	else if (strcmp(plane, "yz")==0) mode=2;
	255	else {
	256	HLTError("invalid plane specification %s, allowed 'xy', 'xz', 'yz'", plane);
	257	return NULL;
	258	}
	259
	260	float maxXAxis=100.0;
	261	float maxYAxis=100.0;
	262	vector<AliHLTUInt32_t> clusterIDs;
	263	GetClusterIDs(clusterIDs);
	264	vector<AliHLTUInt32_t>::iterator clusterID=clusterIDs.begin();
	265	while (clusterID!=clusterIDs.end()) {
	266	vector<AliHLTUInt32_t>::const_iterator element=selection.begin();
	267	for (; element!=selection.end(); element++) {
	268	if (((element)&0xffff0000)==((clusterID)&0xffff0000)) break;
	269	}
	270	if (element==selection.end() && selection.size()>0) {
	271	clusterID=clusterIDs.erase(clusterID);
	272	continue;
	273	}
	274	float XValue=0.0;
	275	float YValue=0.0;
	276	if (mode==0) {
	277	XValue=GetX(*clusterID);
	278	YValue=GetY(*clusterID);
	279	} else if (mode==1) {
	280	XValue=GetX(*clusterID);
	281	YValue=GetZ(*clusterID);
	282	} else if (mode==2) {
	283	XValue=GetY(*clusterID);
	284	YValue=GetZ(*clusterID);
	285	}
	286	if (maxXAxis<XValue) maxXAxis=XValue;
	287	if (maxYAxis<YValue) maxYAxis=YValue;
	288	clusterID++;
	289	}
	290	if (maxXAxis<maxYAxis) maxXAxis=maxYAxis;
	291	else maxYAxis=maxXAxis;
	292
	293	TH2F* histogram=new TH2F(plane, plane, 1000, -maxXAxis, maxXAxis, 1000, -maxYAxis, maxYAxis);
	294	if (!histogram) return NULL;
	295	if (mode==0) {
	296	histogram->GetXaxis()->SetTitle("X [cm]");
	297	histogram->GetYaxis()->SetTitle("Y [cm]");
	298	} else if (mode==1) {
	299	histogram->GetXaxis()->SetTitle("Z [cm]");
	300	histogram->GetYaxis()->SetTitle("X [cm]");
	301	} else if (mode==2) {
	302	histogram->GetXaxis()->SetTitle("Z [cm]");
	303	histogram->GetYaxis()->SetTitle("Y [cm]");
	304	}
	305
	306	for (clusterID=clusterIDs.begin(); clusterID!=clusterIDs.end(); clusterID++) {
	307	float XValue=GetX(*clusterID);
	308	float YValue=GetY(*clusterID);
	309	float ZValue=GetZ(*clusterID);
310	Float_t phi = GetPhi(*clusterID);
311	Float_t cos = TMath::Cos( phi );
312	Float_t sin = TMath::Sin( phi );
313	if (mode==0) {
314	histogram->SetTitle("XY projection");
315	histogram->Fill(cosXValue - sinYValue, sinXValue + cosYValue);
316	} else if (mode==1) {
317	// should be maybe 'ZX' but somehow 'XZ' is more 'lingual convenient'
318	histogram->SetTitle("XZ projection");
319	histogram->Fill(ZValue, cosXValue - sinYValue);
320	} else if (mode==2) {
321	// same comment
322	histogram->SetTitle("YZ projection");
323	histogram->Fill(ZValue, sinXValue + cosYValue);
324	}
325	}
326
327	return histogram;
328	}
329
5f195a83	330	void AliHLTSpacePointContainer::Draw(Option_t *option)
	331	{
	332	/// Inherited from TObject, draw clusters
	333	float scale=250;
	334	float center[2]={0.5,0.5};
54ff4c01	335	int markersize=1;
54ff4c01	336	int markercolor=5;
5f195a83	337
	338	TString strOption(option);
	339	std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
	340	if (!tokens.get()) return;
	341	for (int i=0; i<tokens->GetEntriesFast(); i++) {
	342	if (!tokens->At(i)) continue;
	343	const char* key="";
	344	TString arg=tokens->At(i)->GetName();
	345
	346	key="scale=";
	347	if (arg.BeginsWith(key)) {
	348	arg.ReplaceAll(key, "");
	349	scale=arg.Atof();
	350	}
	351	key="centerx=";
	352	if (arg.BeginsWith(key)) {
	353	arg.ReplaceAll(key, "");
	354	center[0]=arg.Atof();
	355	}
	356	key="centery=";
	357	if (arg.BeginsWith(key)) {
	358	arg.ReplaceAll(key, "");
	359	center[1]=arg.Atof();
	360	}
54ff4c01	361	key="markersize=";
	362	if (arg.BeginsWith(key)) {
	363	arg.ReplaceAll(key, "");
	364	markersize=arg.Atoi();
	365	}
	366	key="markercolor=";
	367	if (arg.BeginsWith(key)) {
	368	arg.ReplaceAll(key, "");
	369	markercolor=arg.Atoi();
	370	}
5f195a83	371	}
	372
	373	vector<AliHLTUInt32_t> clusterIDs;
	374	GetClusterIDs(clusterIDs);
	375	for (vector<AliHLTUInt32_t>::iterator clusterID=clusterIDs.begin();
	376	clusterID!=clusterIDs.end();
	377	clusterID++) {
	378	float clusterphi=GetPhi(*clusterID);
	379	float cosphi=TMath::Cos(clusterphi);
	380	float sinphi=TMath::Sin(clusterphi);
	381	float clusterx=GetX(*clusterID);
	382	float clustery=GetY(*clusterID);
54ff4c01	383	//cout << " " << *clusterID << " " << clusterx << " " << clustery << " " << clusterphi << endl;
5f195a83	384	// rotate
	385	float pointx= clusterxsinphi + clusterycosphi;
	386	float pointy=-clusterxcosphi + clusterysinphi;
	387
	388	TMarker* m=new TMarker(pointx/(2scale)+center[0], pointy/(2scale)+center[1], 3);
54ff4c01	389	m->SetMarkerSize(markersize);
54ff4c01	390	m->SetMarkerColor(markercolor);
5f195a83	391	m->Draw("same");
	392	}
	393	}
	394
6fa220b4	395	TTree* AliHLTSpacePointContainer::FillTree(const char* name, const char* title)
	396	{
	397	// create tree and fill the space point coordinates
	398	TString treename=name;
	399	TString treetitle=title;
	400	if (treename.IsNull()) treename="spacepoints";
	401	if (treetitle.IsNull()) treetitle="HLT space point coordinates";
	402	std::auto_ptr<TTree> tree(new TTree(treename, treetitle));
	403	if (!tree.get()) return NULL;
	404
37c3ba28	405	const unsigned dimension=8;
6fa220b4	406	float values[dimension];
	407	memset(values, 0, sizeof(values));
	408	const char* names[dimension]={
37c3ba28	409	"x", "y", "z", "sigmaY2", "sigmaZ2", "charge", "qmax", "alpha"
6fa220b4	410	};
	411
	412	for (unsigned i=0; i<dimension; i++) {
	413	TString type=names[i]; type+="/F";
	414	tree->Branch(names[i], &values[i], type);
	415	}
	416
	417	const vector<AliHLTUInt32_t>* clusterIDs=GetClusterIDs(0xffffffff);
	418	for (vector<AliHLTUInt32_t>::const_iterator clusterID=clusterIDs->begin();
	419	clusterID!=clusterIDs->end();
	420	clusterID++) {
	421	unsigned pos=0;
	422	values[pos++]=GetX(*clusterID);
	423	values[pos++]=GetY(*clusterID);
	424	values[pos++]=GetZ(*clusterID);
37c3ba28	425	values[pos++]=GetYWidth(*clusterID);
37c3ba28	426	values[pos++]=GetZWidth(*clusterID);
6fa220b4	427	values[pos++]=GetCharge(*clusterID);
37c3ba28	428	values[pos++]=GetMaxSignal(*clusterID);
6fa220b4	429	values[pos++]=GetPhi(*clusterID);
	430
	431	tree->Fill();
	432	}
	433
	434	return tree.release();
	435	}
	436
86e69fe5	437	ClassImp(AliHLTSpacePointContainer::AliHLTSpacePointGrid)
	438
	439	AliHLTSpacePointContainer::AliHLTSpacePointGrid::AliHLTSpacePointGrid(float maxX, float stepX,
	440	float maxY, float stepY,
	441	float maxZ, float stepZ,
	442	int initialDataSize)
	443	: fMaxX(maxX)
	444	, fStepX(stepX)
	445	, fMaxY(maxY)
	446	, fStepY(stepY)
	447	, fMaxZ(maxZ)
	448	, fStepZ(stepZ)
	449	, fDimX(0)
	450	, fDimY(0)
	451	, fDimZ(0)
	452	, fCells(NULL)
	453	, fCellDimension(0)
	454	, fData(NULL)
	455	, fDataDimension(initialDataSize)
	456	, fCount(0)
	457	, fIterator()
	458	, fIteratorEnd()
	459	{
	460	// constructor
	461	if (fMaxX>0. && fMaxY>0. && fMaxZ>0 &&
	462	fStepX>0. && fStepY>0. && fStepZ>0) {
	463	fDimX=ceil(fMaxX/fStepX);
	464	fDimY=ceil(fMaxY/fStepY);
	465	fDimZ=ceil(fMaxZ/fStepZ);
	466
	467	fCellDimension=fDimXfDimYfDimZ;
	468	fCells=new AliHLTSpacePointCell[fCellDimension];
fc7ea273	469	if (fDataDimension<0) fDataDimension=fgkDefaultDataSize;
86e69fe5	470	fData=new AliHLTSpacePointGrid::ValueType[fDataDimension];
	471	Clear();
	472	}
	473	}
	474
fc7ea273	475	const int AliHLTSpacePointContainer::AliHLTSpacePointGrid::fgkDefaultDataSize=10000;
fc7ea273	476
86e69fe5	477	AliHLTSpacePointContainer::AliHLTSpacePointGrid::~AliHLTSpacePointGrid()
	478	{
	479	// destructor
	480	if (fData) delete [] fData;
	481	if (fCells) delete [] fCells;
	482	}
	483
	484	int AliHLTSpacePointContainer::AliHLTSpacePointGrid::CountSpacePoint(float x, float y, float z)
	485	{
	486	// increment counter of the cell where the spacepoint is
	487	int cell=GetCellIndex(x, y, z);
	488	if (cell<0 \|\| !fCells \|\| cell>=fCellDimension) return -EFAULT;
	489	if (fCells[cell].fCount<0) fCells[cell].fCount=1;
	490	else fCells[cell].fCount++;
	491	return 0;
	492	}
	493
	494	int AliHLTSpacePointContainer::AliHLTSpacePointGrid::IndexCells()
	495	{
	496	// set the start index for data of every cell based on the counts
	497	if (!fCells \|\| fCellDimension<=0) return -ENOBUFS;
	498	int offset=0;
	499	int cell=0;
	500	for (; cell<fCellDimension; cell++) {
	501	if (fCells[cell].fCount<0) continue;
	502	fCells[cell].fStartIndex=offset;
	503	offset+=fCells[cell].fCount;
	504	fCells[cell].fFilled=0;
	505	}
	506
	507	if (offset>fDataDimension) {
	508	// grow the data array
	509	auto_ptr<AliHLTSpacePointGrid::ValueType> newArray(new AliHLTSpacePointGrid::ValueType[offset]);
	510	if (newArray.get()) {
	511	memcpy(newArray.get(), fData, fDataDimension);
fc7ea273	512	memset(newArray.get()+fDataDimension, 0, (offset-fDataDimension)*sizeof(AliHLTSpacePointGrid::ValueType));
86e69fe5	513	delete fData;
	514	fData=newArray.release();
	515	fDataDimension=offset;
	516	} else {
	517	for (cell=0; cell<fCellDimension; cell++) {
	518	fCells[cell].fStartIndex=-1;
	519	}
	520	}
	521	}
	522	return 0;
	523	}
	524
	525	int AliHLTSpacePointContainer::AliHLTSpacePointGrid::AddSpacePoint(AliHLTSpacePointContainer::AliHLTSpacePointGrid::ValueType t,
	526	float x, float y, float z)
	527	{
	528	// add spacepoint, all spacepoints must have been counted before
	529	int cell=GetCellIndex(x, y, z);
	530	if (cell<0 \|\| !fCells \|\| cell>=fCellDimension) return -EFAULT;
	531	if (fCells[cell].fFilled==fCells[cell].fCount) return -ENOSPC;
	532	if (fCells[cell].fStartIndex<0 && IndexCells()<0) return -EACCES;
	533	int offset=fCells[cell].fStartIndex+fCells[cell].fFilled;
	534	fData[offset]=t;
	535	fCells[cell].fFilled++;
	536	fCount++;
	537	return 0;
	538	}
	539
	540	void AliHLTSpacePointContainer::AliHLTSpacePointGrid::Clear(const char* /option/)
	541	{
	542	// clear internal data
	543	if (fCells) memset(fCells, 0xff, fCellDimension*sizeof(AliHLTSpacePointCell));
	544	if (fData) memset(fData, 0, fDataDimension*sizeof(AliHLTSpacePointGrid::ValueType));
	545	fCount=0;
	546	}
	547
	548	void AliHLTSpacePointContainer::AliHLTSpacePointGrid::Print(const char* /option/)
	549	{
	550	// print info
	551	bool bPrintEmpty=false;
	552	cout << "AliHLTSpacePointGrid: " << (fCells?fCellDimension:0) << " cells" << endl;
	553	cout << " x: " << fDimX << " [0," << fMaxX << "]" << endl;
	554	cout << " y: " << fDimY << " [0," << fMaxY << "]" << endl;
	555	cout << " z: " << fDimZ << " [0," << fMaxZ << "]" << endl;
	556	cout << " " << GetNumberOfSpacePoints(0, fCellDimension) << " point(s)" << endl;
	557	if (fCells) {
	558	for (int i=0; i<fCellDimension; i++) {
	559	if (!bPrintEmpty && fCells[i].fCount<=0) continue;
	560	cout << " " << setfill(' ') << setw(7) << setprecision(0) << i << " ("
	561	<< " " << setw(3) << GetLowerBoundX(i)
	562	<< " " << setw(3) << GetLowerBoundY(i)
	563	<< " " << setw(4) << GetLowerBoundZ(i)
	564	<< "): ";
	565	cout << setw(3) << fCells[i].fCount << " entries, " << setw(3) << fCells[i].fFilled << " filled";
	566	cout << " start index " << setw(5) << fCells[i].fStartIndex;
	567	cout << endl;
	568	if (fCells[i].fCount>0) {
	569	cout << " ";
	570	for (iterator id=begin(GetLowerBoundX(i), GetLowerBoundY(i), GetLowerBoundZ(i));
	571	id!=end(); id++) {
	572	cout << " 0x" << hex << setw(8) << setfill('0') << id.Data();
	573	}
	574	cout << dec << endl;
	575	}
	576	}
577	}
578	}
579
488caaa0	580	ostream& operator<<(ostream &out, const AliHLTSpacePointContainer& c)
	581	{
	582	c.Print(out);
	583	return out;
	584	}