TRD/AliTRDalignment.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // An AliTRDalignment object contains the alignment data (3 shifts and 3     //
  20 // tilts) for all the alignable volumes of the TRD, i.e. for 18 supermodules //
  21 // and 540 chambers. The class provides simple tools for reading and writing //
  22 // these data in different formats, and for generating fake data that can be //
  23 // used to simulate misalignment.                                            //
  24 // The six alignment variables have the following meaning:                   //
  25 // shift in rphi                                                             //
  26 // shift in z                                                                //
  27 // shift in r                                                                //
  28 // tilt around rphi                                                          //
  29 // tilt around z                                                             //
  30 // tilt around r                                                             //
  31 // The shifts are in cm and the tilts are in degrees.                        //
  32 // The currently supported formats are:                                      //
  33 // - ascii                                                                   //
  34 // - root file containing a TClonesArray of alignment objects                //
  35 // - offline conditions database                                             //
  36 // - OCDB-like root file                                                     //
  37 // - geometry file (like misaligned_geometry.root)                           //
  38 //                                                                           //
  39 // Some examples of usage (in an aliroot session):                           //
  40 // AliTRDalignment a,b,c,d,e;                                                //
  41 // double xsm[]={0,0,0,-70,0,0};                                             //
  42 // double xch[]={0,0,-50,0,0,0};                                             //
  43 // a.SetSm(4,xsm);                                                           //
  44 // a.SetCh(120,xch);                                                         //
  45 // a.WriteAscii("kuku.dat");                                                 //
  46 // TGeoManager::Import("geometry.root"); a.WriteRoot("kuku.root");           //
  47 // TGeoManager::Import("geometry.root"); a.WriteDB("kukudb.root",0,0);       //
  48 // TGeoManager::Import("geometry.root");                                     //
  49 // a.WriteDB("local://$ALICE_ROOT", "TRD/Align/Data", 0,0);                  //
  50 // TGeoManager::Import("geometry.root"); a.WriteGeo("kukugeometry.root");    //
  51 //                                                                           //
  52 // b.ReadAscii("kuku.dat");                                                  //
  53 // TGeoManager::Import("geometry.root"); c.ReadRoot("kuku.root");            //
  54 // TGeoManager::Import("geometry.root"); d.ReadDB("kukudb.root");            //
  55 // TGeoManager::Import("kukugeometry.root"); e.ReadCurrentGeo();             //
  56 //                                                                           //
  57 // e.PrintSm(4);                                                             //
  58 // e.PrintCh(120);                                                           //
  59 // a.PrintRMS();                                                             //
  60 // b.PrintRMS();                                                             //
  61 // e.PrintRMS();                                                             //
  62 //                                                                           //
  63 //                                                                           //
  64 // D.Miskowiec, November 2006                                                //
  65 //                                                                           //
  66 ///////////////////////////////////////////////////////////////////////////////
  67
  68 #include <iostream>
  69 #include <fstream>
  70 #include <string>
  71
  72 #include "TMath.h"
  73 #include "TFile.h"
  74 #include "TGeoManager.h"
  75 #include "TGeoPhysicalNode.h"
  76 #include "TClonesArray.h"
  77 #include "TString.h"
  78 #include "TFitter.h"
  79 #include "TMinuit.h"
  80
  81 #include "AliLog.h"
  82 #include "AliAlignObj.h"
  83 #include "AliAlignObjAngles.h"
  84 #include "AliCDBManager.h"
  85 #include "AliCDBStorage.h"
  86 #include "AliCDBMetaData.h"
  87 #include "AliCDBEntry.h"
  88 #include "AliCDBId.h"
  89
  90 #include "AliTRDalignment.h"
  91
  92 void trdAlignmentFcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
  93
  94 ClassImp(AliTRDalignment)
  95
  96 //_____________________________________________________________________________
  97 AliTRDalignment::AliTRDalignment()
  98   :TObject()
  99   ,fComment()
 100   ,fRan(0)
 101 {
 102   //
 103   // constructor
 104   //
 105
 106   SetZero();
 107
 108   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++) {
 109     fSurveyX[i][j][k][l] = 0.0;
 110     fSurveyY[i][j][k][l] = 0.0;
 111     fSurveyZ[i][j][k][l] = 0.0;
 112     fSurveyE[i][j][k][l] = 0.0;
 113   }
 114
 115   // Initialize the nominal positions of the survey points
 116   // in the local frame of supermodule (where y is the long side,
 117   // z corresponds to the radius in lab, and x to the phi in lab).
 118   // Four survey marks are on each z-side of the supermodule.
 119   //               A           B
 120   //           ----o-----------o----        x |
 121   //           \                   /          |
 122   //            \                 /           |
 123   //             \               /            |
 124   //              \             /             |
 125   //               ---o-----o---              -------------->
 126   //                  C     D                              y
 127   //
 128   // For the purpose of this explanation lets define the origin such that
 129   // the supermodule occupies 0 < x < 77.9 cm. Then the coordinates (x,y)
 130   // are (in cm)
 131   // A (76.2,-30.25)
 132   // B (76.2,+30.25)
 133   // C ( 2.2,-22.5 )
 134   // D ( 2.2,+22.5 )
 135   //
 136
 137   double x[2] = {22.5,30.25};                   // lab phi, or tracking-y
 138   double y[2] = {353.0, -353.0};                // lab z; inc. 2 cm survey target offset
 139   double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5}; // lab r, or better tracking-x
 140
 141   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
 142     fSurveyX0[j][k][l] = -TMath::Power(-1,l) * x[k];
 143     fSurveyY0[j][k][l] = y[j];
 144     fSurveyZ0[j][k][l] = z[k];
 145   }
 146
 147 }
 148
 149 //_____________________________________________________________________________
 150 AliTRDalignment::AliTRDalignment(const AliTRDalignment& source)
 151   :TObject(source)
 152   ,fComment(source.fComment)
 153   ,fRan(source.fRan)
 154 {
 155   //
 156   // copy constructor
 157   //
 158
 159   for (int i=0; i<18; i++) SetSm(i,source.fSm[i]);
 160   for (int i=0; i<540; i++) SetCh(i,source.fCh[i]);
 161   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++) {
 162     fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
 163     fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
 164     fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
 165     fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
 166   }
 167   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
 168     fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
 169     fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
 170     fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
 171   }
 172
 173 }
 174
 175 //_____________________________________________________________________________
 176 AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source)
 177 {
 178   //
 179   // assignment operator
 180   //
 181
 182   if (this != &source) {
 183     for (int i = 0; i <  18; i++) SetSm(i,source.fSm[i]);
 184     for (int i = 0; i < 540; i++) SetCh(i,source.fCh[i]);
 185     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++) {
 186       fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
 187       fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
 188       fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
 189       fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
 190     }
 191     for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
 192       fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
 193       fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
 194       fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
 195     }
 196     fComment = source.fComment;
 197   }
 198
 199   return *this;
 200
 201 }
 202
 203 //_____________________________________________________________________________
 204 AliTRDalignment& AliTRDalignment::operator*=(double fac)
 205 {
 206   //
 207   // multiplication operator
 208   //
 209
 210   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] *= fac;
 211   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] *= fac;
 212
 213   return *this;
 214
 215 }
 216
 217 //_____________________________________________________________________________
 218 AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source)
 219 {
 220   //
 221   // addition operator
 222   //
 223
 224   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] += source.fSm[i][j];
 225   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] += source.fCh[i][j];
 226
 227   return *this;
 228
 229 }
 230
 231 //_____________________________________________________________________________
 232 AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source)
 233 {
 234   //
 235   // subtraction operator
 236   //
 237
 238   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) fSm[i][j] -= source.fSm[i][j];
 239   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) fCh[i][j] -= source.fCh[i][j];
 240
 241   return *this;
 242
 243 }
 244
 245 //_____________________________________________________________________________
 246 Bool_t AliTRDalignment::operator==(const AliTRDalignment &source) const
 247 {
 248   //
 249   // comparison operator
 250   //
 251
 252   Bool_t areEqual = 1;
 253
 254   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) areEqual &= (fSm[i][j] == source.fSm[i][j]);
 255   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) areEqual &= (fCh[i][j] == source.fCh[i][j]);
 256
 257   return areEqual;
 258
 259 }
 260
 261 //_____________________________________________________________________________
 262 void AliTRDalignment::SetSmZero()
 263 {
 264   //
 265   // reset to zero supermodule data
 266   //
 267
 268   memset(&fSm[0][0],0,sizeof(fSm));
 269
 270 }
 271
 272 //_____________________________________________________________________________
 273 void AliTRDalignment::SetChZero()
 274 {
 275   //
 276   // reset to zero chamber data
 277   //
 278
 279   memset(&fCh[0][0],0,sizeof(fCh));
 280
 281 }
 282
 283 //_____________________________________________________________________________
 284 void AliTRDalignment::SetSmRandom(double a[6])
 285 {
 286   //
 287   // generate random gaussian supermodule data with sigmas a
 288   //
 289
 290   double x[6];
 291
 292   for (int i = 0; i < 18; i++) {
 293     fRan.Rannor(x[0],x[1]);
 294     fRan.Rannor(x[2],x[3]);
 295     fRan.Rannor(x[4],x[5]);
 296     for (int j = 0; j < 6; j++) x[j] *= a[j];
 297     SetSm(i,x);
 298     //PrintSm(i);
 299   }
 300
 301 }
 302
 303 //_____________________________________________________________________________
 304 void AliTRDalignment::SetChRandom(double a[6])
 305 {
 306   //
 307   // generate random gaussian chamber data with sigmas a
 308   //
 309
 310   double x[6];
 311
 312   for (int i = 0; i < 540; i++) {
 313     fRan.Rannor(x[0],x[1]);
 314     fRan.Rannor(x[2],x[3]);
 315     fRan.Rannor(x[4],x[5]);
 316     for (int j = 0; j < 6; j++) x[j] *= a[j];
 317     SetCh(i,x);
 318     //PrintCh(i);
 319   }
 320
 321 }
 322
 323 //_____________________________________________________________________________
 324 void AliTRDalignment::SetSmFull()
 325 {
 326   //
 327   // generate random gaussian supermodule data similar to the misalignment
 328   // expected from the mechanical precision
 329   //
 330
 331   double a[6];
 332
 333   a[0] = 0.3; // phi
 334   a[1] = 0.3; // z
 335   a[2] = 0.3; // r
 336   a[3] = 0.4/1000.0 / TMath::Pi()*180.0; // phi
 337   a[4] = 2.0/1000.0 / TMath::Pi()*180.0; // z
 338   a[5] = 0.4/1000.0 / TMath::Pi()*180.0; // r
 339
 340   SetSmRandom(a);
 341
 342 }
 343
 344 //_____________________________________________________________________________
 345 void AliTRDalignment::SetChFull()
 346 {
 347   //
 348   // generate random gaussian chamber data similar to the misalignment
 349   // expected from the mechanical precision
 350   //
 351
 352   double a[6];
 353
 354   a[0] = 0.1; // phi
 355   a[1] = 0.1; // z
 356   a[2] = 0.1; // r
 357   a[3] = 1.0/1000.0 / TMath::Pi()*180.0; // phi
 358   a[4] = 1.0/1000.0 / TMath::Pi()*180.0; // z
 359   a[5] = 0.7/1000.0 / TMath::Pi()*180.0; // r
 360
 361   SetChRandom(a);
 362
 363 }
 364
 365 //_____________________________________________________________________________
 366 void AliTRDalignment::SetSmResidual()
 367 {
 368   //
 369   // generate random gaussian supermodule data similar to the misalignment
 370   // remaining after full calibration
 371   // I assume that it will be negligible
 372   //
 373
 374   SetSmZero();
 375
 376 }
 377
 378 //_____________________________________________________________________________
 379 void AliTRDalignment::SetChResidual()
 380 {
 381   //
 382   // generate random gaussian chamber data similar to the misalignment
 383   // remaining after full calibration
 384   //
 385
 386   double a[6];
 387
 388   a[0] = 0.002; // phi
 389   a[1] = 0.003; // z
 390   a[2] = 0.007; // r
 391   a[3] = 0.3/1000.0 / TMath::Pi()*180.0; // phi
 392   a[4] = 0.3/1000.0 / TMath::Pi()*180.0; // z
 393   a[5] = 0.1/1000.0 / TMath::Pi()*180.0; // r
 394
 395   SetChRandom(a);
 396
 397 }
 398
 399 //_____________________________________________________________________________
 400 void AliTRDalignment::PrintSm(int i, FILE *fp) const
 401 {
 402   //
 403   // print the supermodule data
 404   //
 405
 406   fprintf(fp,"%4d   %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f   %6d  %s\n"
 407          ,i,fSm[i][0],fSm[i][1],fSm[i][2],fSm[i][3],fSm[i][4],fSm[i][5]
 408          ,0,GetSmName(i));
 409
 410 }
 411
 412 //_____________________________________________________________________________
 413 void AliTRDalignment::PrintCh(int i, FILE *fp) const
 414 {
 415   //
 416   // print the chamber data
 417   //
 418
 419   fprintf(fp,"%4d   %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f   %6d  %s\n"
 420          ,i,fCh[i][0],fCh[i][1],fCh[i][2],fCh[i][3],fCh[i][4],fCh[i][5]
 421          ,GetVoi(i),GetChName(i));
 422
 423 }
 424
 425 //_____________________________________________________________________________
 426 void AliTRDalignment::ReadAscii(char *filename)
 427 {
 428   //
 429   // read the alignment data from ascii file
 430   //
 431
 432   double x[6];      // alignment data
 433   int volid;        // volume id
 434   std::string syna; // symbolic name
 435   int j;            // dummy index
 436
 437   fstream fi(filename,fstream::in);
 438   if (!fi) {
 439     AliError(Form("cannot open input file %s",filename));
 440     return;
 441   }
 442
 443   // supermodules
 444
 445   for (int i = 0; i < 18; i++) {
 446     fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
 447     if (j != i) AliError(Form("sm %d expected, %d found",i,j));
 448     if (volid != 0) AliError(Form("sm %d volume id %d expected, %d found",i,0,volid));
 449     std::string symnam = GetSmName(i);
 450     if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data()));
 451     SetSm(i,x);
 452   }
 453
 454   // chambers
 455
 456   for (int i = 0; i < 540; i++) {
 457     fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
 458     if (j != i) AliError(Form("ch %d expected, %d found",i,j));
 459     if (volid != GetVoi(i)) AliError(Form("ch %d volume id %d expected, %d found",i,GetVoi(i),volid));
 460     std::string symnam = GetChName(i);
 461     if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data()));
 462     SetCh(i,x);
 463   }
 464
 465   fi.close();
 466
 467 }
 468
 469 //_____________________________________________________________________________
 470 void AliTRDalignment::ReadCurrentGeo()
 471 {
 472   //
 473   // use currently loaded geometry to determine misalignment by comparing
 474   // original and misaligned matrix of the last node
 475   // Now, original, does not mean "ideal". It is the matrix before the alignment.
 476   // So, if alignment was applied more than once, the numbers extracted will
 477   // represent just the last alignment. -- check this!
 478   //
 479
 480   TGeoHMatrix *ideSm[18];  // ideal
 481   TGeoHMatrix *misSm[18];  // misaligned
 482   for (int i = 0; i < 18; i++) {
 483
 484     // read misaligned and original matrices
 485
 486     TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetSmName(i));
 487     if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
 488     if (!pne) continue;
 489     TGeoPhysicalNode *node = pne->GetPhysicalNode();
 490     if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
 491     if (!node) continue;
 492     misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
 493     ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
 494
 495     // calculate the local misalignment matrices as inverse misaligned times ideal
 496
 497     TGeoHMatrix mat(ideSm[i]->Inverse());
 498     mat.Multiply(misSm[i]);
 499     double *tra = mat.GetTranslation();
 500     double *rot = mat.GetRotationMatrix();
 501     double pars[6];
 502     pars[0] = tra[0];
 503     pars[1] = tra[1];
 504     pars[2] = tra[2];
 505     if (TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) AliError("Failed to extract roll-pitch-yall angles!");
 506     double raddeg = TMath::RadToDeg();
 507     pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
 508     pars[4] = raddeg * TMath::ASin(rot[2]);
 509     pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
 510     SetSm(i,pars);
 511
 512     // cleanup
 513
 514     delete ideSm[i];
 515     delete misSm[i];
 516   }
 517
 518   TGeoHMatrix *ideCh[540]; // ideal
 519   TGeoHMatrix *misCh[540]; // misaligned
 520   for (int i = 0; i < 540; i++) {
 521
 522     // read misaligned and original matrices
 523
 524     TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetChName(i));
 525     if (!pne) AliError(Form("no such physical node entry: %s",GetChName(i)));
 526     if (!pne) continue;
 527     TGeoPhysicalNode *node = pne->GetPhysicalNode();
 528     if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i)));
 529     if (!node) continue;
 530     misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
 531     ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
 532
 533     // calculate the local misalignment matrices as inverse misaligned times ideal
 534
 535     TGeoHMatrix mat(ideCh[i]->Inverse());
 536     mat.Multiply(misCh[i]);
 537     double *tra = mat.GetTranslation();
 538     double *rot = mat.GetRotationMatrix();
 539     double pars[6];
 540     pars[0] = tra[0];
 541     pars[1] = tra[1];
 542     pars[2] = tra[2];
 543     if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
 544       AliError("Failed to extract roll-pitch-yall angles!");
 545       return;
 546     }
 547     double raddeg = TMath::RadToDeg();
 548     pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
 549     pars[4] = raddeg * TMath::ASin(rot[2]);
 550     pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
 551     SetCh(i,pars);
 552
 553     // cleanup
 554     delete ideCh[i];
 555     delete misCh[i];
 556   }
 557
 558   return;
 559
 560 }
 561
 562 //_____________________________________________________________________________
 563 void AliTRDalignment::ReadRoot(char *filename)
 564 {
 565   //
 566   // read the alignment data from root file
 567   //
 568
 569   TFile fi(filename,"READ");
 570
 571   if (fi.IsOpen()) {
 572     TClonesArray *ar = (TClonesArray*) fi.Get("TRDAlignObjs");
 573     ArToNumbers(ar);
 574     fi.Close();
 575   }
 576   else AliError(Form("cannot open input file %s",filename));
 577
 578   return;
 579
 580 }
 581
 582 //_____________________________________________________________________________
 583 void AliTRDalignment::ReadDB(char *filename)
 584 {
 585   //
 586   // read the alignment data from database file
 587   //
 588
 589   TFile fi(filename,"READ");
 590
 591   if (fi.IsOpen()) {
 592     AliCDBEntry  *e  = (AliCDBEntry *) fi.Get("AliCDBEntry");
 593     e->PrintMetaData();
 594     fComment.SetString(e->GetMetaData()->GetComment());
 595     TClonesArray *ar = (TClonesArray *) e->GetObject();
 596     ArToNumbers(ar);
 597     fi.Close();
 598   }
 599   else AliError(Form("cannot open input file %s",filename));
 600
 601   return;
 602
 603 }
 604
 605 //_____________________________________________________________________________
 606 void AliTRDalignment::ReadDB(char *db, char *path, int run
 607                            , int version, int subversion)
 608 {
 609   //
 610   // read the alignment data from database
 611   //
 612
 613   AliCDBManager *cdb     = AliCDBManager::Instance();
 614   AliCDBStorage *storLoc = cdb->GetStorage(db);
 615   AliCDBEntry   *e       = storLoc->Get(path,run,version,subversion);
 616   if (e) {
 617     e->PrintMetaData();
 618     fComment.SetString(e->GetMetaData()->GetComment());
 619     TClonesArray  *ar      =  (TClonesArray *) e->GetObject();
 620     ArToNumbers(ar);
 621   }
 622 }
 623
 624 //_____________________________________________________________________________
 625 void AliTRDalignment::ReadSurveyReport(char *filename)
 626 {
 627   //
 628   // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ,
 629   // and fSurveyE.  Store the survey info in the fComment.
 630   // Each supermodule has 8 survey points. The point names look like
 631   // TRD_sm08ah0 and have the following meaning.
 632   //
 633   // sm00..17 mean supermodule 0 through 17, following the phi.
 634   // Supermodule 00 is between phi=0 and phi=20 degrees.
 635   //
 636   // a or c denotes the anticlockwise and clockwise end of the supermodule
 637   // in z. Clockwise end is where z is negative and where the muon arm sits.
 638   //
 639   // l or h denote low radius and high radius holes
 640   //
 641   // 0 or 1 denote the hole at smaller and at larger phi, respectively.
 642   //
 643
 644   // read the survey file
 645
 646   fstream in(filename,fstream::in);
 647   if (!in) {
 648     AliError(Form("cannot open input file %s",filename));
 649     return;
 650   }
 651
 652   // loop through the lines of the file until the beginning of data
 653
 654   TString title,date,subdetector,url,version,observations,system,units;
 655   while (1) {
 656     char pee=in.peek();
 657     if (pee==EOF) break;
 658     TString line;
 659     line.ReadLine(in);
 660     if (line.Contains("Title:"))        title.ReadLine(in);
 661     if (line.Contains("Date:"))         date.ReadLine(in);
 662     if (line.Contains("Subdetector:"))  subdetector.ReadLine(in);
 663     if (line.Contains("URL:"))          url.ReadLine(in);
 664     if (line.Contains("Version:"))      version.ReadLine(in);
 665     if (line.Contains("Observations:")) observations.ReadLine(in);
 666     if (line.Contains("System:"))       system.ReadLine(in);
 667     if (line.Contains("Units:"))        units.ReadLine(in);
 668     if (line.Contains("Data:"))         break;
 669   }
 670
 671   // check what we found so far (watch out, they have \r at the end)
 672
 673   std::cout<<"title .........."<<title<<std::endl;
 674   std::cout<<"date ..........."<<date<<std::endl;
 675   std::cout<<"subdetector ...."<<subdetector<<std::endl;
 676   std::cout<<"url ............"<<url<<std::endl;
 677   std::cout<<"version ........"<<version<<std::endl;
 678   std::cout<<"observations ..."<<observations<<std::endl;
 679   std::cout<<"system ........."<<system<<std::endl;
 680   std::cout<<"units .........."<<units<<std::endl;
 681
 682   if (!subdetector.Contains("TRD")) {
 683     AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
 684     return;
 685   }
 686   double tocm = 0; // we want to have it in cm
 687   if (units.Contains("mm"))      tocm = 0.1;
 688   else if (units.Contains("cm")) tocm = 1.0;
 689   else if (units.Contains("m"))  tocm = 100.0;
 690   else if (units.Contains("pc")) tocm = 3.24078e-15;
 691   else {
 692     AliError(Form("unexpected units: %s",units.Data()));
 693     return;
 694   }
 695   if (!system.Contains("ALICEPH")) {
 696     AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
 697     return;
 698   }
 699
 700   // scan the rest of the file which should contain list of surveyed points
 701   // for every point, decode the point name and store the numbers in the right
 702   // place in the arrays fSurveyX etc.
 703
 704   while (1) {
 705     TString pna; // point name
 706     char type;
 707     double x,y,z,precision;
 708     in >> pna >> x >> y >> z >> type >> precision;
 709     if (in.fail()) break;
 710     if (pna(0,6)!="TRD_sm") {
 711       AliError(Form("unexpected point name: %s",pna.Data()));
 712       break;
 713     }
 714     int i = atoi(pna(6,0).Data()); // supermodule number
 715     int j = -1;
 716     if (pna(8) == 'a') j=0; // anticlockwise, positive z
 717     if (pna(8) == 'c') j=1; // clockwise, negative z
 718     int k = -1;
 719     if (pna(9) == 'l') k=0; // low radius
 720     if (pna(9) == 'h') k=1; // high radius
 721     int l = atoi(pna(10,0).Data()); // phi within supermodule
 722     if (i>=0 && i<18 && j>=0 && j<2 && k>=0 && k<2 && l>=0 && l<2) {
 723       fSurveyX[i][j][k][l] = tocm*x;
 724       fSurveyY[i][j][k][l] = tocm*y;
 725       fSurveyZ[i][j][k][l] = tocm*z;
 726       fSurveyE[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
 727       std::cout << "decoded "<<pna<<" "
 728                 <<fSurveyX[i][j][k][l]<<" "
 729                 <<fSurveyY[i][j][k][l]<<" "
 730                 <<fSurveyZ[i][j][k][l]<<" "
 731                 <<fSurveyE[i][j][k][l]<<std::endl;
 732     } else AliError(Form("cannot decode point name: %s",pna.Data()));
 733   }
 734   in.close();
 735   TString info = "Survey "+title+" "+date+" "+url+" "+version+" "+observations;
 736   info.ReplaceAll("\r","");
 737   fComment.SetString(info.Data());
 738
 739 }
 740
 741 //_____________________________________________________________________________
 742 double AliTRDalignment::SurveyChi2(int i, double *a) {
 743
 744   //
 745   // Compare the survey results to the ideal positions of the survey marks
 746   // in the local frame of supermodule. When transforming, use the alignment
 747   // parameters a[6]. Return chi-squared.
 748   //
 749
 750   if (!IsGeoLoaded()) return 0;
 751   printf("Survey of supermodule %d\n",i);
 752   AliAlignObjAngles al(GetSmName(i),0,a[0],a[1],a[2],a[3],a[4],a[5],0);
 753   TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetSmName(i));
 754   if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
 755   TGeoPhysicalNode *node = pne->GetPhysicalNode();
 756   if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
 757
 758   //  al.ApplyToGeometry();
 759   //  node = pne->GetPhysicalNode(); // changed in the meantime
 760   //  TGeoHMatrix *ma = node->GetMatrix();
 761
 762   // a less destructive method (it does not modify geometry), gives the same result:
 763
 764   TGeoHMatrix *ma = new TGeoHMatrix();
 765   al.GetLocalMatrix(*ma);
 766   ma->MultiplyLeft(node->GetMatrix()); // global trafo, modified by a[]
 767
 768   double chi2=0;
 769   printf("              sm   z   r  phi    x (lab phi)  y (lab z)   z (lab r)   all in cm\n");
 770   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
 771     if (fSurveyE[i][j][k][l] == 0.0) continue; // no data for this survey point
 772     double master[3] = {fSurveyX[i][j][k][l],fSurveyY[i][j][k][l],fSurveyZ[i][j][k][l]};
 773     double local[3];
 774     ma->MasterToLocal(master,local);
 775     double dx = local[0]-fSurveyX0[j][k][l];
 776     double dy = local[1]-fSurveyY0[j][k][l];
 777     double dz = local[2]-fSurveyZ0[j][k][l];
 778     chi2 += (dx*dx+dy*dy+dz*dz)/fSurveyE[i][j][k][l]/fSurveyE[i][j][k][l];
 779     printf("local survey %3d %3d %3d %3d %12.3f %12.3f %12.3f\n",i,j,k,l,local[0],local[1],local[2]);
 780     printf("local ideal                  %12.3f %12.3f %12.3f\n",fSurveyX0[j][k][l],
 781            fSurveyY0[j][k][l],fSurveyZ0[j][k][l]);
 782     printf("difference                   %12.3f %12.3f %12.3f\n",dx,dy,dz);
 783   }
 784   printf("chi2 = %.2f\n",chi2);
 785   return chi2;
 786 }
 787
 788 //_____________________________________________________________________________
 789 void trdAlignmentFcn(int &npar, double *g, double &f, double *par, int iflag) {
 790
 791   //
 792   // Standard function as needed by Minuit-like minimization procedures.
 793   // For the set of parameters par calculates and returns chi-squared.
 794   //
 795
 796   // smuggle a C++ object into a C function
 797   AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit();
 798
 799   f = alignment->SurveyChi2(par);
 800   if (iflag==3);
 801   if (npar);
 802   if (g); // no warnings about unused stuff...
 803
 804 }
 805
 806 //_____________________________________________________________________________
 807 void AliTRDalignment::SurveyToAlignment(int i,char *flag) {
 808
 809   //
 810   // Find the supermodule alignment parameters needed to make the survey
 811   // results coincide with the ideal positions of the survey marks.
 812   // The string flag should look like "101000"; the six characters corresponds
 813   // to the six alignment parameters and 0/1 mean that the parameter should
 814   // be fixed/released in the fit.
 815
 816   if (strlen(flag)!=6) {
 817     AliError(Form("unexpected flag: %s",flag));
 818     return;
 819   }
 820
 821   printf("Finding alignment matrix for supermodule %d\n",i);
 822   fIbuffer[0] = i; // store the sm number in the buffer so minuit can see it
 823
 824   TFitter fitter(100);
 825   gMinuit->SetObjectFit(this);
 826   fitter.SetFCN(trdAlignmentFcn);
 827   fitter.SetParameter(0,"dx",0,0.5,0,0);
 828   fitter.SetParameter(1,"dy",0,0.5,0,0);
 829   fitter.SetParameter(2,"dz",0,0.5,0,0);
 830   fitter.SetParameter(3,"rx",0,0.1,0,0);
 831   fitter.SetParameter(4,"ry",0,0.1,0,0);
 832   fitter.SetParameter(5,"rz",0,0.1,0,0);
 833
 834   for (int j=0; j<6; j++) if (flag[j]=='0') fitter.FixParameter(j);
 835
 836   double arglist[100];
 837   arglist[0] = 2;
 838   fitter.ExecuteCommand("SET PRINT", arglist, 1);
 839   fitter.ExecuteCommand("SET ERR", arglist, 1);
 840   arglist[0]=50;
 841   //fitter.ExecuteCommand("SIMPLEX", arglist, 1);
 842   fitter.ExecuteCommand("MINIMIZE", arglist, 1);
 843   fitter.ExecuteCommand("CALL 3", arglist,0);
 844   double a[6];
 845   for (int j=0; j<6; j++) a[j] = fitter.GetParameter(j);
 846   SetSm(i,a);
 847   for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParameter(j));
 848   printf("\n");
 849   for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParError(j));
 850   printf("\n");
 851
 852 }
 853
 854 //_____________________________________________________________________________
 855 void AliTRDalignment::ReadAny(char *filename)
 856 {
 857   //
 858   // read the alignment data from any kind of file
 859   //
 860
 861   TString fist(filename);
 862   if (fist.EndsWith(".txt")) ReadAscii(filename);
 863   if (fist.EndsWith(".dat")) ReadAscii(filename);
 864   if (fist.EndsWith(".root")) {
 865     if (fist.Contains("Run")) ReadDB(filename);
 866     else ReadRoot(filename);
 867   }
 868
 869 }
 870
 871 //_____________________________________________________________________________
 872 void AliTRDalignment::WriteAscii(char *filename) const
 873 {
 874   //
 875   // store the alignment data on ascii file
 876   //
 877
 878   FILE *fp = fopen(filename, "w");
 879   if (!fp) {
 880     AliError(Form("cannot open output file %s",filename));
 881     return;
 882   }
 883
 884   PrintSm(fp);
 885   PrintCh(fp);
 886
 887   fclose(fp);
 888
 889 }
 890
 891 //_____________________________________________________________________________
 892 void AliTRDalignment::WriteRoot(char *filename)
 893 {
 894   //
 895   // store the alignment data on root file
 896   //
 897
 898   TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
 899   NumbersToAr(ar);
 900   TFile fo(filename,"RECREATE");
 901   if (fo.IsOpen()) {
 902     fo.cd();
 903     fo.WriteObject(ar,"TRDAlignObjs","kSingleKey");
 904     fo.Close();
 905   }
 906   else AliError(Form("cannot open output file %s",filename));
 907
 908   delete ar;
 909
 910 }
 911
 912 //_____________________________________________________________________________
 913 void AliTRDalignment::WriteDB(char *filename, int run0, int run1)
 914 {
 915   //
 916   // dumping on a DB-like file
 917   //
 918
 919   TClonesArray   *ar = new TClonesArray("AliAlignObjAngles",10000);
 920   NumbersToAr(ar);
 921   char *path = "TRD/Align/Data";
 922   AliCDBId id(path,run0,run1);
 923   AliCDBMetaData *md = new AliCDBMetaData();
 924   md->SetResponsible("Dariusz Miskowiec");
 925   md->SetComment(fComment.GetString().Data());
 926   AliCDBEntry    *e  = new AliCDBEntry(ar, id, md);
 927   TFile fi(filename,"RECREATE");
 928   if (fi.IsOpen()) {
 929     e->Write();
 930     fi.Close();
 931   }
 932   else AliError(Form("cannot open input file %s",filename));
 933
 934   delete e;
 935   delete md;
 936   delete ar;
 937
 938   return;
 939
 940 }
 941
 942 //_____________________________________________________________________________
 943 void AliTRDalignment::WriteDB(char *db, char *path, int run0, int run1)
 944 {
 945   //
 946   // store the alignment data in database
 947   //
 948
 949   TClonesArray   *ar      = new TClonesArray("AliAlignObjAngles",10000);
 950   NumbersToAr(ar);
 951   AliCDBManager  *cdb     = AliCDBManager::Instance();
 952   AliCDBStorage  *storLoc = cdb->GetStorage(db);
 953   AliCDBMetaData *md      = new AliCDBMetaData();
 954   md->SetResponsible("Dariusz Miskowiec");
 955   md->SetComment(fComment.GetString().Data());
 956   AliCDBId id(path,run0,run1);
 957   storLoc->Put(ar,id,md);
 958   md->Delete();
 959   delete ar;
 960
 961 }
 962
 963 //_____________________________________________________________________________
 964 void AliTRDalignment::WriteGeo(char *filename)
 965 {
 966   //
 967   // apply misalignment to current geometry and store the
 968   // resulting geometry on a root file
 969   //
 970
 971   TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
 972   NumbersToAr(ar);
 973   delete ar;
 974   gGeoManager->Export(filename);
 975
 976 }
 977
 978 //_____________________________________________________________________________
 979 double AliTRDalignment::GetSmRMS(int xyz) const
 980 {
 981   //
 982   // rms fSm[][xyz]
 983   //
 984
 985   double s1 = 0.0;
 986   double s2 = 0.0;
 987   for (int i = 0; i < 18; i++) {
 988     s1 += fSm[i][xyz];
 989     s2 += fSm[i][xyz]*fSm[i][xyz];
 990   }
 991   double rms2 = s2/18.0 - s1*s1/18.0/18.0;
 992
 993   return rms2>0 ? sqrt(rms2) : 0.0;
 994
 995 }
 996
 997 //_____________________________________________________________________________
 998 double AliTRDalignment::GetChRMS(int xyz) const
 999 {
1000   //
1001   // rms fCh[][xyz]
1002   //
1003
1004   double s1 =0.0;
1005   double s2 =0.0;
1006   for (int i = 0; i < 540; i++) {
1007     s1 += fCh[i][xyz];
1008     s2 += fCh[i][xyz]*fCh[i][xyz];
1009   }
1010   double rms2 = s2/540.0 - s1*s1/540.0/540.0;
1011
1012   return rms2>0 ? sqrt(rms2) : 0.0;
1013
1014 }
1015
1016 //_____________________________________________________________________________
1017 void AliTRDalignment::PrintSmRMS() const
1018 {
1019   //
1020   // dump rms of fSm
1021   //
1022
1023   printf("       %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f  supermodule rms\n"
1024         ,GetSmRMS(0),GetSmRMS(1),GetSmRMS(2),GetSmRMS(3),GetSmRMS(4),GetSmRMS(5));
1025
1026 }
1027
1028 //_____________________________________________________________________________
1029 void AliTRDalignment::PrintChRMS() const
1030 {
1031   //
1032   // dump rms of fCh
1033   //
1034
1035   printf("       %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f  chamber rms\n"
1036         ,GetChRMS(0),GetChRMS(1),GetChRMS(2),GetChRMS(3),GetChRMS(4),GetChRMS(5));
1037
1038 }
1039
1040 //_____________________________________________________________________________
1041 void AliTRDalignment::ArToNumbers(TClonesArray *ar)
1042 {
1043   //
1044   // for each of the alignment objects in array ar extract the six local
1045   // alignment parameters; recognize by name to which supermodule or chamber
1046   // the alignment object pertains; set the respective fSm or fCh
1047   //
1048
1049   ar->Sort();
1050   if (!IsGeoLoaded()) return;
1051   for (int i = 0; i < ar->GetEntries(); i++) {
1052     AliAlignObj *aao = (AliAlignObj *) ar->At(i);
1053     aao->ApplyToGeometry();
1054   }
1055   SetZero();
1056   ReadCurrentGeo();
1057
1058 }
1059
1060 //_____________________________________________________________________________
1061 void AliTRDalignment::NumbersToAr(TClonesArray *ar)
1062 {
1063   //
1064   // build array of AliAlignObj objects based on fSm and fCh data
1065   // at the same time, apply misalignment to the currently loaded geometry
1066   // it is important to apply misalignment of supermodules before creating
1067   // alignment objects for chambers
1068   //
1069
1070   if (!IsGeoLoaded()) return;
1071   TClonesArray &alobj = *ar;
1072   int nobj = 0;
1073   for (int i = 0; i <  18; i++) {
1074       new(alobj[nobj]) AliAlignObjAngles(GetSmName(i)
1075                                         ,0
1076                                         ,fSm[i][0],fSm[i][1],fSm[i][2]
1077                                         ,fSm[i][3],fSm[i][4],fSm[i][5]
1078                                         ,0);
1079     ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1080     nobj++;
1081   }
1082
1083   for (int i = 0; i < 540; i++) {
1084     new(alobj[nobj]) AliAlignObjAngles(GetChName(i)
1085                                       ,GetVoi(i)
1086                                       ,fCh[i][0],fCh[i][1],fCh[i][2]
1087                                       ,fCh[i][3],fCh[i][4],fCh[i][5]
1088                                       ,0);
1089     ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1090     nobj++;
1091   }
1092   AliInfo("current geometry modified");
1093
1094 }
1095
1096 //_____________________________________________________________________________
1097 int AliTRDalignment::IsGeoLoaded()
1098 {
1099   //
1100   // check whether a geometry is loaded
1101   // issue a warning if geometry is not ideal
1102   //
1103
1104   if (gGeoManager) {
1105     if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) AliWarning("current geometry is not ideal");
1106     return 1;
1107   } else {
1108     AliError("first load geometry by calling TGeoManager::Import(filename)");
1109     return 0;
1110   }
1111
1112 }
1113
1114 //_____________________________________________________________________________