]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDalignment.cxx
temporary bug fix (speed issues)
[u/mrichter/AliRoot.git] / TRD / AliTRDalignment.cxx
index 797ed2724fb8b92e82f345a5c7362daece6e38aa..e41c6bec3a804580280f9fec681172fe4a70340e 100644 (file)
@@ -46,7 +46,7 @@
 // TGeoManager::Import("geometry.root"); a.WriteRoot("kuku.root");           //
 // TGeoManager::Import("geometry.root"); a.WriteDB("kukudb.root",0,0);       //
 // TGeoManager::Import("geometry.root");                                     //
-// a.WriteDB("local://$ALICE_ROOT", "TRD/Align/Data", 0,0);                  //
+// a.WriteDB("local://$ALICE_ROOT/OCDB", "TRD/Align/Data", 0,0);                  //
 // TGeoManager::Import("geometry.root"); a.WriteGeo("kukugeometry.root");    //
 //                                                                           //
 // b.ReadAscii("kuku.dat");                                                  //
@@ -67,7 +67,7 @@
 
 #include <iostream>
 #include <fstream>
-#include <string>
+//#include <string>
 
 #include "TMath.h"
 #include "TFile.h"
 
 #include "AliLog.h"
 #include "AliAlignObj.h"
-#include "AliAlignObjAngles.h"
+#include "AliAlignObjParams.h"
 #include "AliCDBManager.h"
 #include "AliCDBStorage.h"
 #include "AliCDBMetaData.h"
 #include "AliCDBEntry.h"
 #include "AliCDBId.h"
+#include "AliSurveyObj.h"
+#include "AliSurveyPoint.h"
 
 #include "AliTRDalignment.h"
 
@@ -109,7 +111,9 @@ AliTRDalignment::AliTRDalignment()
     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;
+    fSurveyEX[i][j][k][l] = 0.0;
+    fSurveyEY[i][j][k][l] = 0.0;
+    fSurveyEZ[i][j][k][l] = 0.0;
   }
 
   // Initialize the nominal positions of the survey points 
@@ -162,7 +166,9 @@ AliTRDalignment::AliTRDalignment(const AliTRDalignment& source)
     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];
+    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];
@@ -186,7 +192,9 @@ AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source)
       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];
+      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];
@@ -288,12 +296,12 @@ 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++) {
-    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];
+    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);
     //PrintSm(i);
   }
@@ -477,15 +485,13 @@ void AliTRDalignment::ReadCurrentGeo()
   // represent just the last alignment. -- check this!
   //
 
+  TGeoPNEntry *pne;
   TGeoHMatrix *ideSm[18];  // ideal
   TGeoHMatrix *misSm[18];  // misaligned
-  for (int i = 0; i < 18; i++) {
+  for (int i = 0; i < 18; i++) if ((pne = gGeoManager->GetAlignableEntry(GetSmName(i)))) {
 
     // read misaligned and original matrices
 
-    TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetSmName(i));
-    if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
-    if (!pne) continue;
     TGeoPhysicalNode *node = pne->GetPhysicalNode();
     if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
     if (!node) continue;
@@ -517,13 +523,10 @@ void AliTRDalignment::ReadCurrentGeo()
 
   TGeoHMatrix *ideCh[540]; // ideal
   TGeoHMatrix *misCh[540]; // misaligned
-  for (int i = 0; i < 540; i++) {
+  for (int i = 0; i < 540; i++) if ((pne = gGeoManager->GetAlignableEntry(GetChName(i)))) {
 
     // read misaligned and original matrices
 
-    TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetChName(i));
-    if (!pne) AliError(Form("no such physical node entry: %s",GetChName(i)));
-    if (!pne) continue;
     TGeoPhysicalNode *node = pne->GetPhysicalNode();
     if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i)));
     if (!node) continue;
@@ -621,6 +624,30 @@ void AliTRDalignment::ReadDB(char *db, char *path, int run
   }
 }
 
