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