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