1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
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: //
28 // tilt around rphi //
31 // The shifts are in cm and the tilts are in degrees. //
32 // The currently supported formats are: //
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) //
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}; //
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"); //
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(); //
64 // D.Miskowiec, November 2006 //
66 ///////////////////////////////////////////////////////////////////////////////
74 #include "TGeoManager.h"
75 #include "TGeoPhysicalNode.h"
76 #include "TClonesArray.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"
90 #include "AliTRDalignment.h"
92 void trdAlignmentFcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
94 ClassImp(AliTRDalignment)
96 //_____________________________________________________________________________
97 AliTRDalignment::AliTRDalignment()
108 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++) {
109 fSurveyX[i][j][k][l] = 0.0;
110 fSurveyY[i][j][k][l] = 0.0;
111 fSurveyZ[i][j][k][l] = 0.0;
112 fSurveyE[i][j][k][l] = 0.0;
115 // Initialize the nominal positions of the survey points
116 // in the local frame of supermodule (where y is the long side,
117 // z corresponds to the radius in lab, and x to the phi in lab).
118 // Four survey marks are on each z-side of the supermodule.
120 // ----o-----------o---- x |
125 // ---o-----o--- -------------->
128 // For the purpose of this explanation lets define the origin such that
129 // the supermodule occupies 0 < x < 77.9 cm. Then the coordinates (x,y)
137 double x[2] = {22.5,30.25}; // lab phi, or tracking-y
138 double y[2] = {353.0, -353.0}; // lab z; inc. 2 cm survey target offset
139 double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5}; // lab r, or better tracking-x
141 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
142 fSurveyX0[j][k][l] = -TMath::Power(-1,l) * x[k];
143 fSurveyY0[j][k][l] = y[j];
144 fSurveyZ0[j][k][l] = z[k];
149 //_____________________________________________________________________________
150 AliTRDalignment::AliTRDalignment(const AliTRDalignment& source)
152 ,fComment(source.fComment)
159 for (int i=0; i<18; i++) SetSm(i,source.fSm[i]);
160 for (int i=0; i<540; i++) SetCh(i,source.fCh[i]);
161 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++) {
162 fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
163 fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
164 fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
165 fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
167 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
168 fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
169 fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
170 fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
175 //_____________________________________________________________________________
176 AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source)
179 // assignment operator
182 if (this != &source) {
183 for (int i = 0; i < 18; i++) SetSm(i,source.fSm[i]);
184 for (int i = 0; i < 540; i++) SetCh(i,source.fCh[i]);
185 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++) {
186 fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
187 fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
188 fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
189 fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
191 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
192 fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
193 fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
194 fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
196 fComment = source.fComment;
203 //_____________________________________________________________________________
204 AliTRDalignment& AliTRDalignment::operator*=(double fac)
207 // multiplication operator
210 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] *= fac;
211 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] *= fac;
217 //_____________________________________________________________________________
218 AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source)
224 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] += source.fSm[i][j];
225 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] += source.fCh[i][j];
231 //_____________________________________________________________________________
232 AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source)
235 // subtraction operator
238 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) fSm[i][j] -= source.fSm[i][j];
239 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) fCh[i][j] -= source.fCh[i][j];
245 //_____________________________________________________________________________
246 Bool_t AliTRDalignment::operator==(const AliTRDalignment &source) const
249 // comparison operator
254 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) areEqual &= (fSm[i][j] == source.fSm[i][j]);
255 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) areEqual &= (fCh[i][j] == source.fCh[i][j]);
261 //_____________________________________________________________________________
262 void AliTRDalignment::SetSmZero()
265 // reset to zero supermodule data
268 memset(&fSm[0][0],0,sizeof(fSm));
272 //_____________________________________________________________________________
273 void AliTRDalignment::SetChZero()
276 // reset to zero chamber data
279 memset(&fCh[0][0],0,sizeof(fCh));
283 //_____________________________________________________________________________
284 void AliTRDalignment::SetSmRandom(double a[6])
287 // generate random gaussian supermodule data with sigmas a
292 for (int i = 0; i < 18; i++) {
293 fRan.Rannor(x[0],x[1]);
294 fRan.Rannor(x[2],x[3]);
295 fRan.Rannor(x[4],x[5]);
296 for (int j = 0; j < 6; j++) x[j] *= a[j];
303 //_____________________________________________________________________________
304 void AliTRDalignment::SetChRandom(double a[6])
307 // generate random gaussian chamber data with sigmas a
312 for (int i = 0; i < 540; i++) {
313 fRan.Rannor(x[0],x[1]);
314 fRan.Rannor(x[2],x[3]);
315 fRan.Rannor(x[4],x[5]);
316 for (int j = 0; j < 6; j++) x[j] *= a[j];
323 //_____________________________________________________________________________
324 void AliTRDalignment::SetSmFull()
327 // generate random gaussian supermodule data similar to the misalignment
328 // expected from the mechanical precision
336 a[3] = 0.4/1000.0 / TMath::Pi()*180.0; // phi
337 a[4] = 2.0/1000.0 / TMath::Pi()*180.0; // z
338 a[5] = 0.4/1000.0 / TMath::Pi()*180.0; // r
344 //_____________________________________________________________________________
345 void AliTRDalignment::SetChFull()
348 // generate random gaussian chamber data similar to the misalignment
349 // expected from the mechanical precision
357 a[3] = 1.0/1000.0 / TMath::Pi()*180.0; // phi
358 a[4] = 1.0/1000.0 / TMath::Pi()*180.0; // z
359 a[5] = 0.7/1000.0 / TMath::Pi()*180.0; // r
365 //_____________________________________________________________________________
366 void AliTRDalignment::SetSmResidual()
369 // generate random gaussian supermodule data similar to the misalignment
370 // remaining after full calibration
371 // I assume that it will be negligible
378 //_____________________________________________________________________________
379 void AliTRDalignment::SetChResidual()
382 // generate random gaussian chamber data similar to the misalignment
383 // remaining after full calibration
391 a[3] = 0.3/1000.0 / TMath::Pi()*180.0; // phi
392 a[4] = 0.3/1000.0 / TMath::Pi()*180.0; // z
393 a[5] = 0.1/1000.0 / TMath::Pi()*180.0; // r
399 //_____________________________________________________________________________
400 void AliTRDalignment::PrintSm(int i, FILE *fp) const
403 // print the supermodule data
406 fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
407 ,i,fSm[i][0],fSm[i][1],fSm[i][2],fSm[i][3],fSm[i][4],fSm[i][5]
412 //_____________________________________________________________________________
413 void AliTRDalignment::PrintCh(int i, FILE *fp) const
416 // print the chamber data
419 fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
420 ,i,fCh[i][0],fCh[i][1],fCh[i][2],fCh[i][3],fCh[i][4],fCh[i][5]
421 ,GetVoi(i),GetChName(i));
425 //_____________________________________________________________________________
426 void AliTRDalignment::ReadAscii(char *filename)
429 // read the alignment data from ascii file
432 double x[6]; // alignment data
433 int volid; // volume id
434 std::string syna; // symbolic name
435 int j; // dummy index
437 fstream fi(filename,fstream::in);
439 AliError(Form("cannot open input file %s",filename));
445 for (int i = 0; i < 18; i++) {
446 fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
447 if (j != i) AliError(Form("sm %d expected, %d found",i,j));
448 if (volid != 0) AliError(Form("sm %d volume id %d expected, %d found",i,0,volid));
449 std::string symnam = GetSmName(i);
450 if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data()));
456 for (int i = 0; i < 540; i++) {
457 fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
458 if (j != i) AliError(Form("ch %d expected, %d found",i,j));
459 if (volid != GetVoi(i)) AliError(Form("ch %d volume id %d expected, %d found",i,GetVoi(i),volid));
460 std::string symnam = GetChName(i);
461 if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data()));
469 //_____________________________________________________________________________
470 void AliTRDalignment::ReadCurrentGeo()
473 // use currently loaded geometry to determine misalignment by comparing
474 // original and misaligned matrix of the last node
475 // Now, original, does not mean "ideal". It is the matrix before the alignment.
476 // So, if alignment was applied more than once, the numbers extracted will
477 // represent just the last alignment. -- check this!
480 TGeoHMatrix *ideSm[18]; // ideal
481 TGeoHMatrix *misSm[18]; // misaligned
482 for (int i = 0; i < 18; i++) {
484 // read misaligned and original matrices
486 TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetSmName(i));
487 if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
489 TGeoPhysicalNode *node = pne->GetPhysicalNode();
490 if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
492 misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
493 ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
495 // calculate the local misalignment matrices as inverse misaligned times ideal
497 TGeoHMatrix mat(ideSm[i]->Inverse());
498 mat.Multiply(misSm[i]);
499 double *tra = mat.GetTranslation();
500 double *rot = mat.GetRotationMatrix();
505 if (TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) AliError("Failed to extract roll-pitch-yall angles!");
506 double raddeg = TMath::RadToDeg();
507 pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
508 pars[4] = raddeg * TMath::ASin(rot[2]);
509 pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
518 TGeoHMatrix *ideCh[540]; // ideal
519 TGeoHMatrix *misCh[540]; // misaligned
520 for (int i = 0; i < 540; i++) {
522 // read misaligned and original matrices
524 TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetChName(i));
525 if (!pne) AliError(Form("no such physical node entry: %s",GetChName(i)));
527 TGeoPhysicalNode *node = pne->GetPhysicalNode();
528 if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i)));
530 misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
531 ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
533 // calculate the local misalignment matrices as inverse misaligned times ideal
535 TGeoHMatrix mat(ideCh[i]->Inverse());
536 mat.Multiply(misCh[i]);
537 double *tra = mat.GetTranslation();
538 double *rot = mat.GetRotationMatrix();
543 if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
544 AliError("Failed to extract roll-pitch-yall angles!");
547 double raddeg = TMath::RadToDeg();
548 pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
549 pars[4] = raddeg * TMath::ASin(rot[2]);
550 pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
562 //_____________________________________________________________________________
563 void AliTRDalignment::ReadRoot(char *filename)
566 // read the alignment data from root file
569 TFile fi(filename,"READ");
572 TClonesArray *ar = (TClonesArray*) fi.Get("TRDAlignObjs");
576 else AliError(Form("cannot open input file %s",filename));
582 //_____________________________________________________________________________
583 void AliTRDalignment::ReadDB(char *filename)
586 // read the alignment data from database file
589 TFile fi(filename,"READ");
592 AliCDBEntry *e = (AliCDBEntry *) fi.Get("AliCDBEntry");
594 fComment.SetString(e->GetMetaData()->GetComment());
595 TClonesArray *ar = (TClonesArray *) e->GetObject();
599 else AliError(Form("cannot open input file %s",filename));
605 //_____________________________________________________________________________
606 void AliTRDalignment::ReadDB(char *db, char *path, int run
607 , int version, int subversion)
610 // read the alignment data from database
613 AliCDBManager *cdb = AliCDBManager::Instance();
614 AliCDBStorage *storLoc = cdb->GetStorage(db);
615 AliCDBEntry *e = storLoc->Get(path,run,version,subversion);
618 fComment.SetString(e->GetMetaData()->GetComment());
619 TClonesArray *ar = (TClonesArray *) e->GetObject();
624 //_____________________________________________________________________________
625 void AliTRDalignment::ReadSurveyReport(char *filename)
628 // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ,
629 // and fSurveyE. Store the survey info in the fComment.
630 // Each supermodule has 8 survey points. The point names look like
631 // TRD_sm08ah0 and have the following meaning.
633 // sm00..17 mean supermodule 0 through 17, following the phi.
634 // Supermodule 00 is between phi=0 and phi=20 degrees.
636 // a or c denotes the anticlockwise and clockwise end of the supermodule
637 // in z. Clockwise end is where z is negative and where the muon arm sits.
639 // l or h denote low radius and high radius holes
641 // 0 or 1 denote the hole at smaller and at larger phi, respectively.
644 // read the survey file
646 fstream in(filename,fstream::in);
648 AliError(Form("cannot open input file %s",filename));
652 // loop through the lines of the file until the beginning of data
654 TString title,date,subdetector,url,version,observations,system,units;
660 if (line.Contains("Title:")) title.ReadLine(in);
661 if (line.Contains("Date:")) date.ReadLine(in);
662 if (line.Contains("Subdetector:")) subdetector.ReadLine(in);
663 if (line.Contains("URL:")) url.ReadLine(in);
664 if (line.Contains("Version:")) version.ReadLine(in);
665 if (line.Contains("Observations:")) observations.ReadLine(in);
666 if (line.Contains("System:")) system.ReadLine(in);
667 if (line.Contains("Units:")) units.ReadLine(in);
668 if (line.Contains("Data:")) break;
671 // check what we found so far (watch out, they have \r at the end)
673 std::cout<<"title .........."<<title<<std::endl;
674 std::cout<<"date ..........."<<date<<std::endl;
675 std::cout<<"subdetector ...."<<subdetector<<std::endl;
676 std::cout<<"url ............"<<url<<std::endl;
677 std::cout<<"version ........"<<version<<std::endl;
678 std::cout<<"observations ..."<<observations<<std::endl;
679 std::cout<<"system ........."<<system<<std::endl;
680 std::cout<<"units .........."<<units<<std::endl;
682 if (!subdetector.Contains("TRD")) {
683 AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
686 double tocm = 0; // we want to have it in cm
687 if (units.Contains("mm")) tocm = 0.1;
688 else if (units.Contains("cm")) tocm = 1.0;
689 else if (units.Contains("m")) tocm = 100.0;
690 else if (units.Contains("pc")) tocm = 3.24078e-15;
692 AliError(Form("unexpected units: %s",units.Data()));
695 if (!system.Contains("ALICEPH")) {
696 AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
700 // scan the rest of the file which should contain list of surveyed points
701 // for every point, decode the point name and store the numbers in the right
702 // place in the arrays fSurveyX etc.
705 TString pna; // point name
707 double x,y,z,precision;
708 in >> pna >> x >> y >> z >> type >> precision;
709 if (in.fail()) break;
710 if (pna(0,6)!="TRD_sm") {
711 AliError(Form("unexpected point name: %s",pna.Data()));
714 int i = atoi(pna(6,0).Data()); // supermodule number
716 if (pna(8) == 'a') j=0; // anticlockwise, positive z
717 if (pna(8) == 'c') j=1; // clockwise, negative z
719 if (pna(9) == 'l') k=0; // low radius
720 if (pna(9) == 'h') k=1; // high radius
721 int l = atoi(pna(10,0).Data()); // phi within supermodule
722 if (i>=0 && i<18 && j>=0 && j<2 && k>=0 && k<2 && l>=0 && l<2) {
723 fSurveyX[i][j][k][l] = tocm*x;
724 fSurveyY[i][j][k][l] = tocm*y;
725 fSurveyZ[i][j][k][l] = tocm*z;
726 fSurveyE[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
727 std::cout << "decoded "<<pna<<" "
728 <<fSurveyX[i][j][k][l]<<" "
729 <<fSurveyY[i][j][k][l]<<" "
730 <<fSurveyZ[i][j][k][l]<<" "
731 <<fSurveyE[i][j][k][l]<<std::endl;
732 } else AliError(Form("cannot decode point name: %s",pna.Data()));
735 TString info = "Survey "+title+" "+date+" "+url+" "+version+" "+observations;
736 info.ReplaceAll("\r","");
737 fComment.SetString(info.Data());
741 //_____________________________________________________________________________
742 double AliTRDalignment::SurveyChi2(int i, double *a) {
745 // Compare the survey results to the ideal positions of the survey marks
746 // in the local frame of supermodule. When transforming, use the alignment
747 // parameters a[6]. Return chi-squared.
750 if (!IsGeoLoaded()) return 0;
751 printf("Survey of supermodule %d\n",i);
752 AliAlignObjParams al(GetSmName(i),0,a[0],a[1],a[2],a[3],a[4],a[5],0);
753 TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetSmName(i));
754 if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
755 TGeoPhysicalNode *node = pne->GetPhysicalNode();
756 if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
758 // al.ApplyToGeometry();
759 // node = pne->GetPhysicalNode(); // changed in the meantime
760 // TGeoHMatrix *ma = node->GetMatrix();
762 // a less destructive method (it does not modify geometry), gives the same result:
764 TGeoHMatrix *ma = new TGeoHMatrix();
765 al.GetLocalMatrix(*ma);
766 ma->MultiplyLeft(node->GetMatrix()); // global trafo, modified by a[]
769 printf(" sm z r phi x (lab phi) y (lab z) z (lab r) all in cm\n");
770 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
771 if (fSurveyE[i][j][k][l] == 0.0) continue; // no data for this survey point
772 double master[3] = {fSurveyX[i][j][k][l],fSurveyY[i][j][k][l],fSurveyZ[i][j][k][l]};
774 ma->MasterToLocal(master,local);
775 double dx = local[0]-fSurveyX0[j][k][l];
776 double dy = local[1]-fSurveyY0[j][k][l];
777 double dz = local[2]-fSurveyZ0[j][k][l];
778 chi2 += (dx*dx+dy*dy+dz*dz)/fSurveyE[i][j][k][l]/fSurveyE[i][j][k][l];
779 printf("local survey %3d %3d %3d %3d %12.3f %12.3f %12.3f\n",i,j,k,l,local[0],local[1],local[2]);
780 printf("local ideal %12.3f %12.3f %12.3f\n",fSurveyX0[j][k][l],
781 fSurveyY0[j][k][l],fSurveyZ0[j][k][l]);
782 printf("difference %12.3f %12.3f %12.3f\n",dx,dy,dz);
784 printf("chi2 = %.2f\n",chi2);
788 //_____________________________________________________________________________
789 void trdAlignmentFcn(int &npar, double *g, double &f, double *par, int iflag) {
792 // Standard function as needed by Minuit-like minimization procedures.
793 // For the set of parameters par calculates and returns chi-squared.
796 // smuggle a C++ object into a C function
797 AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit();
799 f = alignment->SurveyChi2(par);
802 if (g); // no warnings about unused stuff...
806 //_____________________________________________________________________________
807 void AliTRDalignment::SurveyToAlignment(int i,char *flag) {
810 // Find the supermodule alignment parameters needed to make the survey
811 // results coincide with the ideal positions of the survey marks.
812 // The string flag should look like "101000"; the six characters corresponds
813 // to the six alignment parameters and 0/1 mean that the parameter should
814 // be fixed/released in the fit.
816 if (strlen(flag)!=6) {
817 AliError(Form("unexpected flag: %s",flag));
821 printf("Finding alignment matrix for supermodule %d\n",i);
822 fIbuffer[0] = i; // store the sm number in the buffer so minuit can see it
825 gMinuit->SetObjectFit(this);
826 fitter.SetFCN(trdAlignmentFcn);
827 fitter.SetParameter(0,"dx",0,0.5,0,0);
828 fitter.SetParameter(1,"dy",0,0.5,0,0);
829 fitter.SetParameter(2,"dz",0,0.5,0,0);
830 fitter.SetParameter(3,"rx",0,0.1,0,0);
831 fitter.SetParameter(4,"ry",0,0.1,0,0);
832 fitter.SetParameter(5,"rz",0,0.1,0,0);
834 for (int j=0; j<6; j++) if (flag[j]=='0') fitter.FixParameter(j);
838 fitter.ExecuteCommand("SET PRINT", arglist, 1);
839 fitter.ExecuteCommand("SET ERR", arglist, 1);
841 //fitter.ExecuteCommand("SIMPLEX", arglist, 1);
842 fitter.ExecuteCommand("MINIMIZE", arglist, 1);
843 fitter.ExecuteCommand("CALL 3", arglist,0);
845 for (int j=0; j<6; j++) a[j] = fitter.GetParameter(j);
847 for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParameter(j));
849 for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParError(j));
854 //_____________________________________________________________________________
855 void AliTRDalignment::ReadAny(char *filename)
858 // read the alignment data from any kind of file
861 TString fist(filename);
862 if (fist.EndsWith(".txt")) ReadAscii(filename);
863 if (fist.EndsWith(".dat")) ReadAscii(filename);
864 if (fist.EndsWith(".root")) {
865 if (fist.Contains("Run")) ReadDB(filename);
866 else ReadRoot(filename);
871 //_____________________________________________________________________________
872 void AliTRDalignment::WriteAscii(char *filename) const
875 // store the alignment data on ascii file
878 FILE *fp = fopen(filename, "w");
880 AliError(Form("cannot open output file %s",filename));
891 //_____________________________________________________________________________
892 void AliTRDalignment::WriteRoot(char *filename)
895 // store the alignment data on root file
898 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
900 TFile fo(filename,"RECREATE");
903 fo.WriteObject(ar,"TRDAlignObjs","kSingleKey");
906 else AliError(Form("cannot open output file %s",filename));
912 //_____________________________________________________________________________
913 void AliTRDalignment::WriteDB(char *filename, int run0, int run1)
916 // dumping on a DB-like file
919 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
921 char *path = "TRD/Align/Data";
922 AliCDBId id(path,run0,run1);
923 AliCDBMetaData *md = new AliCDBMetaData();
924 md->SetResponsible("Dariusz Miskowiec");
925 md->SetComment(fComment.GetString().Data());
926 AliCDBEntry *e = new AliCDBEntry(ar, id, md);
927 TFile fi(filename,"RECREATE");
932 else AliError(Form("cannot open input file %s",filename));
942 //_____________________________________________________________________________
943 void AliTRDalignment::WriteDB(char *db, char *path, int run0, int run1)
946 // store the alignment data in database
949 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
951 AliCDBManager *cdb = AliCDBManager::Instance();
952 AliCDBStorage *storLoc = cdb->GetStorage(db);
953 AliCDBMetaData *md = new AliCDBMetaData();
954 md->SetResponsible("Dariusz Miskowiec");
955 md->SetComment(fComment.GetString().Data());
956 AliCDBId id(path,run0,run1);
957 storLoc->Put(ar,id,md);
963 //_____________________________________________________________________________
964 void AliTRDalignment::WriteGeo(char *filename)
967 // apply misalignment to current geometry and store the
968 // resulting geometry on a root file
971 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
974 gGeoManager->Export(filename);
978 //_____________________________________________________________________________
979 double AliTRDalignment::GetSmRMS(int xyz) const
987 for (int i = 0; i < 18; i++) {
989 s2 += fSm[i][xyz]*fSm[i][xyz];
991 double rms2 = s2/18.0 - s1*s1/18.0/18.0;
993 return rms2>0 ? sqrt(rms2) : 0.0;
997 //_____________________________________________________________________________
998 double AliTRDalignment::GetChRMS(int xyz) const
1006 for (int i = 0; i < 540; i++) {
1008 s2 += fCh[i][xyz]*fCh[i][xyz];
1010 double rms2 = s2/540.0 - s1*s1/540.0/540.0;
1012 return rms2>0 ? sqrt(rms2) : 0.0;
1016 //_____________________________________________________________________________
1017 void AliTRDalignment::PrintSmRMS() const
1023 printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f supermodule rms\n"
1024 ,GetSmRMS(0),GetSmRMS(1),GetSmRMS(2),GetSmRMS(3),GetSmRMS(4),GetSmRMS(5));
1028 //_____________________________________________________________________________
1029 void AliTRDalignment::PrintChRMS() const
1035 printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f chamber rms\n"
1036 ,GetChRMS(0),GetChRMS(1),GetChRMS(2),GetChRMS(3),GetChRMS(4),GetChRMS(5));
1040 //_____________________________________________________________________________
1041 void AliTRDalignment::ArToNumbers(TClonesArray *ar)
1044 // for each of the alignment objects in array ar extract the six local
1045 // alignment parameters; recognize by name to which supermodule or chamber
1046 // the alignment object pertains; set the respective fSm or fCh
1050 if (!IsGeoLoaded()) return;
1051 for (int i = 0; i < ar->GetEntries(); i++) {
1052 AliAlignObj *aao = (AliAlignObj *) ar->At(i);
1053 aao->ApplyToGeometry();
1060 //_____________________________________________________________________________
1061 void AliTRDalignment::NumbersToAr(TClonesArray *ar)
1064 // build array of AliAlignObj objects based on fSm and fCh data
1065 // at the same time, apply misalignment to the currently loaded geometry
1066 // it is important to apply misalignment of supermodules before creating
1067 // alignment objects for chambers
1070 if (!IsGeoLoaded()) return;
1071 TClonesArray &alobj = *ar;
1073 for (int i = 0; i < 18; i++) {
1074 new(alobj[nobj]) AliAlignObjParams(GetSmName(i)
1076 ,fSm[i][0],fSm[i][1],fSm[i][2]
1077 ,fSm[i][3],fSm[i][4],fSm[i][5]
1079 ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1083 for (int i = 0; i < 540; i++) {
1084 new(alobj[nobj]) AliAlignObjParams(GetChName(i)
1086 ,fCh[i][0],fCh[i][1],fCh[i][2]
1087 ,fCh[i][3],fCh[i][4],fCh[i][5]
1089 ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1092 AliInfo("current geometry modified");
1096 //_____________________________________________________________________________
1097 int AliTRDalignment::IsGeoLoaded()
1100 // check whether a geometry is loaded
1101 // issue a warning if geometry is not ideal
1105 if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) AliWarning("current geometry is not ideal");
1108 AliError("first load geometry by calling TGeoManager::Import(filename)");
1114 //_____________________________________________________________________________