Non-implemented private copy constructor and assignment operator
[u/mrichter/AliRoot.git] / 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 //