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