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