FMD/AliFMDSurveyToAlignObjs.cxx

   1 //
   2 // Class to take survey data and
   3 // transform that to alignment objects.
   4 //
   5 // FMD
   6 //
   7 #include "AliFMDSurveyToAlignObjs.h"
   8 #include "AliLog.h"
   9 #include "AliSurveyPoint.h"
  10 #include <TGraph2DErrors.h>
  11 #include <TF2.h>
  12 #include <TVector3.h>
  13 #include <iostream>
  14 #include <iomanip>
  15 #include <TMath.h>
  16 #include <TRotation.h>
  17 #include <TGeoMatrix.h>
  18 #include <TGeoManager.h>
  19 #include <TGeoPhysicalNode.h>
  20 #include "AliFMDGeometry.h"
  21
  22 //____________________________________________________________________
  23 Double_t
  24 AliFMDSurveyToAlignObjs::GetUnitFactor() const
  25 {
  26   // Returns the conversion factor from the measured values to
  27   // centimeters.
  28   if (!fSurveyObj) return 0;
  29   TString units(fSurveyObj->GetUnits());
  30   if      (units.CompareTo("mm", TString::kIgnoreCase) == 0) return .1;
  31   else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) return 1.;
  32   else if (units.CompareTo("m",  TString::kIgnoreCase) == 0) return 100.;
  33   return 1;
  34 }
  35
  36 //____________________________________________________________________
  37 Bool_t
  38 AliFMDSurveyToAlignObjs::GetPoint(const char* name,
  39                                   TVector3&   point,
  40                                   TVector3&   error) const
  41 {
  42   // Get named point.   On return, point will contain the point
  43   // coordinates in centimeters, and error will contain the
  44   // meassurement errors in centimeters too.  If no point is found,
  45   // returns false, otherwise true.
  46   if (!fSurveyPoints) return kFALSE;
  47
  48   Double_t unit = GetUnitFactor();
  49   if (unit == 0) return kFALSE;
  50
  51   TObject* obj  = fSurveyPoints->FindObject(name);
  52   if (!obj) return kFALSE;
  53
  54   AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
  55   point.SetXYZ(unit * p->GetX(),
  56                unit * p->GetY(),
  57                unit * p->GetZ());
  58   error.SetXYZ(unit * p->GetPrecisionX(),
  59                unit * p->GetPrecisionY(),
  60                unit * p->GetPrecisionZ());
  61   return kTRUE;
  62 }
  63
  64 //____________________________________________________________________
  65 Bool_t
  66 AliFMDSurveyToAlignObjs::CalculatePlane(const     TVector3& a,
  67                                         const     TVector3& b,
  68                                         const     TVector3& c,
  69                                         Double_t  depth,
  70                                         Double_t* trans,
  71                                         Double_t* rot) const
  72 {
  73   //
  74   // Calculate the plane translation and rotation from 3 survey points
  75   //
  76   // Parameters:
  77   //    a     1st Survey point
  78   //    b     2nd Survey point
  79   //    c     3rd Survey point
  80   //    trans Translation vector
  81   //    rot   Rotation matrix (direction cosines)
  82   //
  83   // Return:
  84   //
  85   //
  86
  87   // Vector a->b, b->c, and normal to plane defined by these two
  88   // vectors.
  89   TVector3 ab(b-a), bc(c-b);
  90
  91   // Normal vector to the plane of the fiducial marks obtained
  92   // as cross product of the two vectors on the plane d0^d1
  93   TVector3 nn(ab.Cross(bc));
  94   if (nn.Mag() < 1e-8) {
  95     Info("CalculatePlane", "Normal vector is null vector");
  96     return kFALSE;
  97   }
  98
  99   // We express the plane in Hessian normal form.
 100   //
 101   //   n x = -p,
 102   //
 103   // where n is the normalised normal vector given by
 104   //
 105   //   n_x = a / l,   n_y = b / l,   n_z = c / l,  p = d / l
 106   //
 107   // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the
 108   // normal plane equation
 109   //
 110   //  ax + by + cz + d = 0
 111   //
 112   // Normalize
 113   TVector3 n(nn.Unit());
 114   // Double_t p = - (n * a);
 115
 116   // The center of the square with the fiducial marks as the
 117   // corners.  The mid-point of one diagonal - md.  Used to get the
 118   // center of the surveyd box.
 119   TVector3 md(a + c);
 120   md *= 1/2.;
 121
 122   // The center of the box.
 123   TVector3 orig(md - depth * n);
 124   trans[0] = orig[0];
 125   trans[1] = orig[1];
 126   trans[2] = orig[2];
 127
 128   // Normalize the spanning vectors
 129   TVector3 uab(ab.Unit());
 130   TVector3 ubc(bc.Unit());
 131
 132   for (size_t i = 0; i < 3; i++) {
 133     rot[i * 3 + 0] = ubc[i];
 134     rot[i * 3 + 1] = uab[i];
 135     rot[i * 3 + 2] = n[i];
 136   }
 137   return kTRUE;
 138 }
 139
 140 //____________________________________________________________________
 141 Bool_t
 142 AliFMDSurveyToAlignObjs::FitPlane(const TObjArray& points,
 143                                   const TObjArray& errors,
 144                                   Double_t         /* depth */,
 145                                   Double_t*        trans,
 146                                   Double_t*        rot) const
 147 {
 148   //
 149   // Calculate the plane rotation and translation by doing a fit of
 150   // the plane equation to the surveyed points.  At least 4 points
 151   // must be passed in the @a points array with corresponding errors
 152   // in the array @a errors.  The arrays are assumed to contain
 153   // TVector3 objects.
 154   //
 155   // Parameters:
 156   //    points Array surveyed positions
 157   //    errors Array of errors corresponding to @a points
 158   //    depth  Survey targets depth (perpendicular to the plane)
 159   //    trans  On return, translation of the plane
 160   //    rot    On return, rotation (direction cosines) of the plane
 161   //
 162   // Return:
 163   //    @c true on success, @c false otherwise
 164   //
 165
 166   Int_t nPoints = points.GetEntries();
 167   if (nPoints < 4) {
 168     AliError(Form("Cannot fit a plane equation to less than 4 survey points, "
 169                   "got only %d", nPoints));
 170     return kFALSE;
 171   }
 172
 173   TGraph2DErrors g;
 174   // Loop and fill graph
 175   for (int i = 0; i < nPoints; i++) {
 176     TVector3* p = static_cast<TVector3*>(points.At(i));
 177     TVector3* e = static_cast<TVector3*>(errors.At(i));
 178
 179     if (!p || !e) continue;
 180
 181     g.SetPoint(i, p->X(), p->Y(), p->Z());
 182     g.SetPointError(i, e->X(), e->Y(), e->Z());
 183   }
 184
 185   // Check that we have enough points
 186   if (g.GetN() < 4) {
 187     AliError(Form("Only got %d survey points - no good for plane fit",
 188                   g.GetN()));
 189     return kFALSE;
 190   }
 191
 192   // Next, declare fitting function and fit to graph.
 193   // Fit to the plane equation:
 194   //
 195   //   ax + by + cz + d = 0
 196   //
 197   // or
 198   //
 199   //   z = - ax/c - by/c - d/c
 200   //
 201   TF2 f("plane", "-[0]*x-[1]*y-[2]",
 202         g.GetXmin(), g.GetXmax(), g.GetYmin(), g.GetYmax());
 203   g.Fit(&f, "Q");
 204
 205   // Now, extract the normal and offset
 206   TVector3 nv(f.GetParameter(0), f.GetParameter(1), 1);
 207   TVector3 n(nv.Unit());
 208   Double_t p = -f.GetParameter(2);
 209
 210   // Create two vectors spanning the plane
 211   TVector3 a(1, 0, f.Eval(1, 0)-p);
 212   TVector3 b(0, -1, f.Eval(0, -1)-p);
 213   TVector3 ua(a.Unit());
 214   TVector3 ub(b.Unit());
 215   // Double_t angAb = ua.Angle(ub);
 216   // PrintVector("ua: ", ua);
 217   // PrintVector("ub: ", ub);
 218   // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
 219
 220   for (size_t i = 0; i < 3; i++) {
 221     rot[i * 3 + 0] = ua[i];
 222     rot[i * 3 + 1] = ub[i];
 223     rot[i * 3 + 2] = n[i];
 224   }
 225
 226   // The intersection of the plane is given by (0, 0, -d/c)
 227   trans[0] = 0;
 228   trans[1] = 0;
 229   trans[2] = p;
 230
 231   return kTRUE;
 232 }
 233
 234
 235 //____________________________________________________________________
 236 Bool_t
 237 AliFMDSurveyToAlignObjs::MakeDelta(const char*  path,
 238                                    const Double_t*    rot,
 239                                    const Double_t*    trans,
 240                                    TGeoHMatrix& delta) const
 241 {
 242   //
 243   // Create a delta transform from a global rotation matrix and
 244   // translation.
 245   //
 246   // Parameters:
 247   //    path   Path of element to transform.
 248   //    rot    Rotation matrix (direction cosines)
 249   //    trans  Translation
 250   //    delta  On return, the delta transform
 251   //
 252   // Return:
 253   //    Newly
 254   //
 255   if (!gGeoManager)           return kFALSE;
 256   if (!gGeoManager->cd(path)) return kFALSE;
 257
 258
 259   TGeoMatrix*  global = gGeoManager->GetCurrentMatrix();
 260 #if 0
 261   PrintRotation(Form("%s rot:", global->GetName()),global->GetRotationMatrix());
 262   PrintVector(Form("%s trans:", global->GetName()),global->GetTranslation());
 263 #endif
 264
 265   return MakeDelta(global, rot, trans, delta);
 266 }
 267
 268 //____________________________________________________________________
 269 Bool_t
 270 AliFMDSurveyToAlignObjs::MakeDelta(const TGeoMatrix*  global,
 271                                    const Double_t*    rot,
 272                                    const Double_t*    trans,
 273                                    TGeoHMatrix& delta) const
 274 {
 275   //
 276   // Create a delta transform from a global rotation matrix and
 277   // translation.
 278   //
 279   // Parameters:
 280   //    global Global matrix of element to transform.
 281   //    rot    Rotation matrix (direction cosines)
 282   //    trans  Translation
 283   //    delta  On return, the delta transform
 284   //
 285   // Return:
 286   //    Newly
 287   //
 288   TGeoHMatrix* geoM = new TGeoHMatrix;
 289   geoM->SetTranslation(trans);
 290   geoM->SetRotation(rot);
 291
 292   delta = global->Inverse();
 293   delta.MultiplyLeft(geoM);
 294
 295   return true;
 296 }
 297
 298 //____________________________________________________________________
 299 Bool_t
 300 AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const
 301 {
 302   //
 303   // Get the FMD1 plane from the survey points
 304   //
 305   // Parameters:
 306   //    rot    Rotation matrix (direction cosines)
 307   //    trans  Translation
 308   //
 309   // Return:
 310   //    @c true on success, @c false otherwise.
 311   //
 312
 313   // The possile survey points
 314   TVector3  icb, ict, ocb, oct, dummy;
 315   Int_t     missing = 0;
 316   if (!GetPoint("V0L_ICB", icb, dummy)) missing++;
 317   if (!GetPoint("V0L_ICT", ict, dummy)) missing++;
 318   if (!GetPoint("V0L_OCB", ocb, dummy)) missing++;
 319   if (!GetPoint("V0L_OCT", oct, dummy)) missing++;
 320
 321   // Check that we have enough points
 322   if (missing > 1) {
 323     AliWarning(Form("Only got %d survey points - no good for FMD1 plane",
 324                     4-missing));
 325     return kFALSE;
 326   }
 327
 328 #if 0
 329   const char* lidN = "FMD1_lid_mat0";
 330   TGeoMatrix* lidM = static_cast<TGeoMatrix*>(gGeoManager->GetListOfMatrices()
 331                                               ->FindObject(lidN));
 332   if (!lidM) {
 333     AliError(Form("Couldn't find FMD1 lid transformation %s", lidN));
 334     return kFALSE;
 335   }
 336
 337   const Double_t* lidT = lidM->GetTranslation();
 338   Double_t        lidZ = lidT[2];
 339   Double_t        off  = lidZ-3.3;
 340 #else
 341   Double_t        off  = 0;
 342 #endif
 343
 344   if (!CalculatePlane(ocb, icb, ict, off, trans, rot)) return kFALSE;
 345   // PrintRotation("FMD1 rotation:",  rot);
 346   // PrintVector("FMD1 translation:", trans);
 347
 348   return kTRUE;
 349 }
 350
 351 //____________________________________________________________________
 352 Bool_t
 353 AliFMDSurveyToAlignObjs::DoFMD1()
 354 {
 355   //
 356   // Do the FMD1 analysis.  We have 4 survey targets on V0-A on the
 357   // C-side.  These are
 358   //
 359   //  - V0A_ICT  In-side, C-side, top.
 360   //  - V0A_ICB  In-side, C-side, bottom.
 361   //  - V0A_OCT  Out-side, C-side, top.
 362   //  - V0A_OCB  Out-side, C-side, bottom.
 363   //
 364   // These 4 survey targets sit 3.3mm over the V0-A C-side surface, or
 365   // 3.3mm over the back surface of FMD1.
 366   //
 367   // Since these are really sitting on a plane, we can use the method
 368   // proposed by the CORE offline.
 369   //
 370   // Return:
 371   //    @c true on success, @c false otherwise.
 372   //
 373
 374   // Do the FMD1 stuff
 375   Double_t rot[9], trans[3];
 376   if (!GetFMD1Plane(rot, trans)) return kFALSE;
 377   // const char* path = "/ALIC_1/F1MT_1/FMD1_lid_0";
 378
 379   // TGeoHMatrix delta;
 380   TGeoTranslation global(0,0,324.670);
 381   if (!MakeDelta(&global, rot, trans, fFMD1Delta))
 382     return kFALSE;
 383
 384   // PrintRotation("FMD1 delta rotation:",  fFMD1Delta.GetRotationMatrix());
 385   // PrintVector("FMD1 delta translation:", fFMD1Delta.GetTranslation());
 386
 387   return kTRUE;
 388 }
 389
 390 //____________________________________________________________________
 391 Bool_t
 392 AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const
 393 {
 394   //
 395   // Get the surveyed plane corresponding to the backside of FMD2.
 396   // The plane is done as a best fit of the plane equation to at least
 397   // 4 of the available survey points.
 398   //
 399   // Parameters:
 400   //    rot    Rotation matrix (direction cosines)
 401   //    trans  Translation vector.
 402   //
 403   // Return:
 404   //    @c true on success, @c false otherwise
 405   //
 406
 407   // The possible survey points
 408   const char*    names[] = { "FMD2_ITOP",  "FMD2_OTOP",
 409                              "FMD2_IBOTM", "FMD2_OBOTM",
 410                              "FMD2_IBOT",  "FMD2_OBOT",
 411                              0 };
 412   const char**   name    = names;
 413
 414   TObjArray points;
 415   TObjArray errors;
 416
 417   // Loop and fill graph
 418   while (*name) {
 419     TVector3 p, e;
 420     if (!GetPoint(*name++, p, e)) continue;
 421
 422     points.Add(new TVector3(p));
 423     errors.Add(new TVector3(e));
 424   }
 425   if (points.GetEntries() < 4) {
 426     AliWarning(Form("Only got %d survey points - no good for FMD2 plane",
 427                     points.GetEntries()));
 428     return kFALSE;
 429   }
 430
 431   return FitPlane(points, errors, 0, trans, rot);
 432 }
 433
 434 #define M(I,J) rot[(J-1) * 3 + (I-1)]
 435 //____________________________________________________________________
 436 Bool_t
 437 AliFMDSurveyToAlignObjs::DoFMD2()
 438 {
 439   //
 440   // Do the FMD2 calculations.  We have 6 survey points of which only
 441   // 5 are normally surveyed.  These are all sittings
 442   //
 443   //  - FMD2_ITOP   - In-side, top
 444   //  - FMD2_IBOTM  - In-side, middle bottom
 445   //  - FMD2_IBOT   - In-side, bottom
 446   //  - FMD2_OTOP   - Out-side, top
 447   //  - FMD2_OBOTM  - Out-side, middle bottom
 448   //  - FMD2_OBOT   - Out-side, bottom
 449   //
 450   // The nominal coordinates of these retro-fitted survey stickers
 451   // isn't known.  Also, these stickers are put on a thin (0.3mm
 452   // thick) carbon cover which flexes quite easily.  This means, that
 453   // to rotations and xy-translation obtained from the survey data
 454   // cannot be used, and left is only the z-translation.
 455   //
 456   // Further more, since FMD2 to is attached to the ITS SPD thermal
 457   // screen, it is questionable if the FMD2 survey will ever be used.
 458   //
 459   // Return:
 460   //    @c true on success, @c false otherwise.
 461   //
 462
 463   // Do the FMD2 stuff
 464   Double_t rot[9], trans[3];
 465   if (!GetFMD2Plane(rot, trans)) return kFALSE;
 466   // PrintRotation("FMD2 rotation:",  rot);
 467   // PrintVector("FMD2 translation:", trans);
 468
 469   for (int i = 0; i < 3; i++) {
 470     for (int j = 0; j < 3; j++) {
 471       rot[i*3+j] = (i == j ? 1 : 0);
 472     }
 473   }
 474   trans[0] = trans[1] = 0;
 475   trans[2] += 0.015;
 476   // PrintRotation("FMD2 rotation:",  rot);
 477   // PrintVector("FMD2 translation:", trans);
 478
 479   // TGeoHMatrix delta;
 480   if (!MakeDelta("/ALIC_1/F2MT_2/FMD2_support_0/FMD2_back_cover_2",
 481                  rot, trans, fFMD2Delta)) return kFALSE;
 482
 483   // PrintRotation("FMD2 delta rotation:",  fFMD2Delta.GetRotationMatrix());
 484   // PrintVector("FMD2 delta translation:", fFMD2Delta.GetTranslation());
 485
 486   return kTRUE;
 487 }
 488
 489 //____________________________________________________________________
 490 void
 491 AliFMDSurveyToAlignObjs::Run()
 492 {
 493   //
 494   // Run the task.
 495   //
 496   //
 497
 498   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 499   geom->Init();
 500   geom->InitTransformations();
 501
 502   DoFMD1();
 503   DoFMD2();
 504 }
 505
 506 //____________________________________________________________________
 507 Bool_t
 508 AliFMDSurveyToAlignObjs::CreateAlignObjs()
 509 {
 510   //
 511   //
 512   // Method to create the alignment objects
 513   //
 514   // Return:
 515   //    @c true on success, @c false otherwise
 516   //
 517   TClonesArray& array = *fAlignObjArray;
 518   Int_t         n     = array.GetEntriesFast();
 519
 520   if (!fFMD1Delta.IsIdentity()) {
 521     new (array[n++]) AliAlignObjParams("FMD1/FMD1_T", 0, fFMD1Delta, kTRUE);
 522     new (array[n++]) AliAlignObjParams("FMD1/FMD1_B", 0, fFMD1Delta, kTRUE);
 523   }
 524   if (!fFMD2Delta.IsIdentity()) {
 525     new (array[n++]) AliAlignObjParams("FMD2/FMD2_T", 0, fFMD2Delta, kTRUE);
 526     new (array[n++]) AliAlignObjParams("FMD2/FMD2_B", 0, fFMD2Delta, kTRUE);
 527   }
 528   // array.Print();
 529
 530   return kTRUE;
 531 }
 532
 533 //____________________________________________________________________
 534 void
 535 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const TVector3& v)
 536 {
 537   //
 538   // Service member function to print a vector
 539   //
 540   // Parameters:
 541   //    text Prefix text
 542   //    v    Vector
 543   //
 544   Double_t va[] = { v.X(), v.Y(), v.Z() };
 545   PrintVector(text, va);
 546 }
 547 //____________________________________________________________________
 548 void
 549 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const Double_t* v)
 550 {
 551   //
 552   // Service member function to print a vector
 553   //
 554   // Parameters:
 555   //    text Prefix text
 556   //    v    Vector (array of 3 doubles)
 557   //
 558   std::cout << text
 559             << std::setw(15) << v[0]
 560             << std::setw(15) << v[1]
 561             << std::setw(15) << v[2]
 562             << std::endl;
 563 }
 564
 565
 566 //____________________________________________________________________
 567 void
 568 AliFMDSurveyToAlignObjs::PrintRotation(const char* text, const Double_t* rot)
 569 {
 570   //
 571   // Service member function to print a rotation matrix
 572   //
 573   // Parameters:
 574   //    text Prefix text
 575   //    v    Matrix (array of 9 doubles)
 576   //
 577
 578   std::cout << text << std::endl;
 579   for (size_t i = 0; i < 3; i++) {
 580     for (size_t j = 0; j < 3; j++)
 581       std::cout << std::setw(15) << rot[i * 3 + j];
 582     std::cout << std::endl;
 583   }
 584 }
 585
 586 //____________________________________________________________________
 587 //
 588 // EOF
 589 //