Conversion of survey data into alignable objects implemented
[u/mrichter/AliRoot.git] / TRD / AliTRDalignment.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 /* $Id$ */
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
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:                   //
25 // shift in rphi                                                             //
26 // shift in z                                                                //
27 // shift in r                                                                //
28 // tilt around rphi                                                          //
29 // tilt around z                                                             //
30 // tilt around r                                                             //
31 // The shifts are in cm and the tilts are in degrees.                        //
32 // The currently supported formats are:                                      //
33 // - ascii                                                                   //
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)                           //
38 //                                                                           //
39 // D.Miskowiec, November 2006                                                //
40 //                                                                           //
41 ///////////////////////////////////////////////////////////////////////////////
42
43 #include <iostream>
44 #include <fstream>
45 #include <string>
46
47 #include "TMath.h"
48 #include "TFile.h"
49 #include "TGeoManager.h"
50 #include "TGeoPhysicalNode.h"
51 #include "TClonesArray.h"
52 #include "TString.h"
53 #include "TFitter.h"
54 #include "TMinuit.h"
55
56 #include "AliLog.h"
57 #include "AliAlignObj.h"
58 #include "AliAlignObjAngles.h"
59 #include "AliCDBManager.h"
60 #include "AliCDBStorage.h"
61 #include "AliCDBMetaData.h"
62 #include "AliCDBEntry.h"
63 #include "AliCDBId.h"
64
65 #include "AliTRDalignment.h"
66
67 void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
68
69 ClassImp(AliTRDalignment)
70
71 //_____________________________________________________________________________
72 AliTRDalignment::AliTRDalignment() 
73   :TObject()
74   ,fComment()
75   ,fRan(0)
76 {
77   //
78   // constructor
79   //
80
81   SetZero();
82
83   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++) {
84     fSurveyX[i][j][k][l] = 0.0;
85     fSurveyY[i][j][k][l] = 0.0;
86     fSurveyZ[i][j][k][l] = 0.0;
87     fSurveyE[i][j][k][l] = 0.0;
88   }
89
90   // Initialize the nominal positions of the survey points 
91   // in the local frame of supermodule (where y is the long side, 
92   // z corresponds to the radius in lab, and x to the phi in lab).
93   // Four survey marks are on each z-side of the supermodule. 
94   //               A           B
95   //           ----o-----------o----        x |
96   //           \                   /          |
97   //            \                 /           |
98   //             \               /            |
99   //              \             /             |
100   //               ---o-----o---              -------------->
101   //                  C     D                              y
102   // 
103   // For the purpose of this explanation lets define the origin such that 
104   // the supermodule occupies 0 < x < 77.9 cm. Then the coordinates (x,y) 
105   // are (in cm) 
106   // A (76.2,-30.25)
107   // B (76.2,+30.25)
108   // C ( 2.2,-22.5 )
109   // D ( 2.2,+22.5 )
110   // 
111
112   double x[2] = {22.5,30.25};                   // lab phi, or tracking-y
113   double y[2] = {353.0, -353.0};                // lab z; inc. 2 cm survey target offset
114   double z[2] = {-(77.9/2.0-2.0),77.9/2.0-1.5}; // lab r, or better tracking-x
115
116   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
117     fSurveyX0[j][k][l] = -pow(-1,l) * x[k];
118     fSurveyY0[j][k][l] = y[j];
119     fSurveyZ0[j][k][l] = z[k];
120   }
121
122 }
123
124 //_____________________________________________________________________________
125 AliTRDalignment::AliTRDalignment(const AliTRDalignment& source) 
126   :TObject(source)
127   ,fComment(source.fComment)
128   ,fRan(source.fRan)
129 {
130   //
131   // copy constructor
132   //
133
134   for (int i=0; i<18; i++) SetSm(i,source.fSm[i]);
135   for (int i=0; i<540; i++) SetCh(i,source.fCh[i]);
136   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++) {
137     fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
138     fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
139     fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
140     fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
141   }
142   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
143     fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
144     fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
145     fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
146   }
147
148 }
149
150 //_____________________________________________________________________________
151 AliTRDalignment& AliTRDalignment::operator=(const AliTRDalignment &source) 
152 {
153   //
154   // assignment operator
155   //
156
157   if (this != &source) {
158     for (int i = 0; i <  18; i++) SetSm(i,source.fSm[i]);
159     for (int i = 0; i < 540; i++) SetCh(i,source.fCh[i]);
160     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++) {
161       fSurveyX[i][j][k][l] = source.fSurveyX[i][j][k][l];
162       fSurveyY[i][j][k][l] = source.fSurveyY[i][j][k][l];
163       fSurveyZ[i][j][k][l] = source.fSurveyZ[i][j][k][l];
164       fSurveyE[i][j][k][l] = source.fSurveyE[i][j][k][l];
165     }
166     for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
167       fSurveyX0[j][k][l] = source.fSurveyX0[j][k][l];
168       fSurveyY0[j][k][l] = source.fSurveyY0[j][k][l];
169       fSurveyZ0[j][k][l] = source.fSurveyZ0[j][k][l];
170     }
171     fComment = source.fComment;
172   }
173
174   return *this;
175
176 }
177
178 //_____________________________________________________________________________
179 AliTRDalignment& AliTRDalignment::operator*=(double fac) 
180 {
181   //
182   // multiplication operator
183   //
184
185   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] *= fac;
186   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] *= fac;
187
188   return *this;
189
190 }
191
192 //_____________________________________________________________________________
193 AliTRDalignment& AliTRDalignment::operator+=(const AliTRDalignment &source) 
194 {
195   //
196   // addition operator
197   //
198
199   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) this->fSm[i][j] += source.fSm[i][j];
200   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) this->fCh[i][j] += source.fCh[i][j];
201
202   return *this;
203
204 }
205
206 //_____________________________________________________________________________
207 AliTRDalignment& AliTRDalignment::operator-=(const AliTRDalignment &source) 
208 {
209   //
210   // subtraction operator
211   //
212
213   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) fSm[i][j] -= source.fSm[i][j];
214   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) fCh[i][j] -= source.fCh[i][j];
215
216   return *this;
217
218 }
219
220 //_____________________________________________________________________________
221 Bool_t AliTRDalignment::operator==(const AliTRDalignment &source) const
222 {
223   //
224   // comparison operator
225   //
226
227   Bool_t areEqual = 1;
228
229   for (int i = 0; i <  18; i++) for (int j = 0; j < 6; j++) areEqual &= (fSm[i][j] == source.fSm[i][j]);
230   for (int i = 0; i < 540; i++) for (int j = 0; j < 6; j++) areEqual &= (fCh[i][j] == source.fCh[i][j]);
231
232   return areEqual;
233
234 }
235
236 //_____________________________________________________________________________
237 void AliTRDalignment::SetSmZero() 
238 {
239   //
240   // reset to zero supermodule data
241   //
242
243   memset(&fSm[0][0],0,sizeof(fSm));
244
245 }
246
247 //_____________________________________________________________________________
248 void AliTRDalignment::SetChZero() 
249 {
250   //
251   // reset to zero chamber data
252   //
253
254   memset(&fCh[0][0],0,sizeof(fCh));
255
256 }
257
258 //_____________________________________________________________________________
259 void AliTRDalignment::SetSmRandom(double a[6]) 
260 {
261   //
262   // generate random gaussian supermodule data with sigmas a
263   //
264
265   double x[6];
266
267   for (int i = 0; i < 18; i++) {
268     fRan.Rannor(x[0],x[1]);
269     fRan.Rannor(x[2],x[3]);
270     fRan.Rannor(x[4],x[5]);
271     for (int j = 0; j < 6; j++) x[j] *= a[j];
272     SetSm(i,x);
273     //PrintSm(i);
274   }
275
276 }
277
278 //_____________________________________________________________________________
279 void AliTRDalignment::SetChRandom(double a[6]) 
280 {
281   //
282   // generate random gaussian chamber data with sigmas a
283   //
284
285   double x[6];
286
287   for (int i = 0; i < 540; i++) {
288     fRan.Rannor(x[0],x[1]);
289     fRan.Rannor(x[2],x[3]);
290     fRan.Rannor(x[4],x[5]);
291     for (int j = 0; j < 6; j++) x[j] *= a[j];
292     SetCh(i,x);
293     //PrintCh(i);
294   }
295
296 }
297
298 //_____________________________________________________________________________
299 void AliTRDalignment::SetSmFull() 
300 {
301   //
302   // generate random gaussian supermodule data similar to the misalignment 
303   // expected from the mechanical precision 
304   //
305
306   double a[6];
307
308   a[0] = 0.3; // phi
309   a[1] = 0.3; // z
310   a[2] = 0.3; // r
311   a[3] = 0.4/1000.0 / TMath::Pi()*180.0; // phi
312   a[4] = 2.0/1000.0 / TMath::Pi()*180.0; // z
313   a[5] = 0.4/1000.0 / TMath::Pi()*180.0; // r
314
315   SetSmRandom(a);
316
317 }
318
319 //_____________________________________________________________________________
320 void AliTRDalignment::SetChFull() 
321 {
322   //
323   // generate random gaussian chamber data similar to the misalignment 
324   // expected from the mechanical precision 
325   //
326
327   double a[6];
328
329   a[0] = 0.1; // phi
330   a[1] = 0.1; // z
331   a[2] = 0.1; // r
332   a[3] = 1.0/1000.0 / TMath::Pi()*180.0; // phi
333   a[4] = 1.0/1000.0 / TMath::Pi()*180.0; // z
334   a[5] = 0.7/1000.0 / TMath::Pi()*180.0; // r
335
336   SetChRandom(a);
337
338 }
339
340 //_____________________________________________________________________________
341 void AliTRDalignment::SetSmResidual() 
342 {
343   //
344   // generate random gaussian supermodule data similar to the misalignment 
345   // remaining after full calibration
346   // I assume that it will be negligible
347   //
348
349   SetSmZero();
350
351 }
352
353 //_____________________________________________________________________________
354 void AliTRDalignment::SetChResidual() 
355 {
356   //
357   // generate random gaussian chamber data similar to the misalignment 
358   // remaining after full calibration
359   //
360
361   double a[6];
362
363   a[0] = 0.002; // phi
364   a[1] = 0.003; // z
365   a[2] = 0.007; // r
366   a[3] = 0.3/1000.0 / TMath::Pi()*180.0; // phi
367   a[4] = 0.3/1000.0 / TMath::Pi()*180.0; // z
368   a[5] = 0.1/1000.0 / TMath::Pi()*180.0; // r
369
370   SetChRandom(a);
371
372 }
373
374 //_____________________________________________________________________________
375 void AliTRDalignment::PrintSm(int i, FILE *fp) const 
376 {
377   //
378   // print the supermodule data
379   //
380
381   fprintf(fp,"%4d   %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f   %6d  %s\n"
382          ,i,fSm[i][0],fSm[i][1],fSm[i][2],fSm[i][3],fSm[i][4],fSm[i][5]
383          ,0,GetSmName(i));
384
385 }
386
387 //_____________________________________________________________________________
388 void AliTRDalignment::PrintCh(int i, FILE *fp) const 
389 {
390   //
391   // print the chamber data
392   //
393
394   fprintf(fp,"%4d   %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f   %6d  %s\n"
395          ,i,fCh[i][0],fCh[i][1],fCh[i][2],fCh[i][3],fCh[i][4],fCh[i][5]
396          ,GetVoi(i),GetChName(i));
397
398 }
399
400 //_____________________________________________________________________________
401 void AliTRDalignment::ReadAscii(char *filename) 
402 {
403   //
404   // read the alignment data from ascii file
405   //
406
407   double x[6];      // alignment data
408   int volid;        // volume id
409   std::string syna; // symbolic name
410   int j;            // dummy index
411
412   fstream fi(filename,fstream::in);
413   if (!fi) {
414     AliError(Form("cannot open input file %s",filename));
415     return;
416   }
417
418   // supermodules
419
420   for (int i = 0; i < 18; i++) {
421     fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
422     if (j != i) AliError(Form("sm %d expected, %d found",i,j));
423     if (volid != 0) AliError(Form("sm %d volume id %d expected, %d found",i,0,volid));
424     std::string symnam = GetSmName(i);
425     if (syna != symnam) AliError(Form("sm %d name %s expected, %s found",i,symnam.data(),syna.data()));
426     SetSm(i,x);
427   }
428
429   // chambers
430
431   for (int i = 0; i < 540; i++) {
432     fi>>j>>x[0]>>x[1]>>x[2]>>x[3]>>x[4]>>x[5]>>volid>>syna;
433     if (j != i) AliError(Form("ch %d expected, %d found",i,j));
434     if (volid != GetVoi(i)) AliError(Form("ch %d volume id %d expected, %d found",i,GetVoi(i),volid));
435     std::string symnam = GetChName(i);
436     if (syna != symnam) AliError(Form("ch %d name %s expected, %s found",i,symnam.data(),syna.data()));
437     SetCh(i,x);
438   }
439
440   fi.close();
441
442 }
443
444 //_____________________________________________________________________________
445 void AliTRDalignment::ReadRoot(char *filename) 
446 {
447   //
448   // read the alignment data from root file
449   // here I expect a fixed order and number of elements
450   // it would be much better to identify the alignment objects 
451   // one by one and set the parameters of the corresponding sm or ch
452   //
453
454   TFile fi(filename,"READ");
455
456   if (fi.IsOpen()) {
457     TClonesArray *ar = (TClonesArray*) fi.Get("TRDAlignObjs");
458     ArToNumbers(ar);
459     fi.Close();
460   } 
461   else AliError(Form("cannot open input file %s",filename));
462
463   return;
464
465 }
466
467 //_____________________________________________________________________________
468 void AliTRDalignment::ReadDB(char *filename) 
469 {
470   //
471   // read the alignment data from database file
472   //
473
474   TFile fi(filename,"READ");
475
476   if (fi.IsOpen()) {
477     AliCDBEntry  *e  = (AliCDBEntry *) fi.Get("AliCDBEntry");
478     e->PrintMetaData();
479     fComment.SetString(e->GetMetaData()->GetComment());
480     TClonesArray *ar = (TClonesArray *) e->GetObject();
481     ArToNumbers(ar);
482     fi.Close();
483   } 
484   else AliError(Form("cannot open input file %s",filename));
485
486   return;
487
488 }
489
490 //_____________________________________________________________________________
491 void AliTRDalignment::ReadDB(char *db, char *path, int run
492                            , int version, int subversion)
493 {
494   //
495   // read the alignment data from database
496   //
497
498   AliCDBManager *cdb     = AliCDBManager::Instance();
499   AliCDBStorage *storLoc = cdb->GetStorage(db);
500   AliCDBEntry   *e       = storLoc->Get(path,run,version,subversion);
501   e->PrintMetaData();
502   fComment.SetString(e->GetMetaData()->GetComment());
503   TClonesArray  *ar      =  (TClonesArray *) e->GetObject();
504   ArToNumbers(ar);
505
506 }
507
508 //_____________________________________________________________________________
509 void AliTRDalignment::ReadGeo(char *misaligned) 
510 {
511   //
512   // determine misalignment by comparing original and misaligned matrix 
513   // of the last node on the misaligned_geometry file 
514   // an alternative longer way is in attic.C
515   //
516
517   TGeoHMatrix *ideSm[18];  // ideal
518   TGeoHMatrix *ideCh[540];
519   TGeoHMatrix *misSm[18];  // misaligned
520   TGeoHMatrix *misCh[540];
521
522   // read misaligned and original matrices
523
524   TGeoManager::Import(misaligned);
525   for (int i = 0; i < 18; i++) {
526     TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetSmName(i));
527     if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
528     TGeoPhysicalNode *node = pne->GetPhysicalNode();
529     if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i))); 
530     misSm[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
531     ideSm[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
532   }
533   for (int i = 0; i < 540; i++) {
534     TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetChName(i));
535     if (!pne) AliError(Form("no such physical node entry: %s",GetChName(i))); 
536     TGeoPhysicalNode *node = pne->GetPhysicalNode();
537     if (!node) AliError(Form("physical node entry %s has no physical node",GetChName(i))); 
538     misCh[i] = new TGeoHMatrix(*node->GetNode(node->GetLevel())->GetMatrix());
539     ideCh[i] = new TGeoHMatrix(*node->GetOriginalMatrix());
540   }
541
542   // calculate the local misalignment matrices as inverse misaligned times ideal
543
544   for (int i = 0; i < 18; i++) {
545     TGeoHMatrix mat(ideSm[i]->Inverse()); 
546     mat.Multiply(misSm[i]);
547     double *tra = mat.GetTranslation();
548     double *rot = mat.GetRotationMatrix();
549     double pars[6];
550     pars[0] = tra[0];
551     pars[1] = tra[1];
552     pars[2] = tra[2];
553     if (TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) AliError("Failed to extract roll-pitch-yall angles!");
554     double raddeg = TMath::RadToDeg();
555     pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
556     pars[4] = raddeg * TMath::ASin(rot[2]);
557     pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
558     SetSm(i,pars);
559   }
560
561   for (int i = 0; i < 540; i++) {
562     TGeoHMatrix mat(ideCh[i]->Inverse()); 
563     mat.Multiply(misCh[i]);
564     double *tra = mat.GetTranslation();
565     double *rot = mat.GetRotationMatrix();
566     double pars[6];
567     pars[0] = tra[0];
568     pars[1] = tra[1];
569     pars[2] = tra[2];
570     if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
571       AliError("Failed to extract roll-pitch-yall angles!");
572       return;
573     }
574     double raddeg = TMath::RadToDeg();
575     pars[3] = raddeg * TMath::ATan2(-rot[5],rot[8]);
576     pars[4] = raddeg * TMath::ASin(rot[2]);
577     pars[5] = raddeg * TMath::ATan2(-rot[1],rot[0]);
578     SetCh(i,pars);
579   }
580
581   // cleanup
582   for (int i = 0; i <  18; i++) delete ideSm[i];
583   for (int i = 0; i <  18; i++) delete misSm[i];
584   for (int i = 0; i < 540; i++) delete ideCh[i];
585   for (int i = 0; i < 540; i++) delete misCh[i];
586
587   return;
588
589 }
590
591 //_____________________________________________________________________________
592 void AliTRDalignment::ReadSurveyReport(char *filename) 
593 {
594   //
595   // Read survey report and store the numbers in fSurveyX, fSurveyY, fSurveyZ, 
596   // and fSurveyE.  Store the survey info in the fComment.
597   // Each supermodule has 8 survey points. The point names look like 
598   // TRD_sm08ah0 and have the following meaning. 
599   //
600   // sm00..17 mean supermodule 0 through 17, following the phi.
601   // Supermodule 00 is between phi=0 and phi=20 degrees.
602   //
603   // a or c denotes the anticlockwise and clockwise end of the supermodule
604   // in z. Clockwise end is where z is negative and where the muon arm sits.
605   //
606   // l or h denote low radius and high radius holes
607   //
608   // 0 or 1 denote the hole at smaller and at larger phi, respectively.
609   //
610
611   // read the survey file
612
613   fstream in(filename,fstream::in);
614   if (!in) {
615     AliError(Form("cannot open input file %s",filename));
616     return;
617   }
618
619   // loop through the lines of the file until the beginning of data
620
621   TString title,date,subdetector,url,version,observations,system,units;
622   while (1) {
623     char pee=in.peek();
624     if (pee==EOF) break; 
625     TString line;
626     line.ReadLine(in);
627     if (line.Contains("Title:"))        title.ReadLine(in);
628     if (line.Contains("Date:"))         date.ReadLine(in);
629     if (line.Contains("Subdetector:"))  subdetector.ReadLine(in);
630     if (line.Contains("URL:"))          url.ReadLine(in);
631     if (line.Contains("Version:"))      version.ReadLine(in);
632     if (line.Contains("Observations:")) observations.ReadLine(in);
633     if (line.Contains("System:"))       system.ReadLine(in);
634     if (line.Contains("Units:"))        units.ReadLine(in);
635     if (line.Contains("Data:"))         break;
636   }
637
638   // check what we found so far (watch out, they have \r at the end)
639
640   std::cout<<"title .........."<<title<<std::endl;
641   std::cout<<"date ..........."<<date<<std::endl;
642   std::cout<<"subdetector ...."<<subdetector<<std::endl;
643   std::cout<<"url ............"<<url<<std::endl;
644   std::cout<<"version ........"<<version<<std::endl;
645   std::cout<<"observations ..."<<observations<<std::endl;
646   std::cout<<"system ........."<<system<<std::endl;
647   std::cout<<"units .........."<<units<<std::endl;
648
649   if (!subdetector.Contains("TRD")) {
650     AliWarning(Form("Not a TRD survey file, subdetector = %s",subdetector.Data()));
651     return;
652   }
653   double tocm = 0; // we want to have it in cm
654   if (units.Contains("mm"))      tocm = 0.1;
655   else if (units.Contains("cm")) tocm = 1.0;
656   else if (units.Contains("m"))  tocm = 100.0;
657   else if (units.Contains("pc")) tocm = 3.24078e-15;
658   else {
659     AliError(Form("unexpected units: %s",units.Data()));
660     return;
661   }
662   if (!system.Contains("ALICEPH")) {
663     AliError(Form("wrong system: %s, should be ALICEPH",system.Data()));
664     return;
665   }
666
667   // scan the rest of the file which should contain list of surveyed points
668   // for every point, decode the point name and store the numbers in the right 
669   // place in the arrays fSurveyX etc.
670
671   while (1) {
672     TString pna; // point name
673     char type;
674     double x,y,z,precision;
675     in >> pna >> x >> y >> z >> type >> precision;  
676     if (in.fail()) break;
677     if (pna(0,6)!="TRD_sm") {
678       AliError(Form("unexpected point name: %s",pna.Data()));
679       break;
680     }
681     int i = atoi(pna(6,0).Data()); // supermodule number
682     int j = -1;
683     if (pna(8) == 'a') j=0; // anticlockwise, positive z
684     if (pna(8) == 'c') j=1; // clockwise, negative z
685     int k = -1;
686     if (pna(9) == 'l') k=0; // low radius
687     if (pna(9) == 'h') k=1; // high radius
688     int l = atoi(pna(10,0).Data()); // phi within supermodule
689     if (i>=0 && i<18 && j>=0 && j<2 && k>=0 && k<2 && l>=0 && l<2) {
690       fSurveyX[i][j][k][l] = tocm*x;
691       fSurveyY[i][j][k][l] = tocm*y;
692       fSurveyZ[i][j][k][l] = tocm*z;
693       fSurveyE[i][j][k][l] = precision/10; // "precision" is supposed to be in mm
694       std::cout << "decoded "<<pna<<" "
695                 <<fSurveyX[i][j][k][l]<<" "
696                 <<fSurveyY[i][j][k][l]<<" "
697                 <<fSurveyZ[i][j][k][l]<<" "
698                 <<fSurveyE[i][j][k][l]<<std::endl;      
699     } else AliError(Form("cannot decode point name: %s",pna.Data()));
700   }
701   in.close();
702   TString info = "Survey "+title+" "+date+" "+url+" "+version+" "+observations;
703   info.ReplaceAll("\r","");
704   fComment.SetString(info.Data());
705                          
706 }
707
708 //_____________________________________________________________________________
709 double AliTRDalignment::SurveyChi2(int i, double *a) {
710
711   //
712   // Compare the survey results to the ideal positions of the survey marks
713   // in the local frame of supermodule. When transforming, use the alignment 
714   // parameters a[6]. Return chi-squared.
715   //
716
717   if (!gGeoManager) LoadIdealGeometry();
718   printf("Survey of supermodule %d\n",i);
719   AliAlignObjAngles al(GetSmName(i),0,a[0],a[1],a[2],a[3],a[4],a[5],0);
720   TGeoPNEntry      *pne  = gGeoManager->GetAlignableEntry(GetSmName(i));
721   if (!pne) AliError(Form("no such physical node entry: %s",GetSmName(i)));
722   TGeoPhysicalNode *node = pne->GetPhysicalNode();
723   if (!node) AliError(Form("physical node entry %s has no physical node",GetSmName(i))); 
724
725   //  al.ApplyToGeometry();    
726   //  node = pne->GetPhysicalNode(); // changed in the meantime
727   //  TGeoHMatrix *ma = node->GetMatrix();
728
729   // a less destructive method (it does not modify geometry), gives the same result:
730
731   TGeoHMatrix *ma = new TGeoHMatrix();
732   al.GetLocalMatrix(*ma);
733   ma->MultiplyLeft(node->GetMatrix()); // global trafo, modified by a[]
734
735   double chi2=0;
736   printf("              sm   z   r  phi    x (lab phi)  y (lab z)   z (lab r)   all in cm\n");
737   for (int j=0; j<2; j++) for (int k=0; k<2; k++) for (int l=0; l<2; l++) {
738     if (fSurveyE[i][j][k][l] == 0.0) continue; // no data for this survey point
739     double master[3] = {fSurveyX[i][j][k][l],fSurveyY[i][j][k][l],fSurveyZ[i][j][k][l]};
740     double local[3];
741     ma->MasterToLocal(master,local);
742     double dx = local[0]-fSurveyX0[j][k][l];
743     double dy = local[1]-fSurveyY0[j][k][l];
744     double dz = local[2]-fSurveyZ0[j][k][l];
745     chi2 += (dx*dx+dy*dy+dz*dz)/fSurveyE[i][j][k][l]/fSurveyE[i][j][k][l];
746     printf("local survey %3d %3d %3d %3d %12.3f %12.3f %12.3f\n",i,j,k,l,local[0],local[1],local[2]);
747     printf("local ideal                  %12.3f %12.3f %12.3f\n",fSurveyX0[j][k][l],
748            fSurveyY0[j][k][l],fSurveyZ0[j][k][l]);
749     printf("difference                   %12.3f %12.3f %12.3f\n",dx,dy,dz);
750   }
751   printf("chi2 = %.2f\n",chi2);
752   return chi2;
753 }
754
755 //_____________________________________________________________________________
756 void Fcn(int &npar, double *g, double &f, double *par, int iflag) {
757
758   // 
759   // Standard function as needed by Minuit-like minimization procedures. 
760   // For the set of parameters par calculates and returns chi-squared.
761   //
762
763   // smuggle a C++ object into a C function
764   AliTRDalignment *alignment = (AliTRDalignment*) gMinuit->GetObjectFit(); 
765
766   f = alignment->SurveyChi2(par);
767   if (iflag==3);
768   if (npar); 
769   if (g); // no warnings about unused stuff...
770
771 }
772
773 //_____________________________________________________________________________
774 void AliTRDalignment::SurveyToAlignment(int i,char *flag) {
775
776   //
777   // Find the supermodule alignment parameters needed to make the survey 
778   // results coincide with the ideal positions of the survey marks.
779   // The string flag should look like "101000"; the six characters corresponds 
780   // to the six alignment parameters and 0/1 mean that the parameter should 
781   // be fixed/released in the fit. 
782
783   if (strlen(flag)!=6) {
784     AliError(Form("unexpected flag: %s",flag));
785     return;
786   }
787
788   printf("Finding alignment matrix for supermodule %d\n",i);
789   fIbuffer[0] = i; // store the sm number in the buffer so minuit can see it
790
791   TFitter fitter(100);
792   gMinuit->SetObjectFit(this);
793   fitter.SetFCN(Fcn);
794   fitter.SetParameter(0,"dx",0,0.5,0,0);
795   fitter.SetParameter(1,"dy",0,0.5,0,0);
796   fitter.SetParameter(2,"dz",0,0.5,0,0);
797   fitter.SetParameter(3,"rx",0,0.1,0,0);
798   fitter.SetParameter(4,"ry",0,0.1,0,0);
799   fitter.SetParameter(5,"rz",0,0.1,0,0);
800
801   for (int j=0; j<6; j++) if (flag[j]=='0') fitter.FixParameter(j);
802
803   double arglist[100];
804   arglist[0] = 2;
805   fitter.ExecuteCommand("SET PRINT", arglist, 1);
806   fitter.ExecuteCommand("SET ERR", arglist, 1);
807   arglist[0]=50;
808   //fitter.ExecuteCommand("SIMPLEX", arglist, 1);
809   fitter.ExecuteCommand("MINIMIZE", arglist, 1);
810   fitter.ExecuteCommand("CALL 3", arglist,0);
811   double a[6];
812   for (int j=0; j<6; j++) a[j] = fitter.GetParameter(j);
813   SetSm(i,a);
814   for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParameter(j));   
815   printf("\n");
816   for (int j=0; j<6; j++) printf("%10.3f ",fitter.GetParError(j));
817   printf("\n");
818
819 }
820
821 //_____________________________________________________________________________
822 void AliTRDalignment::ReadAny(char *filename) 
823 {
824   //
825   // read the alignment data from any kind of file
826   //
827
828   TString fist(filename);
829   if (fist.EndsWith(".txt")) ReadAscii(filename);
830   if (fist.EndsWith(".dat")) ReadAscii(filename);
831   if (fist.EndsWith(".root")) {
832     if (fist.Contains("Run")) ReadDB(filename);
833     else ReadRoot(filename);
834   }
835
836 }
837
838 //_____________________________________________________________________________
839 void AliTRDalignment::WriteAscii(char *filename) const
840 {
841   //
842   // store the alignment data on ascii file
843   //
844
845   FILE *fp = fopen(filename, "w");
846   if (!fp) {
847     AliError(Form("cannot open output file %s",filename));
848     return;
849   }
850
851   PrintSm(fp);
852   PrintCh(fp);
853   
854   fclose(fp);
855
856 }
857
858 //_____________________________________________________________________________
859 void AliTRDalignment::WriteRoot(char *filename) 
860 {
861   //
862   // store the alignment data on root file
863   //
864
865   TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
866   NumbersToAr(ar);
867   TFile fo(filename,"RECREATE");
868   if (fo.IsOpen()) {
869     fo.cd();
870     fo.WriteObject(ar,"TRDAlignObjs","kSingleKey");
871     fo.Close();
872   } 
873   else AliError(Form("cannot open output file %s",filename));
874
875   delete ar;
876
877 }
878
879 //_____________________________________________________________________________
880 void AliTRDalignment::WriteDB(char *filename, int run0, int run1) 
881 {
882   //
883   // dumping on a DB-like file
884   //
885
886   TClonesArray   *ar = new TClonesArray("AliAlignObjAngles",10000);
887   NumbersToAr(ar);
888   char *path = "di1/di2/di3";
889   AliCDBId id(path,run0,run1);
890   AliCDBMetaData *md = new AliCDBMetaData();
891   md->SetResponsible("Dariusz Miskowiec");
892   md->SetComment(fComment.GetString().Data());
893   AliCDBEntry    *e  = new AliCDBEntry(ar, id, md);
894   TFile fi(filename,"RECREATE");
895   if (fi.IsOpen()) {
896     e->Write();
897     fi.Close();
898   } 
899   else AliError(Form("cannot open input file %s",filename));
900
901   delete e;
902   delete md;
903   delete ar;
904
905   return;
906
907 }
908
909 //_____________________________________________________________________________
910 void AliTRDalignment::WriteDB(char *db, char *path, int run0, int run1) 
911 {
912   //
913   // store the alignment data in database
914   //
915
916   TClonesArray   *ar      = new TClonesArray("AliAlignObjAngles",10000);
917   NumbersToAr(ar);
918   AliCDBManager  *cdb     = AliCDBManager::Instance();
919   AliCDBStorage  *storLoc = cdb->GetStorage(db);
920   AliCDBMetaData *md      = new AliCDBMetaData();
921   md->SetResponsible("Dariusz Miskowiec");
922   md->SetComment(fComment.GetString().Data());
923   AliCDBId id(path,run0,run1);
924   storLoc->Put(ar,id,md);
925   md->Delete();
926   delete ar;
927
928 }
929
930 //_____________________________________________________________________________
931 void AliTRDalignment::WriteGeo(char *filename) 
932 {
933   //
934   // apply misalignment to (currently loaded ideal) geometry and store the 
935   // resulting geometry on a root file
936   //
937
938   TClonesArray *ar = new TClonesArray("AliAlignObjAngles",10000);
939   NumbersToAr(ar);
940   for (int i = 0; i < ar->GetEntriesFast(); i++) {
941     AliAlignObj *alobj = (AliAlignObj *) ar->UncheckedAt(i);
942     alobj->ApplyToGeometry();
943   }
944   delete ar;
945   gGeoManager->Export(filename);
946
947 }
948
949 //_____________________________________________________________________________
950 double AliTRDalignment::GetSmRMS(int xyz) const 
951 {
952   //
953   // rms fSm[][xyz]
954   //
955
956   double s1 = 0.0;
957   double s2 = 0.0;
958   for (int i = 0; i < 18; i++) {
959     s1 += fSm[i][xyz];
960     s2 += fSm[i][xyz]*fSm[i][xyz];
961   }
962   double rms2 = s2/18.0 - s1*s1/18.0/18.0;
963
964   return rms2>0 ? sqrt(rms2) : 0.0;
965
966 }
967
968 //_____________________________________________________________________________
969 double AliTRDalignment::GetChRMS(int xyz) const
970 {
971   //
972   // rms fCh[][xyz]
973   //
974
975   double s1 =0.0;
976   double s2 =0.0;
977   for (int i = 0; i < 540; i++) {
978     s1 += fCh[i][xyz];
979     s2 += fCh[i][xyz]*fCh[i][xyz];
980   }
981   double rms2 = s2/540.0 - s1*s1/540.0/540.0;
982
983   return rms2>0 ? sqrt(rms2) : 0.0;
984
985 }
986
987 //_____________________________________________________________________________
988 void AliTRDalignment::PrintSmRMS() const
989 {
990   //
991   // dump rms of fSm
992   //
993
994   printf("       %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f  supermodule rms\n"
995         ,GetSmRMS(0),GetSmRMS(1),GetSmRMS(2),GetSmRMS(3),GetSmRMS(4),GetSmRMS(5));
996
997 }
998
999 //_____________________________________________________________________________
1000 void AliTRDalignment::PrintChRMS() const
1001 {
1002   //
1003   // dump rms of fCh
1004   //
1005
1006   printf("       %11.4f %11.4f  %11.4f      %11.5f  %11.5f  %11.5f  chamber rms\n"
1007         ,GetChRMS(0),GetChRMS(1),GetChRMS(2),GetChRMS(3),GetChRMS(4),GetChRMS(5));
1008
1009 }
1010
1011 //_____________________________________________________________________________
1012 void AliTRDalignment::ArToNumbers(TClonesArray *ar) 
1013 {
1014   //
1015   // read numbers from the array of AliAlignObj objects and fill fSm and fCh
1016   //
1017
1018   LoadIdealGeometry();
1019   SetZero();
1020   double pa[6];
1021   for (int i = 0; i <  18; i++) {
1022     AliAlignObj *aao = (AliAlignObj *) ar->At(i);
1023     aao->GetLocalPars(pa,pa+3);
1024     SetSm(i,pa);
1025   }
1026   for (int i = 0; i < 540; i++) {
1027     AliAlignObj *aao = (AliAlignObj *) ar->At(18+i);
1028     aao->GetLocalPars(pa,pa+3);
1029     SetCh(i,pa);
1030   }
1031
1032 }
1033
1034 //_____________________________________________________________________________
1035 void AliTRDalignment::NumbersToAr(TClonesArray *ar) 
1036 {
1037   //
1038   // build array of AliAlignObj objects based on fSm and fCh data
1039   //
1040
1041   LoadIdealGeometry();
1042   TClonesArray &alobj = *ar;
1043   int nobj = 0;
1044   for (int i = 0; i <  18; i++) {      
1045       new(alobj[nobj]) AliAlignObjAngles(GetSmName(i)
1046                                         ,0 
1047                                         ,fSm[i][0],fSm[i][1],fSm[i][2]
1048                                         ,fSm[i][3],fSm[i][4],fSm[i][5]
1049                                         ,0);
1050     nobj++;
1051   }
1052
1053   for (int i = 0; i < 540; i++) {
1054     new(alobj[nobj]) AliAlignObjAngles(GetChName(i)
1055                                       ,GetVoi(i)
1056                                       ,fCh[i][0],fCh[i][1],fCh[i][2]
1057                                       ,fCh[i][3],fCh[i][4],fCh[i][5]
1058                                       ,0);
1059     nobj++;
1060   }
1061
1062 }
1063
1064 //_____________________________________________________________________________
1065 void AliTRDalignment::LoadIdealGeometry(char *filename) 
1066 {
1067   //
1068   // load ideal geometry from filename
1069   // it is needed for operations on AliAlignObj objects
1070   // this needs to be straightened out
1071   // particularly, sequences LoadIdealGeometry("file1"); LoadIdealGeometry("file2"); 
1072   // do not work as one would naturally expect
1073   //
1074
1075   static int attempt = 0; // which reload attempt is it? just to avoid endless loops
1076
1077   if (!gGeoManager) TGeoManager::Import(filename);
1078   if (!gGeoManager) AliError(Form("cannot open geometry file %s",filename));
1079   if (gGeoManager->GetListOfPhysicalNodes()->GetEntries()) {
1080     if (attempt) AliError(Form("geometry on file %s is not ideal",filename));
1081     AliWarning("current geometry is not ideal - it contains physical nodes");
1082     AliWarning(Form("reloading geometry from %s - attempt nr %d",filename,attempt));
1083     gGeoManager = 0;
1084     attempt++;
1085     LoadIdealGeometry(filename);
1086   }
1087
1088   attempt = 0;
1089
1090 }