[u/mrichter/AliRoot.git] / TPC / AliTPCAlign.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.1  2007/10/01 14:12:45  kowal2
Class creating the aligmnent object fro the surveyor measurements.

*/ 

//
//  Creates the TPC align object
//

#include "AliTPCAlign.h"
//
#include "TROOT.h"
#include "Riostream.h"
#include "TFile.h"
#include "TMath.h"
#include "TSystem.h"
#include "AliSurveyObj.h"
#include "AliAlignObjParams.h"
#include "AliCDBStorage.h"
#include <TClonesArray.h>
#include <TFile.h>
#include "AliLog.h"
#include "AliCDBManager.h"

ClassImp(AliTPCAlign)

AliTPCAlign::AliTPCAlign() :
  TObject(),
  fFileLoc(0x0),
  fFileGlob(0x0),
  fTPCAlignObj(0x0),
  fX(),
  fA(),
  fY(),
  fDebug(0)
{
  //
  //  default constructor
  //
}   
//________________________________________________________________________
AliTPCAlign::AliTPCAlign(Int_t reportloc, Int_t reportglob) :
  TObject(),
  fFileLoc(0x0),
  fFileGlob(0x0),
  fTPCAlignObj(0x0),
  fX(6,1),
  fA(24,6),
  fY(24,1),
  fDebug(0)
{
  //
  // constructor - defines data files
  //
  fFileLoc = new Char_t[80];
  fFileGlob = new Char_t[80];
  Char_t path[50];
  sprintf(path,gSystem->Getenv("ALICE_ROOT")); 
  //
  sprintf(fFileLoc,"%s/TPC/Survey_%d_TPC.txt",path,reportloc);
  sprintf(fFileGlob,"%s/TPC/Survey_%d_TPC.txt",path,reportglob);
  //

}
//_________________________________________________________________________
AliTPCAlign::AliTPCAlign(const AliTPCAlign &align) :
  TObject(),
  fFileLoc(0x0),
  fFileGlob(0x0),
  fTPCAlignObj(0x0),
  fX(),
  fA(),
  fY(),
  fDebug(0)
{
  //
  //  copy constructor - dummy
  //
  fDebug = align.fDebug;
}
//__________________________________________________________________________
AliTPCAlign & AliTPCAlign::operator =(const AliTPCAlign & align)
{
  //
  // assignment operator - dummy
  //
  fDebug=align.fDebug;
  return (*this);
}

