[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClusterTransformation.cxx

// $Id: AliHLTTPCClusterTransformation.cxx 41244 2010-05-14 08:13:35Z kkanaki $

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        *
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Kalliopi Kanaki <Kalliopi.Kanaki@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   AliHLTTPCClusterTransformation.cxx
    @author Kalliopi Kanaki, Sergey Gorbubnov
    @date   
    @brief 
*/


#include "AliHLTTPCClusterTransformation.h"
#include "AliHLTTPCTransform.h"

#include "AliTPCcalibDB.h"
#include "AliTPCTransform.h"
#include "AliTPCParam.h"
#include "AliTPCRecoParam.h"
#include "AliGeomManager.h"
#include <iostream>
#include <iomanip>

using namespace std;

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

AliHLTTPCClusterTransformation::AliHLTTPCClusterTransformation()
:
  fOfflineTransform(NULL),
  fOfflineTPCParam( NULL ),
  fLastSector(-1)
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt  
  fAliT[0] = 0.;
  fAliT[1] = 0.;
  fAliT[2] = 0.;
  SetRotationMatrix();
}
AliHLTTPCClusterTransformation::~AliHLTTPCClusterTransformation() 
{ 
  // see header file for class documentation
  //delete fOfflineTransform;
}


int  AliHLTTPCClusterTransformation::Init( double FieldBz, UInt_t TimeStamp )
{
  // Initialisation

  //delete fOfflineTransform;
  fOfflineTransform = 0;
  fOfflineTPCParam = 0;

  AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();

  if(!pCalib ) return -1;

  pCalib->SetExBField(FieldBz);
  
  if(!AliGeomManager::GetGeometry()){
     AliGeomManager::LoadGeometry();
  }

  if( !pCalib->GetTransform() ) return -2; 

  //fOfflineTransform = new AliTPCTransform (*pCalib->GetTransform());
  fOfflineTransform = pCalib->GetTransform();
  fOfflineTransform->SetCurrentRecoParam( AliTPCRecoParam::GetHLTParam() );
  fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
  fOfflineTPCParam = pCalib->GetParameters(); 
  if( !fOfflineTPCParam ) return -3;

  fOfflineTPCParam->Update();
  fOfflineTPCParam->ReadGeoMatrices();  

  fLastSector = -1;

  fAliT[0] = 0.;
  fAliT[1] = 0.;
  fAliT[2] = 0.;
  SetRotationMatrix();

  return 0;
}


void AliHLTTPCClusterTransformation::SetCurrentTimeStamp( UInt_t TimeStamp )
{
  // Set the current time stamp
  if( fOfflineTransform ) fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
}


int  AliHLTTPCClusterTransformation::Transform( int Slice, int Row, float Pad, float Time, float XYZ[] )
{
  // Convert row, pad, time to X Y Z
   	   
  Int_t sector=-99, thisrow=-99;
  AliHLTTPCTransform::Slice2Sector( Slice, Row, sector, thisrow);

  if( !fOfflineTransform ){   	   	   
    AliHLTTPCTransform::Raw2Local( XYZ, sector, thisrow, Pad, Time); 
    if(Slice>17) XYZ[1]= - XYZ[1];	   
    return 0;
  }

  Int_t iSector[1]= {sector};   
  Double_t x[3] = { thisrow, Pad, Time }; 
  fOfflineTransform->Transform(x,iSector,0,1);
	  
  if( sector!= fLastSector ){
    if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
      TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
      if ( alignment ){
	const Double_t *tr = alignment->GetTranslation();
	const Double_t *rot = alignment->GetRotationMatrix();
	if( tr && rot ){
	  for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
	  for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
	}
      }
    } else {
      fAliT[0] = 0.;
      fAliT[1] = 0.;
      fAliT[2] = 0.;
      for( int i=0; i<9; i++ ) fAliR[i] = 0;
      fAliR[0] = 1.;
      fAliR[4] = 1.;
      fAliR[8] = 1.;
    }
    fLastSector = sector;
  }
  // alignment->LocalToMaster( x, y);

  XYZ[0] = fAliT[0] + x[0]*fAliR[0] + x[1]*fAliR[1] + x[2]*fAliR[2];
  XYZ[1] = fAliT[1] + x[0]*fAliR[3] + x[1]*fAliR[4] + x[2]*fAliR[5];
  XYZ[2] = fAliT[2] + x[0]*fAliR[6] + x[1]*fAliR[7] + x[2]*fAliR[8];

  return 0; 
}

