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