//__________________________________________________________________________
AliTPCAlign::~AliTPCAlign(){
  //
  // destructor
  //
  if(fTPCAlignObj) delete fTPCAlignObj;
}
//__________________________________________________________________________
Bool_t AliTPCAlign::LoadSurveyData(){
  //
  // for a time being it loads from the local file the surveyed point
  // and has the ideal points hardwired. I am waiting until Ricardo
  // completes his job
  // 

 AliSurveyObj * s1 = new AliSurveyObj();
 s1->FillFromLocalFile(fFileGlob);
 //
 Int_t numberPoints = 8;
 //
 TString pointNames[8] = {"T1Final_R04241","T1Final_R05241","T1Final_R06241",
                          "T1Final_R07241","T1Final_R08241","T1Final_R10241",
                          "T1Final_R11241","T1Final_R12241"};
 //
 Float_t surveyedPoints[8][3];
 AliSurveyPoint *currPoint;
 //
 for(Int_t i=0;i<numberPoints;i++){
   currPoint=0;
   currPoint = (AliSurveyPoint *) s1->GetData()->FindObject(pointNames[i]);
   //
   if(currPoint){
     surveyedPoints[i][0]=currPoint->GetX();
     surveyedPoints[i][1]=currPoint->GetY();
     surveyedPoints[i][2]=currPoint->GetZ();
     if(fDebug)
     Printf(Form("INFO: Point \"%s\" coordinates read.", pointNames[i].Data()));
   }
   else {
     if(fDebug){
    Printf(Form("ERROR: Essential point missing: \"%s\"", pointNames[i].Data()));
    return 1;
     }
   }  
 }
 //
 //  Ideal points
 //
 Float_t idealPoints[8][3];
 //
 idealPoints[0][0]=0.7973;
 idealPoints[0][1]=-2.5278;
 idealPoints[0][2]=3.0165;
 //
 idealPoints[1][0]=2.5873;
 idealPoints[1][1]=-0.5738;
 idealPoints[1][2]=3.0166;
 //
 idealPoints[2][0]=1.7898;
 idealPoints[2][1]=1.9540;
 idealPoints[2][2]=3.0164;
 //
 idealPoints[3][0]=-0.2237;
 idealPoints[3][1]=0.7055;
 idealPoints[3][2]=3.0161;
 //
 idealPoints[4][0]=-0.7230;
 idealPoints[4][1]=0.1598;
 idealPoints[4][2]=3.0162;
//
 idealPoints[5][0]=0.2225;
 idealPoints[5][1]=-0.7058;
 idealPoints[5][2]=3.0166;
 //
 idealPoints[6][0]=0.7224;
 idealPoints[6][1]=-0.1608;
 idealPoints[6][2]=3.0164;
 //
 idealPoints[7][0]=0.4999;
 idealPoints[7][1]=0.5452;
 idealPoints[7][2]=3.0164;
 //
 // Create and fill matrices a & y
 //
 for(Int_t i = 0;i<numberPoints;i++){
   for(Int_t j=0;j<3;j++){
     fY(i*3+j,0)=surveyedPoints[i][j]-idealPoints[i][j];

   }
 }
 fA.Zero(); 
 //
 //
 // setting matrix a
 //
 for(Int_t i=0;i<numberPoints;i++){
   fA(3*i,0)=    -idealPoints[i][1];
   fA(3*i,1)=    idealPoints[i][2];
   fA(3*i,3)=1.; 
   fA(3*i+1,0)=    idealPoints[i][0];
   fA(3*i+1,2)=    -idealPoints[i][2];
   fA(3*i+1,4)=    1.;
   fA(3*i+2,1)=   -idealPoints[i][0];
   fA(3*i+2,2)=   idealPoints[i][1];
   fA(3*i+2,5)=1.; 
 }
 //
 delete s1;
 //
 return 0;
}
//_________________________________________________________________
Double_t AliTPCAlign::ComputeTransform(){
  //
  // Here simple matrix operations for the linear least square
  // The trigonometric function sin is approximated by the angle
  // and the cos by "1", because angles are very small (Y-convention!)
  // This secures the linearity of the problem
  // This method returns a sum of squares of residuals
  //
  TMatrixD tt1(TMatrixD::kInverted,(TMatrixD(fA,TMatrixD::kTransposeMult,fA)));
  fX=TMatrixD(tt1,TMatrixD::kMult,TMatrixD(fA,TMatrixD::kTransposeMult,fY));
  //
 TMatrixD Unit(24,24);
 Unit.UnitMatrix();
 TMatrixD xxt1(TMatrixD::kInverted,TMatrixD(fA,TMatrixD::kTransposeMult,fA));
 TMatrixD t(fA,TMatrixD::kMult,TMatrixD(xxt1,TMatrixD::kMultTranspose,fA));
 TMatrixD t2(Unit,TMatrixD::kMinus,t);
 TMatrixD chi2(fY,TMatrixD::kTransposeMult,TMatrixD(t2,TMatrixD::kMult,fY));
 //

 return chi2(0,0);
 	    
}
//_______________________________________________________________________
void AliTPCAlign::CreateAlignObj(){
  //
  // This method creates AliAlignObj and fills it with Euler angles
  // and shifts. The units are degrees and cm.
  // 
 fTPCAlignObj = new AliAlignObjParams();
 fTPCAlignObj->SetSymName("ALIC_1/TPC_M_1");
 fTPCAlignObj->SetVolUID(0);
 Double_t raddeg = TMath::RadToDeg();
 Double_t phi,theta,psi;
 //
 phi=fX(0,0)*raddeg;
 theta=fX(1,0)*raddeg;
 psi=fX(2,0)*raddeg;
 //
 Double_t dx,dy,dz;
 dx=fX(3,0)*100.;
 dy=fX(4,0)*100.;
 dz=fX(5,0)*100.;
 fTPCAlignObj->SetRotation(psi,theta,phi);
 fTPCAlignObj->SetTranslation(dx,dy,dz);
 fTPCAlignObj->Print("");
}
//______________________________________________________________________
void AliTPCAlign::Run(){
  //
  // runs the full chain
  //
  SetDebug(0);
  Bool_t flag = LoadSurveyData();
  if(flag) {
    cout<<"Missing points"<<endl;
    return;
  }
  Double_t chi2 = ComputeTransform();
  if(chi2>0.01) return;
  CreateAlignObj();
  StoreAlignObj();
  //
}
//_________________________________________________________________________
void AliTPCAlign::StoreAlignObj(){
  //
AliCDBManager* cdb = AliCDBManager::Instance();
 if(!cdb->IsDefaultStorageSet()) cdb->SetDefaultStorage("local://$ALICE_ROOT");
 //
TClonesArray *array = new TClonesArray("AliAlignObjParams",1);
//
 Double_t shifts[3];
 Double_t rots[3];
 //
 fTPCAlignObj->GetPars(shifts,rots);
 new((*array)[0]) AliAlignObjParams(fTPCAlignObj->GetSymName(),0,shifts[0],
                   shifts[1],shifts[2],rots[0],rots[1],rots[2],kTRUE);

//
// storing either in the OCDB or local file
//
  if( TString(gSystem->Getenv("TOCDB")) != TString("kTRUE") ){
    // save on file
    const char* filename = "TPCSurveyMisalignment.root";
    Char_t fullname[80];
    sprintf(fullname,"%s/TPC/%s",gSystem->Getenv("ALICE_ROOT"),filename);
    TFile *f = new TFile(fullname,"RECREATE");
    if(!f){
      AliError("cannot open file for output\n");
      return;
    }
    AliInfo(Form("Saving alignment objects to the file %s", filename));
    f->cd();
    f->WriteObject(array,"TPCAlignObjs","kSingleKey");
    f->Close();
  }else{
    // save in CDB storage
    AliCDBStorage* storage;
    //
   TString Storage = gSystem->Getenv("STORAGE");
    if(!Storage.BeginsWith("local://") && !Storage.BeginsWith("alien://")) {
      AliError(Form(
      "STORAGE variable set to %s is not valid. Exiting\n",Storage.Data()));
      return;
    }
    storage = cdb->GetStorage(Storage.Data());
    if(!storage){
      AliError(Form("Unable to open storage %s\n",Storage.Data()));
      return;
    }
    //
    AliCDBMetaData* md = new AliCDBMetaData();
    md->SetResponsible("Marek Kowalski");
    md->SetComment("Full misalignment of entire TPC from surveyors");
    AliCDBId id("TPC/Align/Data",0,AliCDBRunRange::Infinity());
    storage->Put(array,id,md);
  }

}
Commit	Line	Data
2574314d	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	$Log$
dd89d4d2	18	Revision 1.1 2007/10/01 14:12:45 kowal2
	19	Class creating the aligmnent object fro the surveyor measurements.
	20
