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/OCDB", "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 ///////////////////////////////////////////////////////////////////////////////
73 #include "TGeoManager.h"
74 #include "TGeoPhysicalNode.h"
75 #include "TClonesArray.h"
81 #include "AliAlignObj.h"
82 #include "AliAlignObjParams.h"
83 #include "AliCDBManager.h"
84 #include "AliCDBStorage.h"
85 #include "AliCDBMetaData.h"
86 #include "AliCDBEntry.h"
87 #include "AliSurveyObj.h"
88 #include "AliSurveyPoint.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 fSurveyEX[i][j][k][l] = 0.0;
113 fSurveyEY[i][j][k][l] = 0.0;
114 fSurveyEZ[i][j][k][l] = 0.0;
117 // Initialize the nominal positions of the survey points
118 // in the local frame of supermodule (where y is the long side,
119 // z corresponds to the radius in lab, and x to the phi in lab).
120 // Four survey marks are on each z-side of the supermodule.
122 // ----o-----------o---- x |
127 // ---o-----o--- -------------->
130 // For the purpose of this explanation lets define the origin such that
131 // the supermodule occupies 0 < x < 77.9 cm. Then the coordinates (x,y)
139 double x[2] = {22.5,30.25}; // lab phi, or tracking-y
140 double y[2] = {353.0, -353.0}; // lab z; inc. 2 cm survey target offset
141 double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5}; // lab r, or better tracking-x
143 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
144 fSurveyX0[j][k][l] = -TMath::Power(-1,l) * x[k];
145 fSurveyY0[j][k][l] = y[j];
146 fSurveyZ0[j][k][l] = z[k];
149 for (int i=0; i<1000; i++) {
156 //_____________________________________________________________________________
157 AliTRDalignment::AliTRDalignment(const AliTRDalignment& source)
159 ,fComment(source.fComment)
166 for (int i=0; i<18; i++) SetSm(i,source.fSm[i]);
167 for (int i=0; i<540; i++) SetCh(i,source.fCh[i]);
168 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++) {
169 fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
170 fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
171 fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
172 fSurveyEX[i][j][k][l] = source.fSurveyEX[i][j][k][l];
173 fSurveyEY[i][j][k][l] = source.fSurveyEY[i][j][k][l];
174 fSurveyEZ[i][j][k][l] = source.fSurveyEZ[i][j][k][l];
176 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
177 fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
178 fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
179 fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
181 for (int i=0; i<1000; i++) {
188 //_____________________________________________________________________________
189 AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source)
192 // assignment operator
195 if (this != &source) {
196 for (int i = 0; i < 18; i++) SetSm(i,source.fSm[i]);
197 for (int i = 0; i < 540; i++) SetCh(i,source.fCh[i]);
198 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++) {
199 fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
200 fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
201 fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
202 fSurveyEX[i][j][k][l] = source.fSurveyEX[i][j][k][l];
203 fSurveyEY[i][j][k][l] = source.fSurveyEY[i][j][k][l];
204 fSurveyEZ[i][j][k][l] = source.fSurveyEZ[i][j][k][l];
206 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
207 fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
208 fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
209 fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
211 fComment = source.fComment;
218 //_____________________________________________________________________________
219 AliTRDalignment& AliTRDalignment::operator*=(double fac)
222 // multiplication operator
225 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] *= fac;
226 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] *= fac;
232 //_____________________________________________________________________________
233 AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source)
239 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] += source.fSm[i][j];
240 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] += source.fCh[i][j];
246 //_____________________________________________________________________________
247 AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source)
250 // subtraction operator
253 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) fSm[i][j] -= source.fSm[i][j];
254 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) fCh[i][j] -= source.fCh[i][j];
260 //_____________________________________________________________________________
261 Bool_t AliTRDalignment::operator==(const AliTRDalignment &source) const
264 // comparison operator
269 for (int i = 0; i < 18; i++) for (int j = 0; j < 6; j++) areEqual &= (fSm[i][j] == source.fSm[i][j]);
270 for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) areEqual &= (fCh[i][j] == source.fCh[i][j]);
276 //_____________________________________________________________________________
277 void AliTRDalignment::SetSmZero()
280 // reset to zero supermodule data
283 memset(&fSm[0][0],0,sizeof(fSm));
287 //_____________________________________________________________________________
288 void AliTRDalignment::SetChZero()
291 // reset to zero chamber data
294 memset(&fCh[0][0],0,sizeof(fCh));
298 //_____________________________________________________________________________
299 void AliTRDalignment::SetSmRandom(double a[6])
302 // generate random gaussian supermodule data with sigmas a
306 double xmax[6]={999, 0.6, 999, 999, 999, 999};
308 for (int i = 0; i < 18; i++) {
309 for (int j = 0; j < 6; j++) {
310 do {x[j] = fRan.Gaus(0,a[j]);} while (TMath::Abs(x[j]) > xmax[j]);
318 //_____________________________________________________________________________
319 void AliTRDalignment::SetChRandom(double a[6])
322 // generate random gaussian chamber data with sigmas a
327 for (int i = 0; i < 540; i++) {
328 fRan.Rannor(x[0],x[1]);
329 fRan.Rannor(x[2],x[3]);
330 fRan.Rannor(x[4],x[5]);
331 for (int j = 0; j < 6; j++) x[j] *= a[j];
338 //_____________________________________________________________________________
339 void AliTRDalignment::SetSmFull()
342 // generate random gaussian supermodule data similar to the misalignment
343 // expected from the mechanical precision
351 a[3] = 0.4/1000.0 / TMath::Pi()*180.0; // phi
352 a[4] = 2.0/1000.0 / TMath::Pi()*180.0; // z
353 a[5] = 0.4/1000.0 / TMath::Pi()*180.0; // r
359 //_____________________________________________________________________________
360 void AliTRDalignment::SetChFull()
363 // generate random gaussian chamber data similar to the misalignment
364 // expected from the mechanical precision
372 a[3] = 1.0/1000.0 / TMath::Pi()*180.0; // phi
373 a[4] = 1.0/1000.0 / TMath::Pi()*180.0; // z
374 a[5] = 0.7/1000.0 / TMath::Pi()*180.0; // r
380 //_____________________________________________________________________________
381 void AliTRDalignment::SetSmResidual()
384 // generate random gaussian supermodule data similar to the misalignment
385 // remaining after full calibration
386 // I assume that it will be negligible
393 //_____________________________________________________________________________
394 void AliTRDalignment::SetChResidual()
397 // generate random gaussian chamber data similar to the misalignment
398 // remaining after full calibration
406 a[3] = 0.3/1000.0 / TMath::Pi()*180.0; // phi
407 a[4] = 0.3/1000.0 / TMath::Pi()*180.0; // z
408 a[5] = 0.1/1000.0 / TMath::Pi()*180.0; // r
414 //_____________________________________________________________________________
415 void AliTRDalignment::PrintSm(int i, FILE * const fp) const
418 // print the supermodule data
421 fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
422 ,i,fSm[i][0],fSm[i][1],fSm[i][2],fSm[i][3],fSm[i][4],fSm[i][5]
427 //_____________________________________________________________________________
428 void AliTRDalignment::PrintCh(int i, FILE * const fp) const
431 // print the chamber data
434 fprintf(fp,"%4d %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f %6d %s\n"
435 ,i,fCh[i][0],fCh[i][1],fCh[i][2],fCh[i][3],fCh[i][4],fCh[i][5]
436 ,GetVoi(i),GetChName(i));
440 //_____________________________________________________________________________
441 void AliTRDalignment::ReadAscii(const char * const filename)
444 // read the alignment data from ascii file
447 double x[6]; // alignment data
448 int volid; // volume id
449 std::string syna; // symbolic name
450 int j; // dummy index
452 fstream fi(filename,fstream::in);
454 AliError(Form("cannot open input file %s",filename));
460 for (int i = 0; i < 18; i++) {
461 fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
462 if (j != i) AliError(Form("sm %d expected, %d found",i,j));
463 if (volid != 0) AliError(Form("sm %d volume id %d expected, %d found",i,0,volid));
464 std::string symnam = GetSmName(i);
465 if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data()));
471 for (int i = 0; i < 540; i++) {
472 fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
473 if (j != i) AliError(Form("ch %d expected, %d found",i,j));
474 if (volid != GetVoi(i)) AliError(Form("ch %d volume id %d expected, %d found",i,GetVoi(i),volid));
475 std::string symnam = GetChName(i);
476 if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data()));
484 //_____________________________________________________________________________
485 void AliTRDalignment::ReadCurrentGeo()
488 // use currently loaded geometry to determine misalignment by comparing
489 // original and misaligned matrix of the last node
490 // Now, original, does not mean "ideal". It is the matrix before the alignment.
491 // So, if alignment was applied more than once, the numbers extracted will
492 // represent just the last alignment. -- check this!
496 TGeoHMatrix *ideSm[18]; // ideal
497 TGeoHMatrix *misSm[18]; // misaligned
498 for (int i = 0; i < 18; i++) if ((pne = gGeoManager->GetAlignableEntry(GetSmName(i)))) {
500 // read misaligned and original matrices
502 TGeoPhysicalNode *node = pne->GetPhysicalNode();
503 if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i)));
505 misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
506 ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
508 // calculate the local misalignment matrices as inverse misaligned times ideal
510 TGeoHMatrix mat(ideSm[i]->Inverse());
511 mat.Multiply(misSm[i]);
512 double *tra = mat.GetTranslation();
513 double *rot = mat.GetRotationMatrix();
518 if (TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) AliError("Failed to extract roll-pitch-yall angles!");
519 double raddeg = TMath::RadToDeg();
520 pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
521 pars[4] = raddeg * TMath::ASin(rot[2]);
522 pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
531 TGeoHMatrix *ideCh[540]; // ideal
532 TGeoHMatrix *misCh[540]; // misaligned
533 for (int i = 0; i < 540; i++) if ((pne = gGeoManager->GetAlignableEntry(GetChName(i)))) {
535 // read misaligned and original matrices
537 TGeoPhysicalNode *node = pne->GetPhysicalNode();
538 if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i)));
540 misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
541 ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
543 // calculate the local misalignment matrices as inverse misaligned times ideal
545 TGeoHMatrix mat(ideCh[i]->Inverse());
546 mat.Multiply(misCh[i]);
547 double *tra = mat.GetTranslation();
548 double *rot = mat.GetRotationMatrix();
553 if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
554 AliError("Failed to extract roll-pitch-yall angles!");
557 double raddeg = TMath::RadToDeg();
558 pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
559 pars[4] = raddeg * TMath::ASin(rot[2]);
560 pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
572 //_____________________________________________________________________________
573 void AliTRDalignment::ReadRoot(const char * const filename)
576 // read the alignment data from root file
579 TFile fi(filename,"READ");
582 TClonesArray *ar = (TClonesArray*) fi.Get("TRDAlignObjs");
586 else AliError(Form("cannot open input file %s",filename));
592 //_____________________________________________________________________________
593 void AliTRDalignment::ReadDB(const char * const filename)
596 // read the alignment data from database file
599 TFile fi(filename,"READ");
602 AliCDBEntry *e = (AliCDBEntry *) fi.Get("AliCDBEntry");
604 fComment.SetString(e->GetMetaData()->GetComment());
605 TClonesArray *ar = (TClonesArray *) e->GetObject();
609 else AliError(Form("cannot open input file %s",filename));
615 //_____________________________________________________________________________
616 void AliTRDalignment::ReadDB(const char * const db, const char * const path,
617 int run, int version, int subversion)
620 // read the alignment data from database
623 AliCDBManager *cdb = AliCDBManager::Instance();
624 AliCDBStorage *storLoc = cdb->GetStorage(db);
625 AliCDBEntry *e = storLoc->Get(path,run,version,subversion);
628 fComment.SetString(e->GetMetaData()->GetComment());
629 TClonesArray *ar = (TClonesArray *) e->GetObject();
634 //_____________________________________________________________________________
635 Bool_t AliTRDalignment::DecodeSurveyPointName(TString pna, Int_t &sm, Int_t &iz,
636 Int_t &ir, Int_t &iphi) {
637 // decode the survey point name and extract the sm, z, r and phi indices
639 if (pna(0,6)!="TRD_sm") {
640 AliError(Form("unexpected point name: %s",pna.Data()));
643 sm = atoi(pna(6,2).Data()); // supermodule number
645 if (pna(8) == 'a') iz=0; // anticlockwise, positive z
646 if (pna(8) == 'c') iz=1; // clockwise, negative z
648 if (pna(9) == 'l') ir=0; // low radius
649 if (pna(9) == 'h') ir=1; // high radius
651 if (pna(10) == '0') iphi = 0; // low phi within supermodule
652 if (pna(10) == '1') iphi = 1; // high phi within supermodule
653 if (sm>=0 && sm<18 && iz>=0 && iz<2 && ir>=0 && ir<2 && iphi>=0 && iphi<2) return kTRUE;
654 AliError(Form("cannot decode point name: %s",pna.Data()));
658 //_____________________________________________________________________________
659 void AliTRDalignment::ReadSurveyReport(const char * const filename)
662 // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ,
663 // and fSurveyE. Store the survey info in the fComment.
664 // Each supermodule has 8 survey points. The point names look like
665 // TRD_sm08ah0 and have the following meaning.
667 // sm00..17 mean supermodule 0 through 17, following the phi.
668 // Supermodule 00 is between phi=0 and phi=20 degrees.
670 // a or c denotes the anticlockwise and clockwise end of the supermodule
671 // in z. Clockwise end is where z is negative and where the muon arm sits.
673 // l or h denote low radius and high radius holes
675 // 0 or 1 denote the hole at smaller and at larger phi, respectively.
678 // read the survey file
680 fstream in(filename,fstream::in);
682 AliError(Form("cannot open input file %s",filename));
686 // loop through the lines of the file until the beginning of data
688 TString title,date,subdetector,url,version,observations,system,units;
694 if (line.Contains("Title:")) title.ReadLine(in);
695 if (line.Contains("Date:")) date.ReadLine(in);
696 if (line.Contains("Subdetector:")) subdetector.ReadLine(in);
697 if (line.Contains("URL:")) url.ReadLine(in);
698 if (line.Contains("Version:")) version.ReadLine(in);
699 if (line.Contains("Observations:")) observations.ReadLine(in);
700 if (line.Contains("System:")) system.ReadLine(in);
701 if (line.Contains("Units:")) units.ReadLine(in);
702 if (line.Contains("Data:")) break;
705 // check what we found so far (watch out, they have \r at the end)
707 std::cout<<"title .........."<<title<<std::endl;
708 std::cout<<"date ..........."<<date<<std::endl;
709 std::cout<<"subdetector ...."<<subdetector<<std::endl;
710 std::cout<<"url ............"<<url<<std::endl;
711 std::cout<<"version ........"<<version<<std::endl;
712 std::cout<<"observations ..."<<observations<<std::endl;
713 std::cout<<"system ........."<<system<<std::endl;
714 std::cout<<"units .........."<<units<<std::endl;
716 if (!subdetector.Contains("TRD")) {
717 AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
720 double tocm = 0; // we want to have it in cm
721 if (units.Contains("mm")) tocm = 0.1;
722 else if (units.Contains("cm")) tocm = 1.0;
723 else if (units.Contains("m")) tocm = 100.0;
724 else if (units.Contains("pc")) tocm = 3.24078e-15;
726 AliError(Form("unexpected units: %s",units.Data()));
729 if (!system.Contains("ALICEPH")) {
730 AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
734 // scan the rest of the file which should contain list of surveyed points
735 // for every point, decode the point name and store the numbers in the right
736 // place in the arrays fSurveyX etc.
739 TString pna; // point name
741 double x,y,z,precision;
743 in >> pna >> x >> y >> z >> type >> target >> precision;
744 if (in.fail()) break;
746 if (DecodeSurveyPointName(pna,i,j,k,l)) {
747 fSurveyX[i][j][k][l] = tocm*x;
748 fSurveyY[i][j][k][l] = tocm*y;
749 fSurveyZ[i][j][k][l] = tocm*z;
750 fSurveyEX[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
751 fSurveyEY[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
752 fSurveyEZ[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
753 // if, at some point, separate precision numbers for x,y,z show up in the
754 // survey reports the function will fail here
755 printf("decoded %s %02d %d %d %d %8.2f %8.2f %8.2f %6.2f %6.2f %6.2f\n",
756 pna.Data(), i, j, k, l,
757 fSurveyX[i][j][k][l], fSurveyY[i][j][k][l], fSurveyZ[i][j][k][l],
758 fSurveyEX[i][j][k][l], fSurveyEY[i][j][k][l], fSurveyEZ[i][j][k][l]);
759 } else AliError(Form("cannot decode point name: %s",pna.Data()));
762 TString info = "Survey "+title+" "+date+" "+url+" "+version+" "+observations;
763 info.ReplaceAll("\r","");
764 fComment.SetString(info.Data());
768 //_____________________________________________________________________________
769 void AliTRDalignment::ReadSurveyReport(const AliSurveyObj * const so)
772 // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ,
773 // and fSurveyE. Store the survey info in the fComment.
774 // Each supermodule has 8 survey points. The point names look like
775 // TRD_sm08ah0 and have the following meaning.
777 // sm00..17 mean supermodule 0 through 17, following the phi.
778 // Supermodule 00 is between phi=0 and phi=20 degrees.
780 // a or c denotes the anticlockwise and clockwise end of the supermodule
781 // in z. Clockwise end is where z is negative and where the muon arm sits.
783 // l or h denote low radius and high radius holes
785 // 0 or 1 denote the hole at smaller and at larger phi, respectively.
788 // read and process the data from the survey object
790 Int_t size = so->GetEntries();
791 printf("-> %d\n", size);
793 TString title = so->GetReportTitle();
794 TString date = so->GetReportDate();
795 TString subdetector = so->GetDetector();
796 TString url = so->GetURL();
797 TString report = so->GetReportNumber();
798 TString version = so->GetReportVersion();
799 TString observations = so->GetObservations();
800 TString system = so->GetCoordSys();
801 TString units = so->GetUnits();
803 // check what we found so far (watch out, they have \r at the end)
805 std::cout<<"title .........."<<title<<std::endl;
806 std::cout<<"date ..........."<<date<<std::endl;
807 std::cout<<"subdetector ...."<<subdetector<<std::endl;
808 std::cout<<"url ............"<<url<<std::endl;
809 std::cout<<"version ........"<<version<<std::endl;
810 std::cout<<"observations ..."<<observations<<std::endl;
811 std::cout<<"system ........."<<system<<std::endl;
812 std::cout<<"units .........."<<units<<std::endl;
814 if (!subdetector.Contains("TRD")) {
815 AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
818 double tocm = 0; // we want to have it in cm
819 if (units.Contains("mm")) tocm = 0.1;
820 else if (units.Contains("cm")) tocm = 1.0;
821 else if (units.Contains("m")) tocm = 100.0;
822 else if (units.Contains("pc")) tocm = 3.24078e-15;
824 AliError(Form("unexpected units: %s",units.Data()));
827 if (!system.Contains("ALICEPH")) {
828 AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
832 // for every survey point, decode the point name and store the numbers in
833 // the right place in the arrays fSurveyX etc.
835 TObjArray *points = so->GetData();
836 for (int ip = 0; ip<points->GetEntries(); ++ip) {
837 AliSurveyPoint *po = (AliSurveyPoint *) points->At(ip);
838 TString pna = po->GetPointName();
840 if (DecodeSurveyPointName(pna,i,j,k,l)) {
841 fSurveyX[i][j][k][l] = tocm*po->GetX();
842 fSurveyY[i][j][k][l] = tocm*po->GetY();
843 fSurveyZ[i][j][k][l] = tocm*po->GetZ();
844 fSurveyEX[i][j][k][l] = po->GetPrecisionX()/10; // "precision" is supposed to be in mm
845 fSurveyEY[i][j][k][l] = po->GetPrecisionY()/10;
846 fSurveyEZ[i][j][k][l] = po->GetPrecisionZ()/10;
847 printf("decoded %s %02d %d %d %d %8.2f %8.2f %8.2f %6.2f %6.2f %6.2f\n",
848 pna.Data(), i, j, k, l,
849 fSurveyX[i][j][k][l], fSurveyY[i][j][k][l], fSurveyZ[i][j][k][l],
850 fSurveyEX[i][j][k][l], fSurveyEY[i][j][k][l], fSurveyEZ[i][j][k][l]);
851 } else AliError(Form("cannot decode point name: %s",pna.Data()));
854 TString info = "Survey "+title+" "+date+" "+url+" "+report+" "+version+" "+observations;
855 info.ReplaceAll("\r","");
856 fComment.SetString(info.Data());
859 //_____________________________________________________________________________
860 double AliTRDalignment::SurveyChi2(int i, const double * const a) {
863 // Compare the survey results to the ideal positions of the survey marks
864 // in the local frame of supermodule. When transforming, use the alignment
865 // parameters a[6]. Return chi-squared.
868 if (!IsGeoLoaded()) return 0;
869 printf("Survey of supermodule %d\n",i);
870 AliAlignObjParams al(GetSmName(i),0,a[0],a[1],a[2],a[3],a[4],a[5],0);
872 TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(GetSmName(i));
873 if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
874 TGeoPhysicalNode *node = pne->GetPhysicalNode();
876 AliWarning(Form("physical node entry %s has no physical node; making a new one",GetSmName(i)));
877 node = gGeoManager->MakeAlignablePN(pne);
880 // al.ApplyToGeometry();
881 // node = pne->GetPhysicalNode(); // changed in the meantime
882 // TGeoHMatrix *ma = node->GetMatrix();
884 // a less destructive method (it does not modify geometry), gives the same result:
886 TGeoHMatrix *ma = new TGeoHMatrix();
887 al.GetLocalMatrix(*ma);
888 ma->MultiplyLeft(node->GetMatrix()); // global trafo, modified by a[]
891 printf(" sm z r phi x (lab phi) y (lab z) z (lab r) all in cm\n");
892 for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
893 if (fSurveyEX[i][j][k][l] == 0.0
894 && fSurveyEY[i][j][k][l] == 0.0
895 && fSurveyEZ[i][j][k][l] == 0.0) continue; // no data for this survey point
896 double master[3] = {fSurveyX[i][j][k][l],fSurveyY[i][j][k][l],fSurveyZ[i][j][k][l]};
898 ma->MasterToLocal(master,local);
899 double dx = local[0]-fSurveyX0[j][k][l];
900 double dy = local[1]-fSurveyY0[j][k][l];
901 double dz = local[2]-fSurveyZ0[j][k][l];
902 chi2 += dx*dx/fSurveyEX[i][j][k][l]/fSurveyEX[i][j][k][l];
903 chi2 += dy*dy/fSurveyEY[i][j][k][l]/fSurveyEY[i][j][k][l];
904 chi2 += dz*dz/fSurveyEZ[i][j][k][l]/fSurveyEZ[i][j][k][l];
905 printf("local survey %3d %3d %3d %3d %12.3f %12.3f %12.3f\n",i,j,k,l,local[0],local[1],local[2]);
906 printf("local ideal %12.3f %12.3f %12.3f\n",fSurveyX0[j][k][l],
907 fSurveyY0[j][k][l],fSurveyZ0[j][k][l]);
908 printf("difference %12.3f %12.3f %12.3f\n",dx,dy,dz);
910 printf("chi2 = %.2f\n",chi2);
914 //_____________________________________________________________________________
915 void trdAlignmentFcn(int &npar, double *g, double &f, double *par, int iflag) {
918 // Standard function as needed by Minuit-like minimization procedures.
919 // For the set of parameters par calculates and returns chi-squared.
922 // smuggle a C++ object into a C function
923 AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit();
925 f = alignment->SurveyChi2(par);
928 if (g) {} // no warnings about unused stuff...
932 //_____________________________________________________________________________
933 void AliTRDalignment::SurveyToAlignment(int i, const char * const flag) {
936 // Find the supermodule alignment parameters needed to make the survey
937 // results coincide with the ideal positions of the survey marks.
938 // The string flag should look like "101000"; the six characters corresponds
939 // to the six alignment parameters and 0/1 mean that the parameter should
940 // be fixed/released in the fit.
942 if (strlen(flag)!=6) {
943 AliError(Form("unexpected flag: %s",flag));
947 printf("Finding alignment matrix for supermodule %d\n",i);
948 fIbuffer[0] = i; // store the sm number in the buffer so minuit can see it
951 gMinuit->SetObjectFit(this);
952 fitter.SetFCN(trdAlignmentFcn);
953 fitter.SetParameter(0,"dx",0,0.5,0,0);
954 fitter.SetParameter(1,"dy",0,0.5,0,0);
955 fitter.SetParameter(2,"dz",0,0.5,0,0);
956 fitter.SetParameter(3,"rx",0,0.1,0,0);
957 fitter.SetParameter(4,"ry",0,0.1,0,0);
958 fitter.SetParameter(5,"rz",0,0.1,0,0);
960 for (int j=0; j<6; j++) if (flag[j]=='0') fitter.FixParameter(j);
964 fitter.ExecuteCommand("SET PRINT", arglist, 1);
965 fitter.ExecuteCommand("SET ERR", arglist, 1);
967 //fitter.ExecuteCommand("SIMPLEX", arglist, 1);
968 fitter.ExecuteCommand("MINIMIZE", arglist, 1);
969 fitter.ExecuteCommand("CALL 3", arglist,0);
971 for (int j=0; j<6; j++) a[j] = fitter.GetParameter(j);
973 for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParameter(j));
975 for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParError(j));
980 //_____________________________________________________________________________
981 void AliTRDalignment::ReadAny(const char * const filename)
984 // read the alignment data from any kind of file
987 TString fist(filename);
988 if (fist.EndsWith(".txt")) ReadAscii(filename);
989 if (fist.EndsWith(".dat")) ReadAscii(filename);
990 if (fist.EndsWith(".root")) {
991 if (fist.Contains("Run")) ReadDB(filename);
992 else ReadRoot(filename);
997 //_____________________________________________________________________________
998 void AliTRDalignment::WriteAscii(const char * const filename) const
1001 // store the alignment data on ascii file
1004 FILE *fp = fopen(filename, "w");
1006 AliError(Form("cannot open output file %s",filename));
1017 //_____________________________________________________________________________
1018 void AliTRDalignment::WriteRoot(const char * const filename)
1021 // store the alignment data on root file
1024 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
1026 TFile fo(filename,"RECREATE");
1029 fo.WriteObject(ar,"TRDAlignObjs","kSingleKey");
1032 else AliError(Form("cannot open output file %s",filename));
1038 //_____________________________________________________________________________
1039 void AliTRDalignment::WriteDB(const char * const filename, int run0, int run1)
1042 // dumping on a DB-like file
1045 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
1047 const Char_t *path = "TRD/Align/Data";
1048 AliCDBId id(path,run0,run1);
1049 AliCDBMetaData *md = new AliCDBMetaData();
1050 md->SetResponsible("Dariusz Miskowiec");
1051 md->SetComment(fComment.GetString().Data());
1052 AliCDBEntry *e = new AliCDBEntry(ar, id, md);
1053 TFile fi(filename,"RECREATE");
1058 else AliError(Form("cannot open input file %s",filename));
1068 //_____________________________________________________________________________
1069 void AliTRDalignment::WriteDB(char * const db, const char * const path, int run0, int run1)
1072 // store the alignment data in database
1075 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
1077 AliCDBManager *cdb = AliCDBManager::Instance();
1078 AliCDBStorage *storLoc = cdb->GetStorage(db);
1079 AliCDBMetaData *md = new AliCDBMetaData();
1080 md->SetResponsible("Dariusz Miskowiec");
1081 md->SetComment(fComment.GetString().Data());
1082 AliCDBId id(path,run0,run1);
1083 storLoc->Put(ar,id,md);
1089 //_____________________________________________________________________________
1090 void AliTRDalignment::WriteGeo(char *filename)
1093 // apply misalignment to current geometry and store the
1094 // resulting geometry on a root file
1097 TClonesArray *ar = new TClonesArray("AliAlignObjParams",10000);
1100 gGeoManager->Export(filename);
1104 //_____________________________________________________________________________
1105 double AliTRDalignment::GetSmRMS(int xyz) const
1113 for (int i = 0; i < 18; i++) {
1115 s2 += fSm[i][xyz]*fSm[i][xyz];
1117 double rms2 = s2/18.0 - s1*s1/18.0/18.0;
1119 return rms2>0 ? sqrt(rms2) : 0.0;
1123 //_____________________________________________________________________________
1124 double AliTRDalignment::GetChRMS(int xyz) const
1132 for (int i = 0; i < 540; i++) {
1134 s2 += fCh[i][xyz]*fCh[i][xyz];
1136 double rms2 = s2/540.0 - s1*s1/540.0/540.0;
1138 return rms2>0 ? sqrt(rms2) : 0.0;
1142 //_____________________________________________________________________________
1143 void AliTRDalignment::PrintSmRMS() const
1149 printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f supermodule rms\n"
1150 ,GetSmRMS(0),GetSmRMS(1),GetSmRMS(2),GetSmRMS(3),GetSmRMS(4),GetSmRMS(5));
1154 //_____________________________________________________________________________
1155 void AliTRDalignment::PrintChRMS() const
1161 printf(" %11.4f %11.4f %11.4f %11.5f %11.5f %11.5f chamber rms\n"
1162 ,GetChRMS(0),GetChRMS(1),GetChRMS(2),GetChRMS(3),GetChRMS(4),GetChRMS(5));
1166 //_____________________________________________________________________________
1167 void AliTRDalignment::ArToNumbers(TClonesArray * const ar)
1170 // for each of the alignment objects in array ar extract the six local
1171 // alignment parameters; recognize by name to which supermodule or chamber
1172 // the alignment object pertains; set the respective fSm or fCh
1176 if (!IsGeoLoaded()) return;
1177 for (int i = 0; i < ar->GetEntries(); i++) {
1178 AliAlignObj *aao = (AliAlignObj *) ar->At(i);
1179 aao->ApplyToGeometry();
1186 //_____________________________________________________________________________
1187 void AliTRDalignment::NumbersToAr(TClonesArray * const ar)
1190 // build array of AliAlignObj objects based on fSm and fCh data
1191 // at the same time, apply misalignment to the currently loaded geometry
1192 // it is important to apply misalignment of supermodules before creating
1193 // alignment objects for chambers
1196 if (!IsGeoLoaded()) return;
1197 TClonesArray &alobj = *ar;
1199 for (int i = 0; i < 18; i++) {
1200 new(alobj[nobj]) AliAlignObjParams(GetSmName(i)
1202 ,fSm[i][0],fSm[i][1],fSm[i][2]
1203 ,fSm[i][3],fSm[i][4],fSm[i][5]
1205 ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1209 for (int i = 0; i < 540; i++) {
1210 if (gGeoManager->GetAlignableEntry(GetChName(i))) {
1211 new(alobj[nobj]) AliAlignObjParams(GetChName(i)
1213 ,fCh[i][0],fCh[i][1],fCh[i][2]
1214 ,fCh[i][3],fCh[i][4],fCh[i][5]
1216 ((AliAlignObj *) alobj[nobj])->ApplyToGeometry();
1220 AliInfo("current geometry modified");
1224 //_____________________________________________________________________________
1225 int AliTRDalignment::IsGeoLoaded()
1228 // check whether a geometry is loaded
1229 // issue a warning if geometry is not ideal
1233 if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) AliWarning("current geometry is not ideal");
1236 AliError("first load geometry by calling TGeoManager::Import(filename)");
1242 //_____________________________________________________________________________