From: cblume Date: Thu, 29 Mar 2007 08:12:26 +0000 (+0000) Subject: Conversion of survey data into alignable objects implemented X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=d15124a9794d1b4f2c76a1db6e70f022d30f96de Conversion of survey data into alignable objects implemented --- diff --git a/TRD/AliTRDalignment.cxx b/TRD/AliTRDalignment.cxx index 3d4786124a1..0c8c9100020 100644 --- a/TRD/AliTRDalignment.cxx +++ b/TRD/AliTRDalignment.cxx @@ -14,33 +14,33 @@ **************************************************************************/ /* $Id$ */ - -////////////////////////////////////////////////////////////////////////////////// -// // -// An AliTRDalignment object contains the alignment data (3 shifts and 3 tilts) // -// for all the alignable volumes of the TRD, i.e. for 18 supermodules and 540 // -// chambers. The class provides simple tools for reading and writing these data // -// in different formats, and for generating fake data that can be used to // -// simulate misalignment. // -// The six alignment variables have the following meaning: // -// shift in rphi // -// shift in z // -// shift in r // -// tilt around rphi // -// tilt around z // -// tilt around r // -// The shifts are in cm and the tilts are in degrees. // -// The currently supported formats are: // -// - ascii // -// - root file containing a TClonesArray of alignment objects // -// - offline conditions database // -// - OCDB-like root file // -// - geometry file (like misaligned_geometry.root) // -// // -// D.Miskowiec, November 2006 // -// // -////////////////////////////////////////////////////////////////////////////////// - +/////////////////////////////////////////////////////////////////////////////// +// // +// An AliTRDalignment object contains the alignment data (3 shifts and 3 // +// tilts) for all the alignable volumes of the TRD, i.e. for 18 supermodules // +// and 540 chambers. The class provides simple tools for reading and writing // +// these data in different formats, and for generating fake data that can be // +// used to simulate misalignment. // +// The six alignment variables have the following meaning: // +// shift in rphi // +// shift in z // +// shift in r // +// tilt around rphi // +// tilt around z // +// tilt around r // +// The shifts are in cm and the tilts are in degrees. // +// The currently supported formats are: // +// - ascii // +// - root file containing a TClonesArray of alignment objects // +// - offline conditions database // +// - OCDB-like root file // +// - geometry file (like misaligned_geometry.root) // +// // +// D.Miskowiec, November 2006 // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include #include #include @@ -49,6 +49,9 @@ #include "TGeoManager.h" #include "TGeoPhysicalNode.h" #include "TClonesArray.h" +#include "TString.h" +#include "TFitter.h" +#include "TMinuit.h" #include "AliLog.h" #include "AliAlignObj.h" @@ -61,11 +64,14 @@ #include "AliTRDalignment.h" +void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag); + ClassImp(AliTRDalignment) //_____________________________________________________________________________ AliTRDalignment::AliTRDalignment() :TObject() + ,fComment() ,fRan(0) { // @@ -74,22 +80,69 @@ AliTRDalignment::AliTRDalignment() 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; + fSurveyE[i][j][k][l] = 0.0; + } + + // Initialize the nominal positions of the survey points + // in the local frame of supermodule (where y is the long side, + // z corresponds to the radius in lab, and x to the phi in lab). + // Four survey marks are on each z-side of the supermodule. + // A B + // ----o-----------o---- x | + // \ / | + // \ / | + // \ / | + // \ / | + // ---o-----o--- --------------> + // C D y + // + // For the purpose of this explanation lets define the origin such that + // the supermodule occupies 0 < x < 77.9 cm. Then the coordinates (x,y) + // are (in cm) + // A (76.2,-30.25) + // B (76.2,+30.25) + // C ( 2.2,-22.5 ) + // D ( 2.2,+22.5 ) + // + + double x[2] = {22.5,30.25}; // lab phi, or tracking-y + double y[2] = {353.0, -353.0}; // lab z; inc. 2 cm survey target offset + double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5}; // lab r, or better tracking-x + + 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] = -pow(-1,l) * x[k]; + fSurveyY0[j][k][l] = y[j]; + fSurveyZ0[j][k][l] = z[k]; + } + } //_____________________________________________________________________________ AliTRDalignment::AliTRDalignment(const AliTRDalignment& source) :TObject(source) - ,fRan(source.fRan) + ,fComment(source.fComment) + ,fRan(source.fRan) { // // copy constructor // - for (int i = 0; i < 18; i++) { - SetSm(i,source.fSm[i]); + 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]; + fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l]; } - for (int i = 0; i < 540; i++) { - SetCh(i,source.fCh[i]); + 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]; } } @@ -104,12 +157,38 @@ 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]; + fSurveyE[i][j][k][l] = source.fSurveyE[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) +{ + // + // multiplication operator + // + + 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) { @@ -117,16 +196,8 @@ AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source) // addition operator // - 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]; - } - } + 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; @@ -139,16 +210,8 @@ AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source) // subtraction operator // - 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]; - } - } + 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; @@ -163,16 +226,8 @@ 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]); - } - } + 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; @@ -201,7 +256,7 @@ void AliTRDalignment::SetChZero() } //_____________________________________________________________________________ -void AliTRDalignment::SetSmRandom(Double_t a[6]) +void AliTRDalignment::SetSmRandom(double a[6]) { // // generate random gaussian supermodule data with sigmas a @@ -209,13 +264,11 @@ void AliTRDalignment::SetSmRandom(Double_t a[6]) double x[6]; - for (Int_t i = 0; i < 18; i++) { + for (int i = 0; i < 18; i++) { fRan.Rannor(x[0],x[1]); fRan.Rannor(x[2],x[3]); fRan.Rannor(x[4],x[5]); - for (Int_t j = 0; j < 6; j++) { - x[j] *= a[j]; - } + for (int j = 0; j < 6; j++) x[j] *= a[j]; SetSm(i,x); //PrintSm(i); } @@ -223,7 +276,7 @@ void AliTRDalignment::SetSmRandom(Double_t a[6]) } //_____________________________________________________________________________ -void AliTRDalignment::SetChRandom(Double_t a[6]) +void AliTRDalignment::SetChRandom(double a[6]) { // // generate random gaussian chamber data with sigmas a @@ -231,13 +284,11 @@ void AliTRDalignment::SetChRandom(Double_t a[6]) double x[6]; - for (Int_t i = 0; i < 540; i++) { + 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_t j = 0; j < 6; j++) { - x[j] *= a[j]; - } + for (int j = 0; j < 6; j++) x[j] *= a[j]; SetCh(i,x); //PrintCh(i); } @@ -252,7 +303,7 @@ void AliTRDalignment::SetSmFull() // expected from the mechanical precision // - Double_t a[6]; + double a[6]; a[0] = 0.3; // phi a[1] = 0.3; // z @@ -273,7 +324,7 @@ void AliTRDalignment::SetChFull() // expected from the mechanical precision // - Double_t a[6]; + double a[6]; a[0] = 0.1; // phi a[1] = 0.1; // z @@ -307,7 +358,7 @@ void AliTRDalignment::SetChResidual() // remaining after full calibration // - Double_t a[6]; + double a[6]; a[0] = 0.002; // phi a[1] = 0.003; // z @@ -321,7 +372,7 @@ void AliTRDalignment::SetChResidual() } //_____________________________________________________________________________ -void AliTRDalignment::PrintSm(Int_t i, FILE *fp) const +void AliTRDalignment::PrintSm(int i, FILE *fp) const { // // print the supermodule data @@ -334,7 +385,7 @@ void AliTRDalignment::PrintSm(Int_t i, FILE *fp) const } //_____________________________________________________________________________ -void AliTRDalignment::PrintCh(Int_t i, FILE *fp) const +void AliTRDalignment::PrintCh(int i, FILE *fp) const { // // print the chamber data @@ -360,23 +411,18 @@ void AliTRDalignment::ReadAscii(char *filename) fstream fi(filename,fstream::in); if (!fi) { - AliFatal(Form("cannot open input file %s",filename)); + AliError(Form("cannot open input file %s",filename)); + return; } // supermodules 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)); - } + 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())); - } + if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data())); SetSm(i,x); } @@ -384,16 +430,10 @@ void AliTRDalignment::ReadAscii(char *filename) 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)); - } + 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())); - } + if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data())); SetCh(i,x); } @@ -418,9 +458,7 @@ void AliTRDalignment::ReadRoot(char *filename) ArToNumbers(ar); fi.Close(); } - else { - AliError(Form("cannot open input file %s",filename)); - } + else AliError(Form("cannot open input file %s",filename)); return; @@ -438,21 +476,20 @@ void AliTRDalignment::ReadDB(char *filename) 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)); - } + else AliError(Form("cannot open input file %s",filename)); return; } //_____________________________________________________________________________ -void AliTRDalignment::ReadDB(char *db, char *path, Int_t run - , Int_t version, Int_t subversion) +void AliTRDalignment::ReadDB(char *db, char *path, int run + , int version, int subversion) { // // read the alignment data from database @@ -462,6 +499,7 @@ void AliTRDalignment::ReadDB(char *db, char *path, Int_t run AliCDBStorage *storLoc = cdb->GetStorage(db); AliCDBEntry *e = storLoc->Get(path,run,version,subversion); e->PrintMetaData(); + fComment.SetString(e->GetMetaData()->GetComment()); TClonesArray *ar = (TClonesArray *) e->GetObject(); ArToNumbers(ar); @@ -486,29 +524,17 @@ void AliTRDalignment::ReadGeo(char *misaligned) TGeoManager::Import(misaligned); for (int i = 0; i < 18; i++) { TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetSmName(i)); - if (!pne) { - AliError(Form("no such physical node entry: %s",GetSmName(i))); - return; - } + if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i))); TGeoPhysicalNode *node = pne->GetPhysicalNode(); - if (!node) { - AliError(Form("physical node entry %s has no physical node",GetSmName(i))); - return; - } + if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i))); misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix()); ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix()); } for (int i = 0; i < 540; i++) { TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetChName(i)); - if (!pne) { - AliError(Form("no such physical node entry: %s",GetChName(i))); - return; - } + if (!pne) AliError(Form("no such physical node entry: %s",GetChName(i))); TGeoPhysicalNode *node = pne->GetPhysicalNode(); - if (!node) { - AliError(Form("physical node entry %s has no physical node",GetChName(i))); - return; - } + if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i))); misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix()); ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix()); } @@ -524,10 +550,7 @@ void AliTRDalignment::ReadGeo(char *misaligned) 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; - } + 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]); @@ -568,14 +591,230 @@ void AliTRDalignment::ReadGeo(char *misaligned) //_____________________________________________________________________________ void AliTRDalignment::ReadSurveyReport(char *filename) { - // read survey report and set the supermodule parameters correspondingly + // + // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ, + // and fSurveyE. Store the survey info in the fComment. + // Each supermodule has 8 survey points. The point names look like + // TRD_sm08ah0 and have the following meaning. + // + // sm00..17 mean supermodule 0 through 17, following the phi. + // Supermodule 00 is between phi=0 and phi=20 degrees. + // + // a or c denotes the anticlockwise and clockwise end of the supermodule + // in z. Clockwise end is where z is negative and where the muon arm sits. + // + // l or h denote low radius and high radius holes + // + // 0 or 1 denote the hole at smaller and at larger phi, respectively. + // - fstream fi(filename,fstream::in); - if (!fi) { - AliFatal(Form("cannot open input file %s",filename)); + // read the survey file + + fstream in(filename,fstream::in); + if (!in) { + AliError(Form("cannot open input file %s",filename)); + return; + } + + // loop through the lines of the file until the beginning of data + + 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; + } + + // check what we found so far (watch out, they have \r at the end) + + std::cout<<"title .........."<> pna >> x >> y >> z >> type >> precision; + if (in.fail()) break; + if (pna(0,6)!="TRD_sm") { + AliError(Form("unexpected point name: %s",pna.Data())); + break; + } + int i = atoi(pna(6,0).Data()); // supermodule number + int j = -1; + if (pna(8) == 'a') j=0; // anticlockwise, positive z + if (pna(8) == 'c') j=1; // clockwise, negative z + int k = -1; + if (pna(9) == 'l') k=0; // low radius + if (pna(9) == 'h') k=1; // high radius + int l = atoi(pna(10,0).Data()); // phi within supermodule + if (i>=0 && i<18 && j>=0 && j<2 && k>=0 && k<2 && l>=0 && l<2) { + fSurveyX[i][j][k][l] = tocm*x; + fSurveyY[i][j][k][l] = tocm*y; + fSurveyZ[i][j][k][l] = tocm*z; + fSurveyE[i][j][k][l] = precision/10; // "precision" is supposed to be in mm + std::cout << "decoded "<GetAlignableEntry(GetSmName(i)); + if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i))); + TGeoPhysicalNode *node = pne->GetPhysicalNode(); + if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i))); + + // al.ApplyToGeometry(); + // node = pne->GetPhysicalNode(); // changed in the meantime + // TGeoHMatrix *ma = node->GetMatrix(); + + // a less destructive method (it does not modify geometry), gives the same result: + + TGeoHMatrix *ma = new TGeoHMatrix(); + al.GetLocalMatrix(*ma); + ma->MultiplyLeft(node->GetMatrix()); // global trafo, modified by a[] + + 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 (fSurveyE[i][j][k][l] == 0.0) continue; // no data for this survey point + 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+dy*dy+dz*dz)/fSurveyE[i][j][k][l]/fSurveyE[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 Fcn(int &npar, double *g, double &f, double *par, int iflag) { + + // + // Standard function as needed by Minuit-like minimization procedures. + // For the set of parameters par calculates and returns chi-squared. + // + + // smuggle a C++ object into a C function + AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit(); + + f = alignment->SurveyChi2(par); + if (iflag==3); + if (npar); + if (g); // no warnings about unused stuff... + +} + +//_____________________________________________________________________________ +void AliTRDalignment::SurveyToAlignment(int i,char *flag) { + + // + // Find the supermodule alignment parameters needed to make the survey + // results coincide with the ideal positions of the survey marks. + // The string flag should look like "101000"; the six characters corresponds + // to the six alignment parameters and 0/1 mean that the parameter should + // be fixed/released in the fit. + + if (strlen(flag)!=6) { + AliError(Form("unexpected flag: %s",flag)); + return; + } + + printf("Finding alignment matrix for supermodule %d\n",i); + fIbuffer[0] = i; // store the sm number in the buffer so minuit can see it + + TFitter fitter(100); + gMinuit->SetObjectFit(this); + fitter.SetFCN(Fcn); + 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("SIMPLEX", arglist, 1); + 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"); } @@ -587,19 +826,11 @@ void AliTRDalignment::ReadAny(char *filename) // TString fist(filename); - if (fist.EndsWith(".txt")) { - ReadAscii(filename); - } - if (fist.EndsWith(".dat")) { - ReadAscii(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); - } + if (fist.Contains("Run")) ReadDB(filename); + else ReadRoot(filename); } } @@ -639,16 +870,14 @@ void AliTRDalignment::WriteRoot(char *filename) fo.WriteObject(ar,"TRDAlignObjs","kSingleKey"); fo.Close(); } - else { - AliError(Form("cannot open output file %s",filename)); - } + else AliError(Form("cannot open output file %s",filename)); delete ar; } //_____________________________________________________________________________ -void AliTRDalignment::WriteDB(char *filename, char *comment, Int_t run0, Int_t run1) +void AliTRDalignment::WriteDB(char *filename, int run0, int run1) { // // dumping on a DB-like file @@ -660,16 +889,14 @@ void AliTRDalignment::WriteDB(char *filename, char *comment, Int_t run0, Int_t r AliCDBId id(path,run0,run1); AliCDBMetaData *md = new AliCDBMetaData(); md->SetResponsible("Dariusz Miskowiec"); - md->SetComment(comment); + md->SetComment(fComment.GetString().Data()); AliCDBEntry *e = new AliCDBEntry(ar, id, md); TFile fi(filename,"RECREATE"); if (fi.IsOpen()) { e->Write(); fi.Close(); } - else { - AliError(Form("cannot open input file %s",filename)); - } + else AliError(Form("cannot open input file %s",filename)); delete e; delete md; @@ -680,7 +907,7 @@ void AliTRDalignment::WriteDB(char *filename, char *comment, Int_t run0, Int_t r } //_____________________________________________________________________________ -void AliTRDalignment::WriteDB(char *db, char *path, char *comment, Int_t run0, Int_t run1) +void AliTRDalignment::WriteDB(char *db, char *path, int run0, int run1) { // // store the alignment data in database @@ -692,7 +919,7 @@ void AliTRDalignment::WriteDB(char *db, char *path, char *comment, Int_t run0, I AliCDBStorage *storLoc = cdb->GetStorage(db); AliCDBMetaData *md = new AliCDBMetaData(); md->SetResponsible("Dariusz Miskowiec"); - md->SetComment(comment); + md->SetComment(fComment.GetString().Data()); AliCDBId id(path,run0,run1); storLoc->Put(ar,id,md); md->Delete(); @@ -720,38 +947,38 @@ void AliTRDalignment::WriteGeo(char *filename) } //_____________________________________________________________________________ -Double_t AliTRDalignment::GetSmRMS(Int_t xyz) const +double AliTRDalignment::GetSmRMS(int xyz) const { // // rms fSm[][xyz] // - Double_t s1 = 0.0; - Double_t s2 = 0.0; + 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_t rms2 = s2/18.0 - s1*s1/18.0/18.0; + double rms2 = s2/18.0 - s1*s1/18.0/18.0; return rms2>0 ? sqrt(rms2) : 0.0; } //_____________________________________________________________________________ -Double_t AliTRDalignment::GetChRMS(Int_t xyz) const +double AliTRDalignment::GetChRMS(int xyz) const { // // rms fCh[][xyz] // - Double_t s1 =0.0; - Double_t s2 =0.0; + 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_t rms2 = s2/540.0 - s1*s1/540.0/540.0; + double rms2 = s2/540.0 - s1*s1/540.0/540.0; return rms2>0 ? sqrt(rms2) : 0.0; @@ -847,16 +1074,10 @@ void AliTRDalignment::LoadIdealGeometry(char *filename) static int attempt = 0; // which reload attempt is it? just to avoid endless loops - if (!gGeoManager) { - TGeoManager::Import(filename); - } - if (!gGeoManager) { - AliFatal(Form("cannot open geometry file %s",filename)); - } + if (!gGeoManager) TGeoManager::Import(filename); + if (!gGeoManager) AliError(Form("cannot open geometry file %s",filename)); if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) { - if (attempt) { - AliFatal(Form("geometry on file %s is not ideal",filename)); - } + if (attempt) AliError(Form("geometry on file %s is not ideal",filename)); AliWarning("current geometry is not ideal - it contains physical nodes"); AliWarning(Form("reloading geometry from %s - attempt nr %d",filename,attempt)); gGeoManager = 0; diff --git a/TRD/AliTRDalignment.h b/TRD/AliTRDalignment.h index 3698734c626..4bdcd6f136c 100644 --- a/TRD/AliTRDalignment.h +++ b/TRD/AliTRDalignment.h @@ -15,6 +15,7 @@ #include #include +#include #include @@ -25,42 +26,44 @@ class AliTRDalignment : public TObject { AliTRDalignment(); AliTRDalignment(const AliTRDalignment& source); AliTRDalignment& operator=(const AliTRDalignment& source); + AliTRDalignment& operator*=(double fac); AliTRDalignment& operator+=(const AliTRDalignment& source); AliTRDalignment& operator-=(const AliTRDalignment& source); Bool_t operator==(const AliTRDalignment& source) const; - virtual ~AliTRDalignment() {}; // destructor + virtual ~AliTRDalignment() {}; // setting void SetSmZero(); // reset to zero supermodule data void SetChZero(); // reset to zero chamber data void SetZero() {SetSmZero(); SetChZero();} // reset to zero both - void SetSm(Int_t sm, const Double_t a[6]) { for (int i = 0; i < 6; i++) fSm[sm][i] = a[i]; } - void SetCh(Int_t ch, const Double_t a[6]) { for (int i = 0; i < 6; i++) fCh[ch][i] = a[i]; } - void SetSmRandom(Double_t a[6]); // generate random gaussians with sigmas a - void SetChRandom(Double_t a[6]); // generate random gaussians with sigmas a + void SetSm(int sm, const double a[6]) {for (int i = 0; i < 6; i++) fSm[sm][i] = a[i];} + void SetCh(int ch, const double a[6]) {for (int i = 0; i < 6; i++) fCh[ch][i] = a[i];} + void SetSmRandom(double a[6]); // generate random gaussians with sigmas a + void SetChRandom(double a[6]); // generate random gaussians with sigmas a void SetSmFull(); // set supermodule data to initial aka full void SetChFull(); // set chamber data to initial aka full void SetSmResidual(); // set supermodule data to final aka residual void SetChResidual(); // set chamber data to final aka residual - void SetFull() {SetSmFull(); SetChFull();} - void SetResidual() {SetSmResidual(); SetChResidual();} + void SetFull() {SetSmFull(); SetChFull();} + void SetResidual() {SetSmResidual(); SetChResidual();} + void SetComment(char *s) {fComment.SetString(s);} // dumping on screen - void PrintSm(Int_t sm, FILE *fp = stdout) const; // print data of a supermodule - void PrintCh(Int_t sm, FILE *fp = stdout) const; // print data of a chamber - void PrintSm(FILE *fp = stdout) const { for (int i = 0; i < 18; i++) PrintSm(i,fp); } - void PrintCh(FILE *fp = stdout) const { for (int i = 0; i < 540; i++) PrintCh(i,fp); } - void Print(FILE *fp = stdout) const { PrintSm(fp); PrintCh(fp); } - void Print(Option_t *) const { Print(); } + void PrintSm(int sm, FILE *fp = stdout) const; // print data of a supermodule + void PrintCh(int sm, FILE *fp = stdout) const; // print data of a chamber + void PrintSm(FILE *fp = stdout) const {for (int i = 0; i < 18; i++) PrintSm(i,fp);} + void PrintCh(FILE *fp = stdout) const {for (int i = 0; i < 540; i++) PrintCh(i,fp);} + void Print(FILE *fp = stdout) const {PrintSm(fp); PrintCh(fp); } + void Print(Option_t *) const {Print(); } // reading-in from file void ReadAscii(char *filename); // read from ascii file void ReadRoot(char *filename); // read from root file void ReadDB(char *filename); // read from DB file - void ReadDB(char *db, char *path, Int_t run, Int_t version=-1, Int_t subversion=-1); + void ReadDB(char *db, char *path, int run, int version=-1, int subversion=-1); void ReadGeo(char *misaligned); // read from misaligned_geometry.root void ReadSurveyReport(char *filename); // read from survey report void ReadAny(char *filename); // read from any kind of file @@ -69,49 +72,71 @@ class AliTRDalignment : public TObject { void WriteAscii(char *filename) const; // store data on ascii file void WriteRoot(char *filename); // store data on root file - void WriteDB(char *filename, char *comment, Int_t run0, Int_t run1); - // store data on a local DB-like file - void WriteDB(char *db, char *path, char *comment, Int_t run0, Int_t run1); + void WriteDB(char *filename, int run0, int run1); // store data on a local DB-like file + void WriteDB(char *db, char *path, int run0, int run1); // store data on DB file void WriteGeo(char *filename); // apply misalignment and store geometry // geometry and symbolic names getters // phi-sector number of chamber ch, 0-17 - Int_t GetSec(Int_t ch) const { return ch/30; } + int GetSec(int ch) const {return ch/30;} // stack number, 0-4 - Int_t GetSta(Int_t ch) const { return ch%30/6; } - // plane number, 0-5 - Int_t GetPla(Int_t ch) const { return ch%30%6; } + int GetSta(int ch) const {return ch%30/6;} + // plane number, 0-5 + int GetPla(int ch) const {return ch%30%6;} // module number, 0-89 - Int_t GetMod(Int_t ch) const { return 5*GetSec(ch)+GetSta(ch); } + int GetMod(int ch) const {return 5*GetSec(ch)+GetSta(ch); } // layer number, 9-14 - Int_t GetLay(Int_t ch) const { return AliAlignObj::kTRD1+GetPla(ch); } + int GetLay(int ch) const {return AliAlignObj::kTRD1+GetPla(ch); } // volume id - UShort_t GetVoi(Int_t ch) const { return AliAlignObj::LayerToVolUID(GetLay(ch),GetMod(ch)); } - char *GetSmName(Int_t sm) const { return Form("TRD/sm%02d",sm); } - char *GetChName(Int_t ch) const { return Form("TRD/sm%02d/st%d/pl%d",GetSec(ch),GetSta(ch),GetPla(ch)); } + UShort_t GetVoi(int ch) const {return AliAlignObj::LayerToVolUID(GetLay(ch),GetMod(ch)); } + // symbolic name of a supermodule + char *GetSmName(int sm) const {return Form("TRD/sm%02d",sm); } + // symbolic name of a chamber + char *GetChName(int ch) const {return Form("TRD/sm%02d/st%d/pl%d",GetSec(ch),GetSta(ch),GetPla(ch)); } // data analysis - Double_t GetSmRMS(Int_t xyz) const; // calculate rms fSm[*][xyz] - Double_t GetChRMS(Int_t xyz) const; // calculate rms fCh[*][xyz] - void PrintSmRMS() const; // print rms of fSm - void PrintChRMS() const; // print rms of fCh - void PrintRMS() const { PrintSmRMS(); PrintChRMS();} + double GetSmRMS(int xyz) const; // calculate rms fSm[*][xyz] + double GetChRMS(int xyz) const; // calculate rms fCh[*][xyz] + void PrintSmRMS() const; // print rms of fSm + void PrintChRMS() const; // print rms of fCh + void PrintRMS() const {PrintSmRMS(); PrintChRMS();} + + double SurveyChi2(int i, double *a); // compare survey with ideal, return chi2 + double SurveyChi2(double *a) {return SurveyChi2(fIbuffer[0],a);} + void SurveyToAlignment(int i, char *flag); // determine alignment of supermodule i based on survey protected: - void ArToNumbers(TClonesArray *ar); // read ar and fill fSm and fCh - void NumbersToAr(TClonesArray *ar); // build ar using fSm and fCh data - void LoadIdealGeometry(char *filename); // load ideal geometry from file - void LoadIdealGeometry() {LoadIdealGeometry("geometry.root");} + void ArToNumbers(TClonesArray *ar); // read ar and fill fSm and fCh + void NumbersToAr(TClonesArray *ar); // build ar using fSm and fCh data + void LoadIdealGeometry(char *filename); // load ideal geometry from file + void LoadIdealGeometry() {LoadIdealGeometry("geometry.root");} protected: - Double_t fSm[18][6]; // supermodule data - Double_t fCh[540][6]; // chamber data - TRandom fRan; // random generator for fake alignment data + double fSm[18][6]; // supermodule data + double fCh[540][6]; // chamber data + TObjString fComment; // info concerning origin of the data etc. + TRandom fRan; // random generator for fake alignment data + + // Temporary storage for ideal position of the survey points and the survey data. + // The survey data are in master frame and in cm. Each supermodule has 8 survey marks. + // The indices are sm number, z-end, radius, phi. + // The ideal positions of survey points are in local frame of supermodule and in cm. + // The indices are z-end, radius, phi. + // The processed survey results are stored in fSm. + double fSurveyX[18][2][2][2]; // supermodule survey point X + double fSurveyY[18][2][2][2]; // supermodule survey point Y + double fSurveyZ[18][2][2][2]; // supermodule survey point Z + double fSurveyE[18][2][2][2]; // supermodule survey point error + double fSurveyX0[2][2][2]; // ideal X position of the survey marks + double fSurveyY0[2][2][2]; // ideal Y position of the survey marks + double fSurveyZ0[2][2][2]; // ideal Z position of the survey marks + int fIbuffer[1000]; // generic buffer for misc. operations + double fDbuffer[1000]; // generic buffer for misc. operations ClassDef(AliTRDalignment,1)