int  AliHLTTPCClusterTransformation::ReverseAlignment( float XYZ[], int slice, int padrow)
{
  // reverse the alignment correction
  Int_t sector=-99, thisrow=-99;
  AliHLTTPCTransform::Slice2Sector( slice, padrow, sector, thisrow);
  if( sector!= fLastSector ){
    if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
      TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
      if ( alignment ){
	const Double_t *tr = alignment->GetTranslation();
	const Double_t *rot = alignment->GetRotationMatrix();
	if(tr){
	  for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
	}
	SetRotationMatrix(rot, true);
      }
    } else {
      fAliT[0] = 0.;
      fAliT[1] = 0.;
      fAliT[2] = 0.;
      SetRotationMatrix(NULL, true);
    }
    fLastSector = sector;
  }

  // correct for alignment: translation
  float xyz[3];
  xyz[0] = XYZ[0] - fAliT[0];
  xyz[1] = XYZ[1] - fAliT[1];
  xyz[2] = XYZ[2] - fAliT[2];

  // correct for alignment: rotation
  XYZ[0]=xyz[0]*fAdjR[0] + xyz[1]*fAdjR[1] + xyz[2]*fAdjR[2];
  XYZ[1]=xyz[0]*fAdjR[3] + xyz[1]*fAdjR[4] + xyz[2]*fAdjR[5];
  XYZ[2]=xyz[0]*fAdjR[6] + xyz[1]*fAdjR[7] + xyz[2]*fAdjR[8];

  return 0;
}

void AliHLTTPCClusterTransformation::SetRotationMatrix(const Double_t *rot, bool bCalcAdjugate)
{
  // set the rotation matrix and calculate the adjugate if requested
  if (rot) {
    for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
    if (bCalcAdjugate) {
      CalcAdjugateRotation();
    }
    return;
  }
  for( int i=0; i<9; i++ ) {fAliR[i] = 0; fAdjR[i] = 0;}
  fAliR[0] = 1.;
  fAliR[4] = 1.;
  fAliR[8] = 1.; 
  fAdjR[0] = 1.;
  fAdjR[4] = 1.;
  fAdjR[8] = 1.; 
}

bool AliHLTTPCClusterTransformation::CalcAdjugateRotation(bool bCheck)
{
  // check rotation matrix and adjugate for consistency
  fAdjR[0]= fAliR[4]*fAliR[8]-fAliR[5]*fAliR[7];
  fAdjR[1]= fAliR[5]*fAliR[6]-fAliR[3]*fAliR[8];
  fAdjR[2]= fAliR[3]*fAliR[7]-fAliR[4]*fAliR[6];

  fAdjR[3]= fAliR[2]*fAliR[7]-fAliR[1]*fAliR[8];
  fAdjR[4]= fAliR[0]*fAliR[8]-fAliR[2]*fAliR[6];
  fAdjR[5]= fAliR[2]*fAliR[6]-fAliR[0]*fAliR[7];

  fAdjR[6]= fAliR[1]*fAliR[5]-fAliR[2]*fAliR[4];
  fAdjR[7]= fAliR[2]*fAliR[3]-fAliR[0]*fAliR[5];
  fAdjR[8]= fAliR[0]*fAliR[4]-fAliR[1]*fAliR[3];

  if (bCheck) {
    for (int r=0; r<3; r++) {
      for (int c=0; c<3; c++) {
	float a=0.;
	float expected=0.;
	if (r==c) expected=1.;
	for (int i=0; i<3; i++) {
	  a+=fAliR[3*r+i]*fAdjR[c+(3*i)];
	}
	if (TMath::Abs(a-expected)>0.00001) {
	  std::cout << "inconsistent adjugate at " << r << c << ": " << a << " " << expected << std::endl;
	  return false;
	}
      }
    }
  }
  return true;
}

