#include <iostream>
#include <fstream>
#include "TMath.h"
#include "TFile.h"
#include "TGeoManager.h"
#include "TGeoPhysicalNode.h"
#include "TClonesArray.h"
#include "TString.h"
#include "TFitter.h"
#include "TMinuit.h"
#include "AliLog.h"
#include "AliAlignObj.h"
#include "AliAlignObjParams.h"
#include "AliCDBManager.h"
#include "AliCDBStorage.h"
#include "AliCDBMetaData.h"
#include "AliCDBEntry.h"
#include "AliSurveyObj.h"
#include "AliSurveyPoint.h"
#include "AliTRDalignment.h"
void trdAlignmentFcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
using std::ostream;
using std::fstream;
ClassImp(AliTRDalignment)
AliTRDalignment::AliTRDalignment()
:TObject()
,fComment()
,fRan(0)
{
SetZero();
for (int i=0; i<18; i++) for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX[i][j][k][l] = 0.0;
fSurveyY[i][j][k][l] = 0.0;
fSurveyZ[i][j][k][l] = 0.0;
fSurveyEX[i][j][k][l] = 0.0;
fSurveyEY[i][j][k][l] = 0.0;
fSurveyEZ[i][j][k][l] = 0.0;
}
double x[2] = {22.5,30.25};
double y[2] = {353.0, -353.0};
double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5};
for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX0[j][k][l] = -TMath::Power(-1,l) * x[k];
fSurveyY0[j][k][l] = y[j];
fSurveyZ0[j][k][l] = z[k];
}
for (int i=0; i<1000; i++) {
fIbuffer[i] = 0;
fDbuffer[i] = 0.0;
}
}
AliTRDalignment::AliTRDalignment(const AliTRDalignment& source)
:TObject(source)
,fComment(source.fComment)
,fRan(source.fRan)
{
for (int i=0; i<18; i++) SetSm(i,source.fSm[i]);
for (int i=0; i<540; i++) SetCh(i,source.fCh[i]);
for (int i=0; i<18; i++) for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
fSurveyEX[i][j][k][l] = source.fSurveyEX[i][j][k][l];
fSurveyEY[i][j][k][l] = source.fSurveyEY[i][j][k][l];
fSurveyEZ[i][j][k][l] = source.fSurveyEZ[i][j][k][l];
}
for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
}
for (int i=0; i<1000; i++) {
fIbuffer[i] = 0;
fDbuffer[i] = 0.0;
}
}
AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source)
{
if (this != &source) {
for (int i = 0; i < 18; i++) SetSm(i,source.fSm[i]);
for (int i = 0; i < 540; i++) SetCh(i,source.fCh[i]);
for (int i=0; i<18; i++) for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
fSurveyEX[i][j][k][l] = source.fSurveyEX[i][j][k][l];
fSurveyEY[i][j][k][l] = source.fSurveyEY[i][j][k][l];
fSurveyEZ[i][j][k][l] = source.fSurveyEZ[i][j][k][l];
}
for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
}
fComment = source.fComment;
}
return *this;
}
AliTRDalignment& AliTRDalignment::operator*=(double fac)
{
for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] *= fac;
for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] *= fac;
return *this;
}
AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source)
{
for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] += source.fSm[i][j];
for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] += source.fCh[i][j];
return *this;
}
AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source)
{
for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) fSm[i][j] -= source.fSm[i][j];
for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) fCh[i][j] -= source.fCh[i][j];
return *this;
}
Bool_t AliTRDalignment::operator==(const AliTRDalignment &source) const
{
Bool_t areEqual = 1;
for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) areEqual &= (fSm[i][j] == source.fSm[i][j]);
for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) areEqual &= (fCh[i][j] == source.fCh[i][j]);
return areEqual;
}
void AliTRDalignment::SetSmZero()
{
memset(&fSm[0][0],0,sizeof(fSm));
}
void AliTRDalignment::SetChZero()
{
memset(&fCh[0][0],0,sizeof(fCh));
}
void AliTRDalignment::SetSmRandom(double a[6])
{
double x[6];
double xmax[6]={999, 0.6, 999, 999, 999, 999};
for (int i = 0; i < 18; i++) {
for (int j = 0; j < 6; j++) {
do {x[j] = fRan.Gaus(0,a[j]);} while (TMath::Abs(x[j]) > xmax[j]);
}
SetSm(i,x);
}
}
void AliTRDalignment::SetChRandom(double a[6])
{
double x[6];
for (int i = 0; i < 540; i++) {
fRan.Rannor(x[0],x[1]);
fRan.Rannor(x[2],x[3]);
fRan.Rannor(x[4],x[5]);
for (int j = 0; j < 6; j++) x[j] *= a[j];
SetCh(i,x);
}
}
void AliTRDalignment::SetSmFull()
{
double a[6];
a[0] = 0.3;
a[1] = 0.3;
a[2] = 0.3;
a[3] = 0.4/1000.0 / TMath::Pi()*180.0;
a[4] = 2.0/1000.0 / TMath::Pi()*180.0;
a[5] = 0.4/1000.0 / TMath::Pi()*180.0;
SetSmRandom(a);
}
void AliTRDalignment::SetChFull()
{
double a[6];
a[0] = 0.1;
a[1] = 0.1;
a[2] = 0.1;
a[3] = 1.0/1000.0 / TMath::Pi()*180.0;
a[4] = 1.0/1000.0 / TMath::Pi()*180.0;
a[5] = 0.7/1000.0 / TMath::Pi()*180.0;
SetChRandom(a);
}
void AliTRDalignment::SetSmResidual()
{
SetSmZero();
}
void AliTRDalignment::SetChResidual()
{
double a[6];
a[0] = 0.002;
a[1] = 0.003;
a[2] = 0.007;
a[3] = 0.3/1000.0 / TMath::Pi()*180.0;
a[4] = 0.3/1000.0 / TMath::Pi()*180.0;
a[5] = 0.1/1000.0 / TMath::Pi()*180.0;
SetChRandom(a);
}
void AliTRDalignment::PrintSm(int i, FILE * const fp) const
{
fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
,i,fSm[i][0],fSm[i][1],fSm[i][2],fSm[i][3],fSm[i][4],fSm[i][5]
,0,GetSmName(i));
}
void AliTRDalignment::PrintCh(int i, FILE * const fp) const
{
fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
,i,fCh[i][0],fCh[i][1],fCh[i][2],fCh[i][3],fCh[i][4],fCh[i][5]
,GetVoi(i),GetChName(i));
}
void AliTRDalignment::ReadAscii(const char * const filename)
{
double x[6];
int volid;
std::string syna;
int j;
fstream fi(filename,fstream::in);
if (!fi) {
AliError(Form("cannot open input file %s",filename));
return;
}
for (int i = 0; i < 18; i++) {
fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
if (j != i) AliError(Form("sm %d expected, %d found",i,j));
if (volid != 0) AliError(Form("sm %d volume id %d expected, %d found",i,0,volid));
std::string symnam = GetSmName(i);
if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data()));
SetSm(i,x);
}
for (int i = 0; i < 540; i++) {
fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
if (j != i) AliError(Form("ch %d expected, %d found",i,j));
if (volid != GetVoi(i)) AliError(Form("ch %d volume id %d expected, %d found",i,GetVoi(i),volid));
std::string symnam = GetChName(i);
if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data()));
SetCh(i,x);
}
fi.close();
}
void AliTRDalignment::ReadCurrentGeo()
{
TGeoPNEntry *pne;
TGeoHMatrix *ideSm[18];
TGeoHMatrix *misSm[18];
for (int i = 0; i < 18; i++) if ((pne = gGeoManager->GetAlignableEntry(GetSmName(i)))) {
TGeoPhysicalNode *node = pne->GetPhysicalNode();
if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
if (!node) continue;
misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
TGeoHMatrix mat(ideSm[i]->Inverse());
mat.Multiply(misSm[i]);
double *tra = mat.GetTranslation();
double *rot = mat.GetRotationMatrix();
double pars[6];
pars[0] = tra[0];
pars[1] = tra[1];
pars[2] = tra[2];
if (TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) AliError("Failed to extract roll-pitch-yall angles!");
double raddeg = TMath::RadToDeg();
pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
pars[4] = raddeg * TMath::ASin(rot[2]);
pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
SetSm(i,pars);
delete ideSm[i];
delete misSm[i];
}
TGeoHMatrix *ideCh[540];
TGeoHMatrix *misCh[540];
for (int i = 0; i < 540; i++) if ((pne = gGeoManager->GetAlignableEntry(GetChName(i)))) {
TGeoPhysicalNode *node = pne->GetPhysicalNode();
if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i)));
if (!node) continue;
misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
TGeoHMatrix mat(ideCh[i]->Inverse());
mat.Multiply(misCh[i]);
double *tra = mat.GetTranslation();
double *rot = mat.GetRotationMatrix();
double pars[6];
pars[0] = tra[0];
pars[1] = tra[1];
pars[2] = tra[2];
if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
AliError("Failed to extract roll-pitch-yall angles!");
return;
}
double raddeg = TMath::RadToDeg();
pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
pars[4] = raddeg * TMath::ASin(rot[2]);
pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
SetCh(i,pars);
delete ideCh[i];
delete misCh[i];
}
return;
}
void AliTRDalignment::ReadRoot(const char * const filename)
{
TFile fi(filename,"READ");
if (fi.IsOpen()) {
TClonesArray *ar = (TClonesArray*) fi.Get("TRDAlignObjs");
ArToNumbers(ar);
fi.Close();
}
else AliError(Form("cannot open input file %s",filename));
return;
}
void AliTRDalignment::ReadDB(const char * const filename)
{
TFile fi(filename,"READ");
if (fi.IsOpen()) {
AliCDBEntry *e = (AliCDBEntry *) fi.Get("AliCDBEntry");
e->PrintMetaData();
fComment.SetString(e->GetMetaData()->GetComment());
TClonesArray *ar = (TClonesArray *) e->GetObject();
ArToNumbers(ar);
fi.Close();
}
else AliError(Form("cannot open input file %s",filename));
return;
}
void AliTRDalignment::ReadDB(const char * const db, const char * const path,
int run, int version, int subversion)
{
AliCDBManager *cdb = AliCDBManager::Instance();
AliCDBStorage *storLoc = cdb->GetStorage(db);
AliCDBEntry *e = storLoc->Get(path,run,version,subversion);
if (e) {
e->PrintMetaData();
fComment.SetString(e->GetMetaData()->GetComment());
TClonesArray *ar = (TClonesArray *) e->GetObject();
ArToNumbers(ar);
}
}
Bool_t AliTRDalignment::DecodeSurveyPointName(TString pna, Int_t &sm, Int_t &iz,
Int_t &ir, Int_t &iphi) {
if (pna(0,6)!="TRD_sm") {
AliError(Form("unexpected point name: %s",pna.Data()));
return kFALSE;
}
sm = atoi(pna(6,2).Data());
iz = -1;
if (pna(8) == 'a') iz=0;
if (pna(8) == 'c') iz=1;
ir = -1;
if (pna(9) == 'l') ir=0;
if (pna(9) == 'h') ir=1;
iphi = -1;
if (pna(10) == '0') iphi = 0;
if (pna(10) == '1') iphi = 1;
if (sm>=0 && sm<18 && iz>=0 && iz<2 && ir>=0 && ir<2 && iphi>=0 && iphi<2) return kTRUE;
AliError(Form("cannot decode point name: %s",pna.Data()));
return kFALSE;
}
void AliTRDalignment::ReadSurveyReport(const char * const filename)
{
fstream in(filename,fstream::in);
if (!in) {
AliError(Form("cannot open input file %s",filename));
return;
}
TString title,date,subdetector,url,version,observations,system,units;
while (1) {
char pee=in.peek();
if (pee==EOF) break;
TString line;
line.ReadLine(in);
if (line.Contains("Title:")) title.ReadLine(in);
if (line.Contains("Date:")) date.ReadLine(in);
if (line.Contains("Subdetector:")) subdetector.ReadLine(in);
if (line.Contains("URL:")) url.ReadLine(in);
if (line.Contains("Version:")) version.ReadLine(in);
if (line.Contains("Observations:")) observations.ReadLine(in);
if (line.Contains("System:")) system.ReadLine(in);
if (line.Contains("Units:")) units.ReadLine(in);
if (line.Contains("Data:")) break;
}
std::cout<<"title .........."<<title<<std::endl;
std::cout<<"date ..........."<<date<<std::endl;
std::cout<<"subdetector ...."<<subdetector<<std::endl;
std::cout<<"url ............"<<url<<std::endl;
std::cout<<"version ........"<<version<<std::endl;
std::cout<<"observations ..."<<observations<<std::endl;
std::cout<<"system ........."<<system<<std::endl;
std::cout<<"units .........."<<units<<std::endl;
if (!subdetector.Contains("TRD")) {
AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
return;
}
double tocm = 0;
if (units.Contains("mm")) tocm = 0.1;
else if (units.Contains("cm")) tocm = 1.0;
else if (units.Contains("m")) tocm = 100.0;
else if (units.Contains("pc")) tocm = 3.24078e-15;
else {
AliError(Form("unexpected units: %s",units.Data()));
return;
}
if (!system.Contains("ALICEPH")) {
AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
return;
}
while (1) {
TString pna;
char type, target;
double x,y,z,precision;
in >> pna >> x >> y >> z >> type >> target >> precision;
if (in.fail()) break;
Int_t i,j,k,l;
if (DecodeSurveyPointName(pna,i,j,k,l)) {
fSurveyX[i][j][k][l] = tocm*x;
fSurveyY[i][j][k][l] = tocm*y;
fSurveyZ[i][j][k][l] = tocm*z;
fSurveyEX[i][j][k][l] = precision/10;
fSurveyEY[i][j][k][l] = precision/10;
fSurveyEZ[i][j][k][l] = precision/10;
printf("decoded %s %02d %d %d %d %8.2f %8.2f %8.2f %6.2f %6.2f %6.2f\n",
pna.Data(), i, j, k, l,
fSurveyX[i][j][k][l], fSurveyY[i][j][k][l], fSurveyZ[i][j][k][l],
fSurveyEX[i][j][k][l], fSurveyEY[i][j][k][l], fSurveyEZ[i][j][k][l]);
} else AliError(Form("cannot decode point name: %s",pna.Data()));
}
in.close();
TString info = "Survey "+title+" "+date+" "+url+" "+version+" "+observations;
info.ReplaceAll("\r","");
fComment.SetString(info.Data());
}
void AliTRDalignment::ReadSurveyReport(const AliSurveyObj * const so)
{
Int_t size = so->GetEntries();
printf("-> %d\n", size);
TString title = so->GetReportTitle();
TString date = so->GetReportDate();
TString subdetector = so->GetDetector();
TString url = so->GetURL();
TString report = so->GetReportNumber();
TString version = so->GetReportVersion();
TString observations = so->GetObservations();
TString system = so->GetCoordSys();
TString units = so->GetUnits();
std::cout<<"title .........."<<title<<std::endl;
std::cout<<"date ..........."<<date<<std::endl;
std::cout<<"subdetector ...."<<subdetector<<std::endl;
std::cout<<"url ............"<<url<<std::endl;
std::cout<<"version ........"<<version<<std::endl;
std::cout<<"observations ..."<<observations<<std::endl;
std::cout<<"system ........."<<system<<std::endl;
std::cout<<"units .........."<<units<<std::endl;
if (!subdetector.Contains("TRD")) {
AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
return;
}
double tocm = 0;
if (units.Contains("mm")) tocm = 0.1;
else if (units.Contains("cm")) tocm = 1.0;
else if (units.Contains("m")) tocm = 100.0;
else if (units.Contains("pc")) tocm = 3.24078e-15;
else {
AliError(Form("unexpected units: %s",units.Data()));
return;
}
if (!system.Contains("ALICEPH")) {
AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
return;
}
TObjArray *points = so->GetData();
for (int ip = 0; ip<points->GetEntries(); ++ip) {
AliSurveyPoint *po = (AliSurveyPoint *) points->At(ip);
TString pna = po->GetPointName();
Int_t i,j,k,l;
if (DecodeSurveyPointName(pna,i,j,k,l)) {
fSurveyX[i][j][k][l] = tocm*po->GetX();
fSurveyY[i][j][k][l] = tocm*po->GetY();
fSurveyZ[i][j][k][l] = tocm*po->GetZ();
fSurveyEX[i][j][k][l] = po->GetPrecisionX()/10;
fSurveyEY[i][j][k][l] = po->GetPrecisionY()/10;
fSurveyEZ[i][j][k][l] = po->GetPrecisionZ()/10;
printf("decoded %s %02d %d %d %d %8.2f %8.2f %8.2f %6.2f %6.2f %6.2f\n",
pna.Data(), i, j, k, l,
fSurveyX[i][j][k][l], fSurveyY[i][j][k][l], fSurveyZ[i][j][k][l],
fSurveyEX[i][j][k][l], fSurveyEY[i][j][k][l], fSurveyEZ[i][j][k][l]);
} else AliError(Form("cannot decode point name: %s",pna.Data()));
}
TString info = "Survey "+title+" "+date+" "+url+" "+report+" "+version+" "+observations;
info.ReplaceAll("\r","");
fComment.SetString(info.Data());
}
double AliTRDalignment::SurveyChi2(int i, const double * const a) {
if (!IsGeoLoaded()) return 0;
printf("Survey of supermodule %d\n",i);
AliAlignObjParams al(GetSmName(i),0,a[0],a[1],a[2],a[3],a[4],a[5],0);
TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetSmName(i));
if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
TGeoPhysicalNode *node = pne->GetPhysicalNode();
if (!node) {
AliWarning(Form("physical node entry %s has no physical node; making a new one",GetSmName(i)));
node = gGeoManager->MakeAlignablePN(pne);
}
TGeoHMatrix *ma = new TGeoHMatrix();
al.GetLocalMatrix(*ma);
ma->MultiplyLeft(node->GetMatrix());
double chi2=0;
printf(" sm z r phi x (lab phi) y (lab z) z (lab r) all in cm\n");
for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
if (fSurveyEX[i][j][k][l] == 0.0
&& fSurveyEY[i][j][k][l] == 0.0
&& fSurveyEZ[i][j][k][l] == 0.0) continue;
double master[3] = {fSurveyX[i][j][k][l],fSurveyY[i][j][k][l],fSurveyZ[i][j][k][l]};
double local[3];
ma->MasterToLocal(master,local);
double dx = local[0]-fSurveyX0[j][k][l];
double dy = local[1]-fSurveyY0[j][k][l];
double dz = local[2]-fSurveyZ0[j][k][l];
chi2 += dx*dx/fSurveyEX[i][j][k][l]/fSurveyEX[i][j][k][l];
chi2 += dy*dy/fSurveyEY[i][j][k][l]/fSurveyEY[i][j][k][l];
chi2 += dz*dz/fSurveyEZ[i][j][k][l]/fSurveyEZ[i][j][k][l];
printf("local survey %3d %3d %3d %3d %12.3f %12.3f %12.3f\n",i,j,k,l,local[0],local[1],local[2]);
printf("local ideal %12.3f %12.3f %12.3f\n",fSurveyX0[j][k][l],
fSurveyY0[j][k][l],fSurveyZ0[j][k][l]);
printf("difference %12.3f %12.3f %12.3f\n",dx,dy,dz);
}
printf("chi2 = %.2f\n",chi2);
return chi2;
}
void trdAlignmentFcn(int &npar, double *g, double &f, double *par, int iflag) {
AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit();
f = alignment->SurveyChi2(par);
if (iflag==3) {}
if (npar) {}
if (g) {}
}
void AliTRDalignment::SurveyToAlignment(int i, const char * const flag) {
if (strlen(flag)!=6) {
AliError(Form("unexpected flag: %s",flag));
return;
}
printf("Finding alignment matrix for supermodule %d\n",i);
fIbuffer[0] = i;
TFitter fitter(100);
gMinuit->SetObjectFit(this);
fitter.SetFCN(trdAlignmentFcn);
fitter.SetParameter(0,"dx",0,0.5,0,0);
fitter.SetParameter(1,"dy",0,0.5,0,0);
fitter.SetParameter(2,"dz",0,0.5,0,0);
fitter.SetParameter(3,"rx",0,0.1,0,0);
fitter.SetParameter(4,"ry",0,0.1,0,0);
fitter.SetParameter(5,"rz",0,0.1,0,0);
for (int j=0; j<6; j++) if (flag[j]=='0') fitter.FixParameter(j);
double arglist[100];
arglist[0] = 2;
fitter.ExecuteCommand("SET PRINT", arglist, 1);
fitter.ExecuteCommand("SET ERR", arglist, 1);
arglist[0]=50;
fitter.ExecuteCommand("MINIMIZE", arglist, 1);
fitter.ExecuteCommand("CALL 3", arglist,0);
double a[6];
for (int j=0; j<6; j++) a[j] = fitter.GetParameter(j);
SetSm(i,a);
for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParameter(j));
printf("\n");
for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParError(j));
printf("\n");
}
void AliTRDalignment::ReadAny(const char * const filename)
{
TString fist(filename);
if (fist.EndsWith(".txt")) ReadAscii(filename);
if (fist.EndsWith(".dat")) ReadAscii(filename);
if (fist.EndsWith(".root")) {
if (fist.Contains("Run")) ReadDB(filename);
else ReadRoot(filename);
}
}
void AliTRDalignment::WriteAscii(const char * const filename) const
{
FILE *fp = fopen(filename, "w");
if (!fp) {
AliError(Form("cannot open output file %s",filename));
return;
}
PrintSm(fp);
PrintCh(fp);
fclose(fp);
}
void AliTRDalignment::WriteRoot(const char * const filename)
{
TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
NumbersToAr(ar);
TFile fo(filename,"RECREATE");
if (fo.IsOpen()) {
fo.cd();
fo.WriteObject(ar,"TRDAlignObjs","kSingleKey");
fo.Close();
}
else AliError(Form("cannot open output file %s",filename));
delete ar;
}
void AliTRDalignment::WriteDB(const char * const filename, int run0, int run1, int ver, int subver)
{
TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
NumbersToAr(ar);
const Char_t *path = "TRD/Align/Data";
AliCDBId id(path,run0,run1);
AliCDBMetaData *md = new AliCDBMetaData();
md->SetResponsible("Dariusz Miskowiec");
md->SetComment(fComment.GetString().Data());
AliCDBEntry *e = new AliCDBEntry(ar, id, md);
e->SetVersion(ver);
e->SetSubVersion(subver);
TFile fi(filename,"RECREATE");
if (fi.IsOpen()) {
e->Write();
fi.Close();
}
else AliError(Form("cannot open input file %s",filename));
delete e;
delete md;
delete ar;
return;
}
void AliTRDalignment::WriteDB(char * const db, const char * const path, int run0, int run1)
{
TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
NumbersToAr(ar);
AliCDBManager *cdb = AliCDBManager::Instance();
AliCDBStorage *storLoc = cdb->GetStorage(db);
AliCDBMetaData *md = new AliCDBMetaData();
md->SetResponsible("Dariusz Miskowiec");
md->SetComment(fComment.GetString().Data());
AliCDBId id(path,run0,run1);
storLoc->Put(ar,id,md);
md->Delete();
delete ar;
}
void AliTRDalignment::WriteGeo(char *filename)
{
TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
NumbersToAr(ar);
delete ar;
gGeoManager->Export(filename);
}
double AliTRDalignment::GetSmRMS(int xyz) const
{
double s1 = 0.0;
double s2 = 0.0;
for (int i = 0; i < 18; i++) {
s1 += fSm[i][xyz];
s2 += fSm[i][xyz]*fSm[i][xyz];
}
double rms2 = s2/18.0 - s1*s1/18.0/18.0;
return rms2>0 ? sqrt(rms2) : 0.0;
}
double AliTRDalignment::GetChRMS(int xyz) const
{
double s1 =0.0;
double s2 =0.0;
for (int i = 0; i < 540; i++) {
s1 += fCh[i][xyz];
s2 += fCh[i][xyz]*fCh[i][xyz];
}
double rms2 = s2/540.0 - s1*s1/540.0/540.0;
return rms2>0 ? sqrt(rms2) : 0.0;
}
void AliTRDalignment::PrintSmRMS() const
{
printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f supermodule rms\n"
,GetSmRMS(0),GetSmRMS(1),GetSmRMS(2),GetSmRMS(3),GetSmRMS(4),GetSmRMS(5));
}
void AliTRDalignment::PrintChRMS() const
{
printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f chamber rms\n"
,GetChRMS(0),GetChRMS(1),GetChRMS(2),GetChRMS(3),GetChRMS(4),GetChRMS(5));
}
void AliTRDalignment::ArToNumbers(TClonesArray * const ar)
{
ar->Sort();
if (!IsGeoLoaded()) return;
for (int i = 0; i < ar->GetEntries(); i++) {
AliAlignObj *aao = (AliAlignObj *) ar->At(i);
aao->ApplyToGeometry();
}
SetZero();
ReadCurrentGeo();
}
void AliTRDalignment::NumbersToAr(TClonesArray * const ar)
{
if (!IsGeoLoaded()) return;
TClonesArray &alobj = *ar;
int nobj = 0;
for (int i = 0; i < 18; i++) {
new(alobj[nobj]) AliAlignObjParams(GetSmName(i)
,0
,fSm[i][0],fSm[i][1],fSm[i][2]
,fSm[i][3],fSm[i][4],fSm[i][5]
,0);
((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
nobj++;
}
for (int i = 0; i < 540; i++) {
if (gGeoManager->GetAlignableEntry(GetChName(i))) {
new(alobj[nobj]) AliAlignObjParams(GetChName(i)
,GetVoi(i)
,fCh[i][0],fCh[i][1],fCh[i][2]
,fCh[i][3],fCh[i][4],fCh[i][5]
,0);
((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
nobj++;
}
}
AliInfo("current geometry modified");
}
int AliTRDalignment::IsGeoLoaded()
{
if (gGeoManager) {
if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) AliWarning("current geometry is not ideal");
return 1;
} else {
AliError("first load geometry by calling TGeoManager::Import(filename)");
return 0;
}
}