Conversion of survey data into alignable objects implemented
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Mar 2007 08:12:26 +0000 (08:12 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Mar 2007 08:12:26 +0000 (08:12 +0000)
TRD/AliTRDalignment.cxx
TRD/AliTRDalignment.h

index 3d47861..0c8c910 100644 (file)
  **************************************************************************/
 
 /* $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 <iostream>
 #include <fstream>
 #include <string>
 
@@ -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"
 
 #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 .........."<<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; // we want to have it in cm
+  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;
   }
 
-  // to be continued...
+  // scan the rest of the file which should contain list of surveyed points
+  // for every point, decode the point name and store the numbers in the right 
+  // place in the arrays fSurveyX etc.
+
+  while (1) {
+    TString pna; // point name
+    char type;
+    double x,y,z,precision;
+    in >> 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 "<<pna<<" "
+               <<fSurveyX[i][j][k][l]<<" "
+               <<fSurveyY[i][j][k][l]<<" "
+               <<fSurveyZ[i][j][k][l]<<" "
+               <<fSurveyE[i][j][k][l]<<std::endl;      
+    } 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());
+                        
+}
+
+//_____________________________________________________________________________
+double AliTRDalignment::SurveyChi2(int i, double *a) {
+
+  //
+  // Compare the survey results to the ideal positions of the survey marks
+  // in the local frame of supermodule. When transforming, use the alignment 
+  // parameters a[6]. Return chi-squared.
+  //
+
+  if (!gGeoManager) LoadIdealGeometry();
+  printf("Survey of supermodule %d\n",i);
+  AliAlignObjAngles 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) 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;
index 3698734..4bdcd6f 100644 (file)
@@ -15,6 +15,7 @@
 
 #include <TObject.h>
 #include <TRandom.h>
+#include <TObjString.h>
 
 #include <AliAlignObj.h>
 
@@ -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)