void AliHLTTPCClusterTransformation::Print(const char* /*option*/) const
{
  // print info
  ios::fmtflags coutflags=std::cout.flags(); // backup cout status flags
  std::cout << "AliHLTTPCClusterTransformation for sector " << fLastSector << std::endl;

  std::cout.setf(ios_base::showpos|ios_base::showpos|ios::right);
  std::cout << "  translation: " << std::endl;
  int r=0;
  for (r=0; r<3; r++) {
    std::cout << setw(7) << fixed << setprecision(2);
    cout << "  " << fAliT[r] << std::endl;
  }
  std::cout << "  rotation and adjugated rotation: " << std::endl;
  for (r=0; r<3; r++) {
    int c=0;
    std::cout << setw(7) << fixed << setprecision(2);
    for (c=0; c<3; c++) std::cout << "  " << fAliR[3*r+c];
    std::cout << "      ";
    for (c=0; c<3; c++) std::cout << "  " << fAdjR[3*r+c];
    std::cout << endl;
  }
  std::cout.flags(coutflags); // restore the original flags
}
Commit	Line	Data
0bbbbd17	1	// $Id: AliHLTTPCClusterTransformation.cxx 41244 2010-05-14 08:13:35Z kkanaki $
	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: Kalliopi Kanaki <Kalliopi.Kanaki@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 AliHLTTPCClusterTransformation.cxx
	20	@author Kalliopi Kanaki, Sergey Gorbubnov
	21	@date
	22	@brief
	23	*/
	24
	25
	26	#include "AliHLTTPCClusterTransformation.h"
	27	#include "AliHLTTPCTransform.h"
	28
	29	#include "AliTPCcalibDB.h"
	30	#include "AliTPCTransform.h"
0bbbbd17	31	#include "AliTPCParam.h"
67e7ad4f	32	#include "AliTPCRecoParam.h"
8f2bd9bd	33	#include "AliGeomManager.h"
48fd12f9	34	#include <iostream>
	35	#include <iomanip>
	36
	37	using namespace std;
0bbbbd17	38
	39	ClassImp(AliHLTTPCClusterTransformation) //ROOT macro for the implementation of ROOT specific class methods
	40
	41	AliHLTTPCClusterTransformation::AliHLTTPCClusterTransformation()
	42	:
	43	fOfflineTransform(NULL),
85b28c8f	44	fOfflineTPCParam( NULL ),
85b28c8f	45	fLastSector(-1)
0bbbbd17	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
85b28c8f	52	fAliT[0] = 0.;
	53	fAliT[1] = 0.;
	54	fAliT[2] = 0.;
48fd12f9	55	SetRotationMatrix();
0bbbbd17	56	}
0bbbbd17	57	AliHLTTPCClusterTransformation::~AliHLTTPCClusterTransformation()
	58	{
	59	// see header file for class documentation
d039e59e	60	//delete fOfflineTransform;
0bbbbd17	61	}
	62
	63
dccee396	64	int AliHLTTPCClusterTransformation::Init( double FieldBz, UInt_t TimeStamp )
0bbbbd17	65	{
	66	// Initialisation
	67
d039e59e	68	//delete fOfflineTransform;
d039e59e	69	fOfflineTransform = 0;
67e7ad4f	70	fOfflineTPCParam = 0;
67e7ad4f	71
0bbbbd17	72	AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();
	73
	74	if(!pCalib ) return -1;
	75
	76	pCalib->SetExBField(FieldBz);
8f2bd9bd	77
	78	if(!AliGeomManager::GetGeometry()){
	79	AliGeomManager::LoadGeometry();
	80	}
0bbbbd17	81
67e7ad4f	82	if( !pCalib->GetTransform() ) return -2;
67e7ad4f	83
d039e59e	84	//fOfflineTransform = new AliTPCTransform (*pCalib->GetTransform());
d039e59e	85	fOfflineTransform = pCalib->GetTransform();
67e7ad4f	86	fOfflineTransform->SetCurrentRecoParam( AliTPCRecoParam::GetHLTParam() );
dccee396	87	fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
0bbbbd17	88	fOfflineTPCParam = pCalib->GetParameters();
	89	if( !fOfflineTPCParam ) return -3;
	90
	91	fOfflineTPCParam->Update();
85b28c8f	92	fOfflineTPCParam->ReadGeoMatrices();
	93
	94	fLastSector = -1;
	95
	96	fAliT[0] = 0.;
	97	fAliT[1] = 0.;
	98	fAliT[2] = 0.;
48fd12f9	99	SetRotationMatrix();
0bbbbd17	100
	101	return 0;
	102	}
	103
	104