+//_____________________________________________________________________________
+Bool_t AliTRDalignment::DecodeSurveyPointName(TString pna, Int_t &sm, Int_t &iz, 
+                                             Int_t &ir, Int_t &iphi) {
+  // decode the survey point name and extract the sm, z, r and phi indices
+  
+  if (pna(0,6)!="TRD_sm") {
+    AliError(Form("unexpected point name: %s",pna.Data()));
+    return kFALSE;
+  }
+  sm = atoi(pna(6,2).Data()); // supermodule number
+  iz = -1;
+  if (pna(8) == 'a') iz=0; // anticlockwise, positive z
+  if (pna(8) == 'c') iz=1; // clockwise, negative z
+  ir = -1;
+  if (pna(9) == 'l') ir=0; // low radius
+  if (pna(9) == 'h') ir=1; // high radius
+  iphi = -1;
+  if (pna(10) == '0') iphi = 0; // low phi within supermodule
+  if (pna(10) == '1') iphi = 1; // high phi within supermodule
+  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(char *filename) 
 {
@@ -703,39 +730,123 @@ void AliTRDalignment::ReadSurveyReport(char *filename)
 
   while (1) {
     TString pna; // point name
-    char type;
+    char type, target;
     double x,y,z,precision;
-    in >> pna >> x >> y >> z >> type >> precision;  
+    
+    in >> pna >> x >> y >> z >> type >> target >> 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) {
+    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;
-      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;      
+      fSurveyEX[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
+      fSurveyEY[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
+      fSurveyEZ[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
+      // if, at some point, separate precision numbers for x,y,z show up in the 
+      // survey reports the function will fail here
+      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(AliSurveyObj *so) 
+{
+  //
+  // 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.
+  //
+
+  // read and process the data from the survey object
+
+  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();
+
+  // 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;
+  }
+
+  // for every survey point, decode the point name and store the numbers in 
+  // the right place in the arrays fSurveyX etc.
+
+  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; // "precision" is supposed to be in mm
+      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());                      
 }
 
 //_____________________________________________________________________________
@@ -749,7 +860,7 @@ double AliTRDalignment::SurveyChi2(int i, double *a) {
 
   if (!IsGeoLoaded()) return 0;
   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);
+  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();
@@ -768,14 +879,18 @@ double AliTRDalignment::SurveyChi2(int i, double *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
+    if (fSurveyEX[i][j][k][l] == 0.0 
+       && fSurveyEY[i][j][k][l] == 0.0 
+       && fSurveyEZ[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];
+    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]);
@@ -797,9 +912,9 @@ 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); // no warnings about unused stuff...
+  if (iflag==3) {}
+  if (npar) {} 
+  if (g) {} // no warnings about unused stuff...
 
 }
 
@@ -895,7 +1010,7 @@ void AliTRDalignment::WriteRoot(char *filename)
   // store the alignment data on root file
   //
 
-  TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
+  TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
   NumbersToAr(ar);
   TFile fo(filename,"RECREATE");
   if (fo.IsOpen()) {
@@ -916,9 +1031,9 @@ void AliTRDalignment::WriteDB(char *filename, int run0, int run1)
   // dumping on a DB-like file
   //
 
-  TClonesArray   *ar = new TClonesArray("AliAlignObjAngles",10000);
+  TClonesArray   *ar = new TClonesArray("AliAlignObjParams",10000);
   NumbersToAr(ar);
-  char *path = "TRD/Align/Data";
+  const Char_t *path = "TRD/Align/Data";
   AliCDBId id(path,run0,run1);
   AliCDBMetaData *md = new AliCDBMetaData();
   md->SetResponsible("Dariusz Miskowiec");
@@ -946,7 +1061,7 @@ void AliTRDalignment::WriteDB(char *db, char *path, int run0, int run1)
   // store the alignment data in database
   //
 
-  TClonesArray   *ar      = new TClonesArray("AliAlignObjAngles",10000);
+  TClonesArray   *ar      = new TClonesArray("AliAlignObjParams",10000);
   NumbersToAr(ar);
   AliCDBManager  *cdb     = AliCDBManager::Instance();
   AliCDBStorage  *storLoc = cdb->GetStorage(db);
@@ -968,7 +1083,7 @@ void AliTRDalignment::WriteGeo(char *filename)
   // resulting geometry on a root file
   //
 
-  TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
+  TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
   NumbersToAr(ar);
   delete ar;
   gGeoManager->Export(filename);
@@ -1071,7 +1186,7 @@ void AliTRDalignment::NumbersToAr(TClonesArray *ar)
   TClonesArray &alobj = *ar;
   int nobj = 0;
   for (int i = 0; i <  18; i++) {      
-      new(alobj[nobj]) AliAlignObjAngles(GetSmName(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]
@@ -1081,13 +1196,15 @@ void AliTRDalignment::NumbersToAr(TClonesArray *ar)
   }
 
   for (int i = 0; i < 540; i++) {
-    new(alobj[nobj]) AliAlignObjAngles(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++;
+    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");