2574314d	21	*/
	22
	23	//
	24	// Creates the TPC align object
	25	//
	26
	27	#include "AliTPCAlign.h"
	28	//
	29	#include "TROOT.h"
	30	#include "Riostream.h"
	31	#include "TFile.h"
	32	#include "TMath.h"
	33	#include "TSystem.h"
	34	#include "AliSurveyObj.h"
	35	#include "AliAlignObjParams.h"
	36	#include "AliCDBStorage.h"
	37	#include <TClonesArray.h>
	38	#include <TFile.h>
	39	#include "AliLog.h"
	40	#include "AliCDBManager.h"
	41
dd89d4d2	42	ClassImp(AliTPCAlign)
2574314d	43
	44	AliTPCAlign::AliTPCAlign() :
	45	TObject(),
	46	fFileLoc(0x0),
	47	fFileGlob(0x0),
	48	fTPCAlignObj(0x0),
	49	fX(),
	50	fA(),
	51	fY(),
	52	fDebug(0)
	53	{
	54	//
	55	// default constructor
	56	//
	57	}
	58	//________________________________________________________________________
	59	AliTPCAlign::AliTPCAlign(Int_t reportloc, Int_t reportglob) :
	60	TObject(),
	61	fFileLoc(0x0),
	62	fFileGlob(0x0),
	63	fTPCAlignObj(0x0),
	64	fX(6,1),
	65	fA(24,6),
	66	fY(24,1),
	67	fDebug(0)
	68	{
	69	//
	70	// constructor - defines data files
	71	//
	72	fFileLoc = new Char_t[80];
	73	fFileGlob = new Char_t[80];
	74	Char_t path[50];
	75	sprintf(path,gSystem->Getenv("ALICE_ROOT"));
	76	//
	77	sprintf(fFileLoc,"%s/TPC/Survey_%d_TPC.txt",path,reportloc);
	78	sprintf(fFileGlob,"%s/TPC/Survey_%d_TPC.txt",path,reportglob);
	79	//
	80
	81	}
	82	//_________________________________________________________________________
	83	AliTPCAlign::AliTPCAlign(const AliTPCAlign &align) :
	84	TObject(),
	85	fFileLoc(0x0),
	86	fFileGlob(0x0),
	87	fTPCAlignObj(0x0),
	88	fX(),
	89	fA(),
	90	fY(),
	91	fDebug(0)
	92	{
	93	//
	94	// copy constructor - dummy
	95	//
	96	fDebug = align.fDebug;
	97	}
	98	//__________________________________________________________________________
	99	AliTPCAlign & AliTPCAlign::operator =(const AliTPCAlign & align)
	100	{
	101	//
	102	// assignment operator - dummy
	103	//
	104	fDebug=align.fDebug;
	105	return (*this);
	106	}