67e7ad4f	105	void AliHLTTPCClusterTransformation::SetCurrentTimeStamp( UInt_t TimeStamp )
0bbbbd17	106	{
0bbbbd17	107	// Set the current time stamp
67e7ad4f	108	if( fOfflineTransform ) fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
0bbbbd17	109	}
	110
	111
	112	int AliHLTTPCClusterTransformation::Transform( int Slice, int Row, float Pad, float Time, float XYZ[] )
	113	{
	114	// Convert row, pad, time to X Y Z
	115
	116	Int_t sector=-99, thisrow=-99;
	117	AliHLTTPCTransform::Slice2Sector( Slice, Row, sector, thisrow);
	118
	119	if( !fOfflineTransform ){
	120	AliHLTTPCTransform::Raw2Local( XYZ, sector, thisrow, Pad, Time);
	121	if(Slice>17) XYZ[1]= - XYZ[1];
	122	return 0;
	123	}
	124
	125	Int_t iSector[1]= {sector};
	126	Double_t x[3] = { thisrow, Pad, Time };
	127	fOfflineTransform->Transform(x,iSector,0,1);
0bbbbd17	128
85b28c8f	129	if( sector!= fLastSector ){
	130	if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
	131	TGeoHMatrix *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
	132	if ( alignment ){
	133	const Double_t *tr = alignment->GetTranslation();
	134	const Double_t *rot = alignment->GetRotationMatrix();
	135	if( tr && rot ){
	136	for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
	137	for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
	138	}
	139	}
	140	} else {
	141	fAliT[0] = 0.;
	142	fAliT[1] = 0.;
	143	fAliT[2] = 0.;
	144	for( int i=0; i<9; i++ ) fAliR[i] = 0;
	145	fAliR[0] = 1.;
	146	fAliR[4] = 1.;
	147	fAliR[8] = 1.;
	148	}
	149	fLastSector = sector;
	150	}
	151	// alignment->LocalToMaster( x, y);
	152
	153	XYZ[0] = fAliT[0] + x[0]fAliR[0] + x[1]fAliR[1] + x[2]*fAliR[2];
	154	XYZ[1] = fAliT[1] + x[0]fAliR[3] + x[1]fAliR[4] + x[2]*fAliR[5];
	155	XYZ[2] = fAliT[2] + x[0]fAliR[6] + x[1]fAliR[7] + x[2]*fAliR[8];
0bbbbd17	156
	157	return 0;
	158	}
