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