10 #include <TGraph2DErrors.h> 16 #include <TRotation.h> 17 #include <TGeoMatrix.h> 18 #include <TGeoManager.h> 19 #include <TGeoPhysicalNode.h> 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.;
40 TVector3& error)
const 49 if (unit == 0)
return kFALSE;
52 if (!obj)
return kFALSE;
55 point.SetXYZ(unit * p->
GetX(),
89 TVector3 ab(b-a), bc(c-a);
93 TVector3 nn(ab.Cross(bc));
94 if (nn.Mag() < 1e-8) {
95 Info(
"CalculatePlane",
"Normal vector is null vector");
113 TVector3 n(nn.Unit());
129 TVector3 orig(md - depth * n);
137 TVector3 uab(ab.Unit());
138 TVector3 ubc(bc.Unit());
140 for (
size_t i = 0; i < 3; i++) {
141 rot[i * 3 + 0] = ubc[i];
142 rot[i * 3 + 1] = uab[i];
145 rot[i * 3 + 2] = n[i];
176 Int_t nPoints = points.GetEntries();
178 AliError(Form(
"Cannot fit a plane equation to less than 4 survey points, " 179 "got only %d", nPoints));
185 for (
int i = 0; i < nPoints; i++) {
186 TVector3*
p =
static_cast<TVector3*
>(points.At(i));
187 TVector3* e =
static_cast<TVector3*
>(errors.At(i));
189 if (!p || !e)
continue;
191 g.SetPoint(i, p->X(), p->Y(), p->Z());
192 g.SetPointError(i, e->X(), e->Y(), e->Z());
198 AliError(Form(
"Only got %d survey points - no good for plane fit",
212 TF2
f(
"plane",
"-[0]*x-[1]*y-[2]",
213 g.GetXmin(), g.GetXmax(), g.GetYmin(), g.GetYmax());
217 TVector3 nv(
f.GetParameter(0),
f.GetParameter(1), 1);
218 TVector3 n(nv.Unit());
219 Double_t
p = -
f.GetParameter(2);
222 TVector3 a(1, 0,
f.Eval(1, 0)-
p);
223 TVector3
b(0, -1,
f.Eval(0, -1)-
p);
224 TVector3 ua(a.Unit());
225 TVector3 ub(
b.Unit());
231 for (
size_t i = 0; i < 3; i++) {
232 rot[i * 3 + 0] = ua[i];
233 rot[i * 3 + 1] = ub[i];
234 rot[i * 3 + 2] = n[i];
250 const Double_t*
trans,
251 TGeoHMatrix& delta)
const 270 TGeoMatrix* global =
gGeoManager->GetCurrentMatrix();
272 PrintRotation(Form(
"%s rot:", global->GetName()),global->GetRotationMatrix());
273 PrintVector(Form(
"%s trans:", global->GetName()),global->GetTranslation());
276 return MakeDelta(global, rot, trans, delta);
283 const Double_t*
trans,
284 TGeoHMatrix& delta)
const 299 TGeoHMatrix* geoM =
new TGeoHMatrix;
300 geoM->SetTranslation(trans);
301 geoM->SetRotation(rot);
307 delta = global->Inverse();
310 delta.MultiplyLeft(geoM);
317 Double_t getFMD1Offset()
319 static Double_t off = 0;
323 if (off != 0)
return off;
325 const char* lidN =
"FMD1_lid_mat0";
326 TGeoMatrix* lidM =
static_cast<TGeoMatrix*
>(
gGeoManager->GetListOfMatrices()
329 Error(
"getFMD1Offset",
"Couldn't find FMD1 lid transformation %s", lidN);
333 const Double_t* lidT = lidM->GetTranslation();
334 Double_t lidZ = lidT[2];
358 TVector3 icb, ict, ocb, oct, eicb, eict, eocb, eoct;
360 if (!
GetPoint(
"V0L_ICB", icb, eicb)) missing++;
361 if (!
GetPoint(
"V0L_ICT", ict, eict)) missing++;
362 if (!
GetPoint(
"V0L_OCB", ocb, eocb)) missing++;
363 if (!
GetPoint(
"V0L_OCT", oct, eoct)) missing++;
367 AliWarning(Form(
"Only got %d survey points - no good for FMD1 plane",
374 points.Add(&icb); errors.Add(&eicb);
375 points.Add(&ict); errors.Add(&eict);
376 points.Add(&oct); errors.Add(&eoct);
377 points.Add(&ocb); errors.Add(&eocb);
379 Bool_t ret =
FitPlane(points, errors, 0, trans, rot);
381 Warning(
"GetFMD1Plane",
"fit to plane failed");
383 for (Int_t i = 0; i < 4; i++) {
384 TVector3* v =
static_cast<TVector3*
>(points.At(i));
385 TVector3* e =
static_cast<TVector3*
>(errors.At(i));
386 Info(
"GetFMD1Plane",
"p%d=(%8f,%8f,%8f)+/-(%8f,%8f,%8f)",
387 i, v->X(), v->Y(), v->Z(), e->X(), e->Y(), e->Z());
390 Double_t off = getFMD1Offset();
391 Info(
"GetFMD1Plane",
"Lid offset is %f", off);
427 Double_t rot[9],
trans[3];
433 Double_t gRot[9], gTrans[3];
434 TVector3 ocb(-127, -220, 324.67);
435 TVector3 oct(-127, +220, 324.67);
436 TVector3 icb(+127, -220, 324.67);
437 TVector3 ict(+127, +220, 324.67);
439 Warning(
"DoFMD1",
"Failed to make reference plane");
444 TGeoRotation ggRot; ggRot.SetMatrix(gRot);
445 TGeoCombiTrans global(gTrans[0], gTrans[1], gTrans[2], &ggRot);
447 Double_t off = getFMD1Offset();
448 Info(
"DoFMD1",
"Lid offset is %f", off);
450 TGeoTranslation global(0,0,324.670-off);
478 const char*
names[] = {
"FMD2_ITOP",
"FMD2_OTOP",
479 "FMD2_IBOTM",
"FMD2_OBOTM",
480 "FMD2_IBOT",
"FMD2_OBOT",
482 const char** name =
names;
498 Warning(
"GetFMD2plane",
"Setting error on %d, %s to 0.4", i, *name);
499 e.SetXYZ(0.4, 0.4, 0.4);
501 points.Add(
new TVector3(p));
502 errors.Add(
new TVector3(e));
506 if (points.GetEntries() < 4) {
507 AliWarning(Form(
"Only got %d survey points - no good for FMD2 plane",
508 points.GetEntries()));
512 return FitPlane(points, errors, 0, trans, rot);
515 #define M(I,J) rot[(J-1) * 3 + (I-1)] 545 Double_t rot[9],
trans[3];
551 for (
int i = 0; i < 3; i++) {
552 for (
int j = 0; j < 3; j++) {
553 rot[i*3+j] = (i == j ? 1 : 0);
557 trans[0] = trans[1] = 0;
563 if (!
MakeDelta(
"/ALIC_1/F2MT_2/FMD2_support_0/FMD2_back_cover_2",
602 const char** file = files;
604 if ((*file)[0] ==
'\0') {
605 Warning(
"Run",
"no file specified");
610 Warning(
"Run",
"Failed to load %s", *file);
615 Int_t d = Int_t(sDet[sDet.Length()-1] -
'0');
616 Info(
"Run",
"Making alignment for %s (%d)", sDet.Data(), d);
619 case 1: ret =
DoFMD1();
break;
620 case 2: ret =
DoFMD2();
break;
622 Warning(
"Run",
"Do not know how to deal with %s", sDet.Data());
626 Warning(
"Run",
"Calculation for %s failed", sDet.Data());
643 id,0,0,0,0,0,0,kTRUE);
645 AliError(Form(
"Failed to create alignment object for %s", path.Data()));
649 AliError(Form(
"Failed to set local transforms on %s", path.Data()));
671 for (
int d = 1; d <= 3; d++) {
672 const char sides[] = {
'T',
'B', 0 };
673 const char* side = sides;
675 TString
path = TString::Format(
"FMD/FMD%d_%c", d, *side);
678 const char halves[] = {
'I', d == 1 ?
'\0' :
'O', 0 };
679 const char* half = halves;
681 int nsec = *half ==
'I' ? 10 : 20;
682 int start = *side ==
'T' ? 0 : nsec/2;
683 int end = *side ==
'T' ? nsec/2 : nsec;
684 for (
int s=start; s < end; s++) {
685 path = TString::Format(
"FMD/FMD%d_%c/FMD%c_%02d",
710 Int_t n = array.GetEntriesFast();
736 Double_t va[] = { v.X(), v.Y(), v.Z() };
751 << std::setw(15) << v[0]
752 << std::setw(15) << v[1]
753 << std::setw(15) << v[2]
770 std::cout << text << std::endl;
771 for (
size_t i = 0; i < 3; i++) {
772 for (
size_t j = 0; j < 3; j++)
773 std::cout << std::setw(15) << rot[i * 3 + j];
774 std::cout << std::endl;
Bool_t MakeDelta(const TGeoMatrix *global, const Double_t *rot, const Double_t *trans, TGeoHMatrix &delta) const
Float_t GetPrecisionY() const
AliAlignObjParams * CreateDefaultAlignObj(const TString &path, Int_t id=0)
Geometry mananger for the FMD.
virtual void InitTransformations(Bool_t force=kFALSE)
AliSurveyObj * fSurveyObj
Double_t GetUnitFactor() const
AliAlignObjParams * FindAlignObj(const TString &path) const
#define AliWarning(message)
bool trans(const AliFMDIndex &x, const AliFMDIndex &y, const AliFMDIndex &z)
Bool_t LoadSurveyFromLocalFile(const char *filename)
Singleton object of FMD geometry descriptions and parameters. This class is a singleton that handles ...
Bool_t GetFMD1Plane(Double_t *rot, Double_t *trans) const
Float_t GetPrecisionZ() const
TGeoManager * gGeoManager
static void PrintRotation(const char *text, const Double_t *rot)
TObjArray * fSurveyPoints
Bool_t GetFMD2Plane(Double_t *rot, Double_t *trans) const
Bool_t FillDefaultAlignObjs()
Bool_t GetPoint(const char *name, TVector3 &p, TVector3 &e) const
const char * GetSymName() const
TString GetDetector() const
Float_t GetPrecisionX() const
#define AliError(message)
virtual Bool_t SetLocalPars(Double_t x, Double_t y, Double_t z, Double_t psi, Double_t theta, Double_t phi)
Bool_t FitPlane(const TObjArray &points, const TObjArray &errors, Double_t depth, Double_t *trans, Double_t *rot) const
static AliFMDGeometry * Instance()
TClonesArray * fAlignObjArray
static void PrintVector(const char *text, const Double_t *v)
Bool_t CalculatePlane(const TVector3 &a, const TVector3 &b, const TVector3 &c, Double_t depth, Double_t *trans, Double_t *rot) const
TClonesArray * GetAlignObjArray() const