48fd12f9	159
	160	int AliHLTTPCClusterTransformation::ReverseAlignment( float XYZ[], int slice, int padrow)
	161	{
	162	// reverse the alignment correction
	163	Int_t sector=-99, thisrow=-99;
	164	AliHLTTPCTransform::Slice2Sector( slice, padrow, sector, thisrow);
	165	if( sector!= fLastSector ){
	166	if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
	167	TGeoHMatrix *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
	168	if ( alignment ){
	169	const Double_t *tr = alignment->GetTranslation();
	170	const Double_t *rot = alignment->GetRotationMatrix();
	171	if(tr){
	172	for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
	173	}
	174	SetRotationMatrix(rot, true);
	175	}
	176	} else {
	177	fAliT[0] = 0.;
	178	fAliT[1] = 0.;
	179	fAliT[2] = 0.;
	180	SetRotationMatrix(NULL, true);
	181	}
	182	fLastSector = sector;
	183	}
	184
	185	// correct for alignment: translation
	186	float xyz[3];
	187	xyz[0] = XYZ[0] - fAliT[0];
	188	xyz[1] = XYZ[1] - fAliT[1];
	189	xyz[2] = XYZ[2] - fAliT[2];
	190
	191	// correct for alignment: rotation
	192	XYZ[0]=xyz[0]fAdjR[0] + xyz[1]fAdjR[1] + xyz[2]*fAdjR[2];
	193	XYZ[1]=xyz[0]fAdjR[3] + xyz[1]fAdjR[4] + xyz[2]*fAdjR[5];
	194	XYZ[2]=xyz[0]fAdjR[6] + xyz[1]fAdjR[7] + xyz[2]*fAdjR[8];
	195
	196	return 0;
	197	}
	198
	199	void AliHLTTPCClusterTransformation::SetRotationMatrix(const Double_t *rot, bool bCalcAdjugate)
	200	{
	201	// set the rotation matrix and calculate the adjugate if requested
	202	if (rot) {
	203	for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
	204	if (bCalcAdjugate) {
	205	CalcAdjugateRotation();
	206	}
	207	return;
	208	}
	209	for( int i=0; i<9; i++ ) {fAliR[i] = 0; fAdjR[i] = 0;}
	210	fAliR[0] = 1.;
	211	fAliR[4] = 1.;
	212	fAliR[8] = 1.;
	213	fAdjR[0] = 1.;
	214	fAdjR[4] = 1.;
	215	fAdjR[8] = 1.;
	216	}
	217
	218	bool AliHLTTPCClusterTransformation::CalcAdjugateRotation(bool bCheck)
	219	{
	220	// check rotation matrix and adjugate for consistency
	221	fAdjR[0]= fAliR[4]fAliR[8]-fAliR[5]fAliR[7];
	222	fAdjR[1]= fAliR[5]fAliR[6]-fAliR[3]fAliR[8];
223	fAdjR[2]= fAliR[3]fAliR[7]-fAliR[4]fAliR[6];
224
225	fAdjR[3]= fAliR[2]fAliR[7]-fAliR[1]fAliR[8];
226	fAdjR[4]= fAliR[0]fAliR[8]-fAliR[2]fAliR[6];
227	fAdjR[5]= fAliR[2]fAliR[6]-fAliR[0]fAliR[7];
228
229	fAdjR[6]= fAliR[1]fAliR[5]-fAliR[2]fAliR[4];
230	fAdjR[7]= fAliR[2]fAliR[3]-fAliR[0]fAliR[5];
231	fAdjR[8]= fAliR[0]fAliR[4]-fAliR[1]fAliR[3];
232
233	if (bCheck) {
234	for (int r=0; r<3; r++) {
235	for (int c=0; c<3; c++) {
236	float a=0.;
237	float expected=0.;
238	if (r==c) expected=1.;
239	for (int i=0; i<3; i++) {
240	a+=fAliR[3r+i]fAdjR[c+(3*i)];
241	}
242	if (TMath::Abs(a-expected)>0.00001) {
243	std::cout << "inconsistent adjugate at " << r << c << ": " << a << " " << expected << std::endl;
244	return false;
245	}
246	}
247	}
248	}
249	return true;
250	}
251
252	void AliHLTTPCClusterTransformation::Print(const char* /option/) const
253	{
254	// print info
255	ios::fmtflags coutflags=std::cout.flags(); // backup cout status flags
256	std::cout << "AliHLTTPCClusterTransformation for sector " << fLastSector << std::endl;
257
258	std::cout.setf(ios_base::showpos\|ios_base::showpos\|ios::right);
259	std::cout << " translation: " << std::endl;
260	int r=0;
261	for (r=0; r<3; r++) {
262	std::cout << setw(7) << fixed << setprecision(2);
263	cout << " " << fAliT[r] << std::endl;
264	}
265	std::cout << " rotation and adjugated rotation: " << std::endl;
266	for (r=0; r<3; r++) {
267	int c=0;
268	std::cout << setw(7) << fixed << setprecision(2);
269	for (c=0; c<3; c++) std::cout << " " << fAliR[3*r+c];
270	std::cout << " ";
271	for (c=0; c<3; c++) std::cout << " " << fAdjR[3*r+c];
272	std::cout << endl;
273	}
274	std::cout.flags(coutflags); // restore the original flags
275	}