107
108	//__________________________________________________________________________
109	AliTPCAlign::~AliTPCAlign(){
110	//
111	// destructor
112	//
113	if(fTPCAlignObj) delete fTPCAlignObj;
114	}
115	//__________________________________________________________________________
116	Bool_t AliTPCAlign::LoadSurveyData(){
117	//
118	// for a time being it loads from the local file the surveyed point
119	// and has the ideal points hardwired. I am waiting until Ricardo
120	// completes his job
121	//
122
123	AliSurveyObj * s1 = new AliSurveyObj();
124	s1->FillFromLocalFile(fFileGlob);
125	//
126	Int_t numberPoints = 8;
127	//
128	TString pointNames[8] = {"T1Final_R04241","T1Final_R05241","T1Final_R06241",
129	"T1Final_R07241","T1Final_R08241","T1Final_R10241",
130	"T1Final_R11241","T1Final_R12241"};
131	//
132	Float_t surveyedPoints[8][3];
133	AliSurveyPoint *currPoint;
134	//
135	for(Int_t i=0;i<numberPoints;i++){
136	currPoint=0;
137	currPoint = (AliSurveyPoint *) s1->GetData()->FindObject(pointNames[i]);
138	//
139	if(currPoint){
140	surveyedPoints[i][0]=currPoint->GetX();
141	surveyedPoints[i][1]=currPoint->GetY();
142	surveyedPoints[i][2]=currPoint->GetZ();
143	if(fDebug)
144	Printf(Form("INFO: Point \"%s\" coordinates read.", pointNames[i].Data()));
145	}
146	else {
147	if(fDebug){
148	Printf(Form("ERROR: Essential point missing: \"%s\"", pointNames[i].Data()));
149	return 1;
150	}
151	}
152	}
153	//
154	// Ideal points
155	//
156	Float_t idealPoints[8][3];
157	//
158	idealPoints[0][0]=0.7973;
159	idealPoints[0][1]=-2.5278;
160	idealPoints[0][2]=3.0165;
161	//
162	idealPoints[1][0]=2.5873;
163	idealPoints[1][1]=-0.5738;
164	idealPoints[1][2]=3.0166;
165	//
166	idealPoints[2][0]=1.7898;
167	idealPoints[2][1]=1.9540;
168	idealPoints[2][2]=3.0164;
169	//
170	idealPoints[3][0]=-0.2237;
171	idealPoints[3][1]=0.7055;
172	idealPoints[3][2]=3.0161;
173	//
174	idealPoints[4][0]=-0.7230;
175	idealPoints[4][1]=0.1598;
176	idealPoints[4][2]=3.0162;
177	//
178	idealPoints[5][0]=0.2225;
179	idealPoints[5][1]=-0.7058;
180	idealPoints[5][2]=3.0166;
181	//
182	idealPoints[6][0]=0.7224;
183	idealPoints[6][1]=-0.1608;
184	idealPoints[6][2]=3.0164;
185	//
186	idealPoints[7][0]=0.4999;
187	idealPoints[7][1]=0.5452;
188	idealPoints[7][2]=3.0164;
189	//
190	// Create and fill matrices a & y
191	//
192	for(Int_t i = 0;i<numberPoints;i++){
193	for(Int_t j=0;j<3;j++){
194	fY(i*3+j,0)=surveyedPoints[i][j]-idealPoints[i][j];
195
196	}
197	}
198	fA.Zero();
199	//
200	//
201	// setting matrix a
202	//
203	for(Int_t i=0;i<numberPoints;i++){
204	fA(3*i,0)= -idealPoints[i][1];
205	fA(3*i,1)= idealPoints[i][2];
206	fA(3*i,3)=1.;
207	fA(3*i+1,0)= idealPoints[i][0];
208	fA(3*i+1,2)= -idealPoints[i][2];
209	fA(3*i+1,4)= 1.;
210	fA(3*i+2,1)= -idealPoints[i][0];
211	fA(3*i+2,2)= idealPoints[i][1];
212	fA(3*i+2,5)=1.;
213	}
214	//
215	delete s1;
216	//
217	return 0;
218	}
219	//_________________________________________________________________
220	Double_t AliTPCAlign::ComputeTransform(){
221	//
222	// Here simple matrix operations for the linear least square
223	// The trigonometric function sin is approximated by the angle
224	// and the cos by "1", because angles are very small (Y-convention!)
225	// This secures the linearity of the problem
226	// This method returns a sum of squares of residuals
227	//
228	TMatrixD tt1(TMatrixD::kInverted,(TMatrixD(fA,TMatrixD::kTransposeMult,fA)));
229	fX=TMatrixD(tt1,TMatrixD::kMult,TMatrixD(fA,TMatrixD::kTransposeMult,fY));
230	//
231	TMatrixD Unit(24,24);
232	Unit.UnitMatrix();
233	TMatrixD xxt1(TMatrixD::kInverted,TMatrixD(fA,TMatrixD::kTransposeMult,fA));
234	TMatrixD t(fA,TMatrixD::kMult,TMatrixD(xxt1,TMatrixD::kMultTranspose,fA));
235	TMatrixD t2(Unit,TMatrixD::kMinus,t);
236	TMatrixD chi2(fY,TMatrixD::kTransposeMult,TMatrixD(t2,TMatrixD::kMult,fY));
237	//
238
239	return chi2(0,0);
240
241	}
242	//_______________________________________________________________________
243	void AliTPCAlign::CreateAlignObj(){
244	//
245	// This method creates AliAlignObj and fills it with Euler angles
246	// and shifts. The units are degrees and cm.
247	//
248	fTPCAlignObj = new AliAlignObjParams();
249	fTPCAlignObj->SetSymName("ALIC_1/TPC_M_1");
250	fTPCAlignObj->SetVolUID(0);
251	Double_t raddeg = TMath::RadToDeg();
252	Double_t phi,theta,psi;
253	//
254	phi=fX(0,0)*raddeg;
255	theta=fX(1,0)*raddeg;
256	psi=fX(2,0)*raddeg;
257	//
258	Double_t dx,dy,dz;
259	dx=fX(3,0)*100.;
260	dy=fX(4,0)*100.;
261	dz=fX(5,0)*100.;
262	fTPCAlignObj->SetRotation(psi,theta,phi);
263	fTPCAlignObj->SetTranslation(dx,dy,dz);
264	fTPCAlignObj->Print("");
265	}
266	//______________________________________________________________________
267	void AliTPCAlign::Run(){
268	//
269	// runs the full chain
270	//
271	SetDebug(0);
272	Bool_t flag = LoadSurveyData();
273	if(flag) {
274	cout<<"Missing points"<<endl;
275	return;
276	}
277	Double_t chi2 = ComputeTransform();
278	if(chi2>0.01) return;
279	CreateAlignObj();
280	StoreAlignObj();
281	//
282	}
283	//_________________________________________________________________________
284	void AliTPCAlign::StoreAlignObj(){
285	//
286	AliCDBManager* cdb = AliCDBManager::Instance();
287	if(!cdb->IsDefaultStorageSet()) cdb->SetDefaultStorage("local://$ALICE_ROOT");
288	//
289	TClonesArray *array = new TClonesArray("AliAlignObjParams",1);
290	//
291	Double_t shifts[3];
292	Double_t rots[3];
293	//
294	fTPCAlignObj->GetPars(shifts,rots);
295	new((*array)[0]) AliAlignObjParams(fTPCAlignObj->GetSymName(),0,shifts[0],
296	shifts[1],shifts[2],rots[0],rots[1],rots[2],kTRUE);
297
298	//
299	// storing either in the OCDB or local file
300	//
301	if( TString(gSystem->Getenv("TOCDB")) != TString("kTRUE") ){
302	// save on file
303	const char* filename = "TPCSurveyMisalignment.root";
304	Char_t fullname[80];
305	sprintf(fullname,"%s/TPC/%s",gSystem->Getenv("ALICE_ROOT"),filename);
306	TFile *f = new TFile(fullname,"RECREATE");
307	if(!f){
308	AliError("cannot open file for output\n");
309	return;
310	}
311	AliInfo(Form("Saving alignment objects to the file %s", filename));
312	f->cd();
313	f->WriteObject(array,"TPCAlignObjs","kSingleKey");
314	f->Close();
315	}else{
316	// save in CDB storage
317	AliCDBStorage* storage;
318	//
319	TString Storage = gSystem->Getenv("STORAGE");
320	if(!Storage.BeginsWith("local://") && !Storage.BeginsWith("alien://")) {
321	AliError(Form(
322	"STORAGE variable set to %s is not valid. Exiting\n",Storage.Data()));
323	return;
324	}
325	storage = cdb->GetStorage(Storage.Data());
326	if(!storage){
327	AliError(Form("Unable to open storage %s\n",Storage.Data()));
328	return;
329	}
330	//
331	AliCDBMetaData* md = new AliCDBMetaData();
332	md->SetResponsible("Marek Kowalski");
333	md->SetComment("Full misalignment of entire TPC from surveyors");
334	AliCDBId id("TPC/Align/Data",0,AliCDBRunRange::Infinity());
335	storage->Put(array,id,md);
336	}
337
338	}