da documentation page is defined in the Link
[u/mrichter/AliRoot.git] / PMD / PMDSurveyToAlignment_v1.C
CommitLineData
c991220d 1void PMDSurveyToAlignment_v1(){
2
3 // Macro to convert survey data into alignment data.
4 // The position of four fiducial marks, sticked on the
5 // preshower plane of PMD is converted into the
6 // global position.
7
8 if(!gGeoManager) AliGeomManager::LoadGeometry("geometry.root");
9
10 TClonesArray *array = new TClonesArray("AliAlignObjMatrix",10);
11 TClonesArray &mobj = *array;
12
13 Double_t l_vect[3]={0.,0.,0.}; // a local vector (the origin)
14
15 Double_t g_vect_1[3]; // vector corresp. to it in global RS for sector-1
16 Double_t g_vect_2[3]; // vector corresp. to it in global RS for sector-2
17 Double_t g_vect_3[3]; // vector corresp. to it in global RS for sector-3
18 Double_t g_vect_4[3]; // vector corresp. to it in global RS for sector-4
19
20 //**************** get global matrix *******************
21
22 TGeoHMatrix *g3_1 = AliGeomManager::GetMatrix("PMD/Sector1");
23 TGeoNode* n3_1 = gGeoManager->GetCurrentNode();
24 TGeoHMatrix* l3_1 = n3_1->GetMatrix();// to get local matrix
25 g3_1->LocalToMaster(l_vect,g_vect_1); // point coordinates in the global RS
26
27 cout<<endl<<"Origin of sector-1 in the global RS: "<<
28 g_vect_1[0]<<" "<<g_vect_1[1]<<" "<<g_vect_1[2]<<" "<<endl;
29
30
31 TGeoHMatrix *g3_2 = AliGeomManager::GetMatrix("PMD/Sector2");
32 TGeoNode* n3_2 = gGeoManager->GetCurrentNode();
33 TGeoHMatrix* l3_2 = n3_2->GetMatrix();
34 g3_2->LocalToMaster(l_vect,g_vect_2);
35
36 cout<<endl<<"Origin of sector-2 in the global RS: "<<
37 g_vect_2[0]<<" "<<g_vect_2[1]<<" "<<g_vect_2[2]<<" "<<endl;
38
39
40 TGeoHMatrix *g3_3 = AliGeomManager::GetMatrix("PMD/Sector3");
41 TGeoNode* n3_3 = gGeoManager->GetCurrentNode();
42 TGeoHMatrix* l3_3 = n3_3->GetMatrix();
43 g3_3->LocalToMaster(l_vect,g_vect_3);
44
45 cout<<endl<<"Origin of the sector-3 in the global RS: "<<
46 g_vect_3[0]<<" "<<g_vect_3[1]<<" "<<g_vect_3[2]<<" "<<endl;
47
48
49 TGeoHMatrix *g3_4 = AliGeomManager::GetMatrix("PMD/Sector4");
50 TGeoNode* n3_4 = gGeoManager->GetCurrentNode();
51 TGeoHMatrix* l3_4 = n3_4->GetMatrix();
52 g3_4->LocalToMaster(l_vect,g_vect_4);
53
54 cout<<endl<<"Origin of the sector-4 in the global RS: "<<
55 g_vect_4[0]<<" "<<g_vect_4[1]<<" "<<g_vect_4[2]<<" "<<endl;
56
57
58 // From the coordinates of fiducial marks derive back the
59 // new global position of the surveyed volume.
60 // What follows is the actual survey-to-alignment procedure which assumes,
61 // 4 fiducial marks at the corners of a rectangle lying on a plane parallel
62 // to a surface of the detector at a certain offset and with x and y sides
63 // parallel to the detector's x and y axes.
64 // It needs as input of four points and the offset from the origin (zdepth).
65 // The algorithm can be easily modified for different placement and/or
66 // cardinality of the fiducial marks.
67
68 // Now specify the path of the module to be misaligned
69 // as followed for the PMD geometry in Geant.
70
71 /*
72 Y|
73 |
74 _____________ D(3)|_____________C(2)
75 | | | | |
76 | 1 | 3 | | |
77 | |________| | |
78 |____|___| | | |
79 | | 2 | | |
80 | 4 | | |_____________|_______
81 |________|____| A(0) B(1) X
82
83 //
84 // Misalignment Matrix is set for sector-1, sector-2, sector-3 and
85 // sector-4 even though sectors 1 and 4 and sectors 2 and 3
86 // will be mounted on the same steel plate.
87 */
88
89// All the units are in milimeter in survey data file.
90
91// Retrieval of real survey data from ALICE Survey Data Depot :
92
93 AliSurveyObj *so = new AliSurveyObj();
94 //so->FillFromLocalFile("Survey_Points_PMD01.txt");//this file will be given
95 //by the surveyer after the installation of PMD.
96
97 so->FillFromLocalFile("PMDGenSurveyPoints_v1.txt");//This file can be generated by the user.
98
99 Int_t size = so->GetEntries();
100 Printf("Title: \"%s\"", so->GetReportTitle().Data());
101 Printf("Date: \"%s\"", so->GetReportDate().Data());
102 Printf("Detector: \"%s\"", so->GetDetector().Data());
103 Printf("URL: \"%s\"", so->GetURL().Data());
104 Printf("Number: \"%d\"", so->GetReportNumber());
105 Printf("Version: \"%d\"", so->GetReportVersion());
106 Printf("Observations: \"%s\"", so->GetObservations().Data());
107 Printf("Coordinate System: \"%s\"", so->GetCoordSys().Data());
108 Printf("Measurement Units: \"%s\"", so->GetUnits().Data());
109 Printf("Nr Columns: \"%d\" \n", so->GetNrColumns());
110
111 TObjArray *colNames = so->GetColumnNames();
112 TObjArray *points = so->GetData();
113
114 Printf("Relevant points to be used for alignment procedure (in mm):");
115 printf("%d \n",points->GetEntries());
116
117 Float_t surveyFidX[28];
118 Float_t surveyFidY[28];
119 Float_t surveyFidZ[28];
120
121 for (Int_t i = 0; i < points->GetEntries(); ++i)
122 {
123
124 surveyFidX[i] = (((AliSurveyPoint *) points->At(i))->GetX());
125 surveyFidY[i] = (((AliSurveyPoint *) points->At(i))->GetY());
126 surveyFidZ[i] = (((AliSurveyPoint *) points->At(i))->GetZ());
127 }
128
129 for (Int_t i=0; i < 28; i++)
130 {
131 // printf("%f %f %f\n",surveyFidX[i],surveyFidY[i],surveyFidZ[i]);
132
133 }
134
135 const Int_t kNDIM1 = 3;
136 const Int_t kNDIM2 = 4;
137 const Double_t kIdealOrig[kNDIM1] = {0., 0., 3645.}; // Geant values
138
139 Int_t i;
140
141 Double_t s1;
142 Double_t ab1[kNDIM1], bc1[kNDIM1], nn1[kNDIM1], plane1[kNDIM1];
143 Double_t nga1[kNDIM2][kNDIM1];
144
145 Double_t s2;
146 Double_t ab2[kNDIM1], bc2[kNDIM1], nn2[kNDIM1], plane2[kNDIM1];
147 Double_t nga2[kNDIM2][kNDIM1];
148
149 Double_t s3;
150 Double_t ab3[kNDIM1], bc3[kNDIM1], nn3[kNDIM1], plane3[kNDIM1];
151 Double_t nga3[kNDIM2][kNDIM1];
152
153 Double_t s4;
154 Double_t ab4[kNDIM1], bc4[kNDIM1], nn4[kNDIM1], plane4[kNDIM1];
155 Double_t nga4[kNDIM2][kNDIM1];
156
157 // These are the sequence of fiducial marks for sector-1,2,3 & 4 - to be
158 // given by the user.
159
160 Int_t surseq1[kNDIM2] = {19, 20, 18, 17};
161 Int_t surseq2[kNDIM2] = {15, 16, 14, 13};
162 Int_t surseq3[kNDIM2] = {11, 12, 06, 05};
163 Int_t surseq4[kNDIM2] = {27, 28, 22, 21};
164
165 for (i = 0; i < 28; i++)
166 {
167 for (Int_t j = 0; j < 4; j++)
168 {
169 if (surseq1[j] == i+1)
170 {
171 nga1[j][0] = surveyFidX[i];
172 nga1[j][1] = surveyFidY[i];
173 nga1[j][2] = surveyFidZ[i];
174 }
175 else if (surseq2[j] == i+1)
176 {
177 nga2[j][0] = surveyFidX[i];
178 nga2[j][1] = surveyFidY[i];
179 nga2[j][2] = surveyFidZ[i];
180 }
181 else if (surseq3[j] == i+1)
182 {
183 nga3[j][0] = surveyFidX[i];
184 nga3[j][1] = surveyFidY[i];
185 nga3[j][2] = surveyFidZ[i];
186 }
187 else if (surseq4[j] == i+1)
188 {
189 nga4[j][0] = surveyFidX[i];
190 nga4[j][1] = surveyFidY[i];
191 nga4[j][2] = surveyFidZ[i];
192 }
193 }
194
195 }
196
197 // First vector on the plane of the fiducial marks.
198
199 for(i = 0; i < kNDIM1; i++)
200 {
201 ab1[i] = nga1[1][i] - nga1[0][i];
202 ab2[i] = nga2[1][i] - nga2[0][i];
203 ab3[i] = nga3[1][i] - nga3[0][i];
204 ab4[i] = nga4[1][i] - nga4[0][i];
205 }
206
207 // Second vector on the plane of the fiducial marks.
208
209 for(i = 0; i < kNDIM1; i++)
210 {
211 bc1[i] = nga1[2][i] - nga1[1][i];
212 bc2[i] = nga2[2][i] - nga2[1][i];
213 bc3[i] = nga3[2][i] - nga3[1][i];
214 bc4[i] = nga4[2][i] - nga4[1][i];
215 }
216
217 // Vector normal to the plane of the fiducial marks obtained
218 // as cross product of the two vectors on the sector-1,2,3 and 4.
219
220 // Vector normal to the sector-1.
221
222 nn1[0] = ab1[1] * bc1[2] - ab1[2] * bc1[1];
223 nn1[1] = ab1[2] * bc1[0] - ab1[0] * bc1[2];
224 nn1[2] = ab1[0] * bc1[1] - ab1[1] * bc1[0];
225
226 // Vector normal to the sector-2.
227
228 nn2[0] = ab2[1] * bc2[2] - ab2[2] * bc2[1];
229 nn2[1] = ab2[2] * bc2[0] - ab2[0] * bc2[2];
230 nn2[2] = ab2[0] * bc2[1] - ab2[1] * bc2[0];
231
232 // Vector normal to the sector-3.
233
234 nn3[0] = ab3[1] * bc3[2] - ab3[2] * bc3[1];
235 nn3[1] = ab3[2] * bc3[0] - ab3[0] * bc3[2];
236 nn3[2] = ab3[0] * bc3[1] - ab3[1] * bc3[0];
237
238 // Vector normal to the sector-4.
239
240 nn4[0] = ab4[1] * bc4[2] - ab4[2] * bc4[1];
241 nn4[1] = ab4[2] * bc4[0] - ab4[0] * bc4[2];
242 nn4[2] = ab4[0] * bc4[1] - ab4[1] * bc4[0];
243
244 Double_t sizen1 = TMath::Sqrt( nn1[0]*nn1[0] + nn1[1]*nn1[1] +
245 nn1[2]*nn1[2] );
246 Double_t sizen2 = TMath::Sqrt( nn2[0]*nn2[0] + nn2[1]*nn2[1] +
247 nn2[2]*nn2[2] );
248 Double_t sizen3 = TMath::Sqrt( nn3[0]*nn3[0] + nn3[1]*nn3[1] +
249 nn3[2]*nn3[2] );
250 Double_t sizen4 = TMath::Sqrt( nn4[0]*nn4[0] + nn4[1]*nn4[1] +
251 nn4[2]*nn4[2] );
252
253 if(sizen1 > 1.e-8)
254 {
255 s1 = Double_t(1.)/sizen1 ; //normalization factor
256 }
257 else
258 {
259 return 0;
260 }
261
262 if(sizen2 > 1.e-8)
263 {
264 s2 = Double_t(1.)/sizen2 ; //normalization factor
265 }
266 else
267 {
268 return 0;
269 }
270
271 if(sizen3 > 1.e-8)
272 {
273 s3 = Double_t(1.)/sizen3 ; //normalization factor
274 }
275 else
276 {
277 return 0;
278 }
279
280 if(sizen4 > 1.e-8)
281 {
282 s4 = Double_t(1.)/sizen4 ; //normalization factor
283 }
284 else
285 {
286 return 0;
287 }
288
289 // Plane expressed in the hessian normal form, see:
290 // http://mathworld.wolfram.com/HessianNormalForm.html
291
292 for(i = 0; i < kNDIM1; i++)
293 {
294 plane1[i] = nn1[i] * s1;
295 plane2[i] = nn2[i] * s2;
296 plane3[i] = nn3[i] * s3;
297 plane4[i] = nn4[i] * s4;
298 }
299
300 cout<<"Unit vector normal to sector-1 : "<<plane1[0]
301 <<" "<<plane1[1]<<" "<<plane1[2]<<endl;
302 cout<<"Unit vector normal to sector-2 : "<<plane2[0]
303 <<" "<<plane2[1]<<" "<<plane2[2]<<endl;
304 cout<<"Unit vector normal to sector-3 : "<<plane3[0]
305 <<" "<<plane3[1]<<" "<<plane3[2]<<endl;
306 cout<<"Unit vector normal to sector-4 : "<<plane4[0]
307 <<" "<<plane4[1]<<" "<<plane4[2]<<endl;
308
309 // The center of the sector with fiducial marks as corners,
310 // is the middle point of one diagonal.
311
312 Double_t orig1[kNDIM1], md1[kNDIM1];
313 Double_t orig2[kNDIM1], md2[kNDIM1];
314 Double_t orig3[kNDIM1], md3[kNDIM1];
315 Double_t orig4[kNDIM1], md4[kNDIM1];
316
317 for(i = 0; i < kNDIM1; i++)
318 {
319 md1[i] = (nga1[0][i] + nga1[1][i]) * 0.5;
320 md2[i] = (nga2[2][i] + nga2[3][i]) * 0.5;
321 md3[i] = (nga3[0][i] + nga3[2][i]) * 0.5;
322 md4[i] = (nga4[0][i] + nga4[2][i]) * 0.5;
323 }
324
325 const Double_t kZdepth = 23.;//(Zdepth+zoffset) is the distance
326 //from surface to centre of the PMD.
327
328 for(i = 0; i < kNDIM1; i++)
329 {
330 orig1[i] = md1[i] - plane1[i]*kZdepth;
331 orig2[i] = md2[i] - plane2[i]*kZdepth;
332 orig3[i] = md3[i] - plane3[i]*kZdepth;
333 orig4[i] = md4[i] - plane4[i]*kZdepth;
334 }
335
336 Double_t PMDorig1[kNDIM1];
337 Double_t PMDorig2[kNDIM1];
338 for(i = 0; i < kNDIM1; i++)
339 {
340 PMDorig1[i] = (orig1[i] + orig2[i]) * 0.5;
341 PMDorig2[i] = (orig3[i] + orig4[i]) * 0.5;
342 }
343
344 // Get x,y local directions needed to write the global rotation matrix
345 // for the surveyed volume by normalising vectors ab1 and bc1
346
347 Double_t sx1 = TMath::Sqrt(ab1[0]*ab1[0] + ab1[1]*ab1[1] +
348 ab1[2]*ab1[2]);
349 if(sx1 > 1.e-8)
350 {
351 for(i = 0; i < 3; i++)
352 {
353 ab1[i] /= sx1;
354 }
355 }
356
357 Double_t sy1 = TMath::Sqrt(bc1[0]*bc1[0] + bc1[1]*bc1[1] +
358 bc1[2]*bc1[2]);
359 if(sy1 > 1.e-8)
360 {
361 for(i = 0; i < kNDIM1; i++)
362 {
363 bc1[i] /= sy1;
364 }
365 }
366
367 Double_t sx2 = TMath::Sqrt(ab2[0]*ab2[0] + ab2[1]*ab2[1] +
368 ab2[2]*ab2[2]);
369 if(sx2 > 1.e-8)
370 {
371 for(i = 0; i < kNDIM1; i++)
372 {
373 ab2[i] /= sx2;
374 }
375 }
376
377 Double_t sy2 = TMath::Sqrt(bc2[0]*bc2[0] + bc2[1]*bc2[1] +
378 bc2[2]*bc2[2]);
379 if(sy2 > 1.e-8)
380 {
381 for(i = 0; i < kNDIM1; i++)
382 {
383 bc2[i] /= sy2;
384 }
385 }
386
387 Double_t sx3 = TMath::Sqrt(ab3[0]*ab3[0] + ab3[1]*ab3[1] +
388 ab3[2]*ab3[2]);
389 if(sx3 > 1.e-8)
390 {
391 for(i = 0; i < kNDIM1; i++)
392 {
393 ab3[i] /= sx3;
394 }
395 }
396
397 Double_t sy3 = TMath::Sqrt(bc3[0]*bc3[0] + bc3[1]*bc3[1] +
398 bc3[2]*bc3[2]);
399 if(sy3 > 1.e-8)
400 {
401 for(i = 0; i < kNDIM1; i++)
402 {
403 bc3[i] /= sy3;
404 }
405 }
406
407 Double_t sx4 = TMath::Sqrt(ab4[0]*ab4[0] + ab4[1]*ab4[1] +
408 ab4[2]*ab4[2]);
409 if(sx4 > 1.e-8)
410 {
411 for(i = 0; i < kNDIM1; i++)
412 {
413 ab4[i] /= sx4;
414 }
415 }
416
417 Double_t sy4 = TMath::Sqrt(bc4[0]*bc4[0] + bc4[1]*bc4[1] +
418 bc4[2]*bc4[2]);
419 if(sy4 > 1.e-8)
420 {
421 for(i = 0; i < kNDIM1; i++)
422 {
423 bc4[i] /= sy4;
424 }
425 }
426
427 // The global matrix for the surveyed volume - ng_1
428
429 Double_t rot1[9] = { ab1[0], bc1[0], plane1[0],
430 ab1[1], bc1[1], plane1[1],
431 ab1[2], bc1[2], plane1[2]};
432
433 TGeoHMatrix ng_1;
434 for(Int_t i=0; i<3; i++) orig1[i]*=0.1;// To convert the surveyed origin
435 //in cm, as in geometry all the distances are in cm.
436
437 orig1[0]=orig1[0]- 0.321402;//offset in X.
438 orig1[1]=orig1[1]- 0.074998;//offset in Y.
439 ng_1.SetTranslation(orig1);
440 ng_1.SetRotation(rot1);
441
442 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-1 **\n";
443 //ng_1.Print();
444
445 Double_t rot2[9] = {ab2[0], bc2[0], plane2[0],
446 ab2[1], bc2[1], plane2[1],
447 ab2[2], bc2[2], plane2[2]};
448
449 TGeoHMatrix ng_2;
450 for(Int_t i=0; i<3; i++) orig2[i]*=0.1;
451 orig2[0]=orig2[0]+ 0.321402;//offset in X.
452 orig2[1]=orig2[1]+ 0.074998;//offset in Y.
453 ng_2.SetTranslation(orig2);
454 ng_2.SetRotation(rot2);
455
456 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-2 **\n";
457 // ng_2.Print();
458
459 Double_t rot3[9] = {ab3[0], bc3[0], plane3[0],
460 ab3[1], bc3[1], plane3[1],
461 ab3[2], bc3[2], plane3[2]};
462
463 TGeoHMatrix ng_3;
464 for(Int_t i=0; i<3; i++) orig3[i]*=0.1;
465 orig3[0]=orig3[0]+ 0.291602;//offset in X
466 orig3[1]=orig3[1]+ 0.100001;//offset in Y
467 ng_3.SetTranslation(orig3);
468 ng_3.SetRotation(rot3);
469
470 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-3 **\n";
471 //ng_3.Print();
472
473
474 Double_t rot4[9] = {ab4[0], bc4[0], plane4[0],
475 ab4[1], bc4[1], plane4[1],
476 ab4[2], bc4[2], plane4[2]};
477
478 TGeoHMatrix ng_4;
479 for(Int_t i=0; i<3; i++) orig4[i]*=0.1;
480 orig4[0]=orig4[0]- 0.291602;//offset in X
481 orig4[1]=orig4[1]- 0.100001;//offset in Y
482 ng_4.SetTranslation(orig4);
483 ng_4.SetRotation(rot4);
484
485 //cout<<"\n** Global Matrix from surveyed fiducial marks for sector-4 **\n";
486 //ng_4.Print();
487
488
489 // To produce the alignment object for the given volume you would
490 // then do something like this:
491 // Calculate the global delta transformation as ng * g3^-1
492
493 TGeoPhysicalNode* pn3_1 = gGeoManager->MakePhysicalNode("ALIC_1/EPM1_1");
494 TGeoHMatrix gdelta_1 = g3_1->Inverse(); //now equal to the inverse of g3
495 gdelta_1.MultiplyLeft(&ng_1);
496 Int_t index = 0;
497 Printf("The global delta transformation for sector-1");
498 gdelta_1.Print();
499
500 TGeoPhysicalNode* pn3_2 = gGeoManager->MakePhysicalNode("ALIC_1/EPM2_1");
501 TGeoHMatrix gdelta_2 = g3_2->Inverse();
502 gdelta_2.MultiplyLeft(&ng_2);
503 Int_t index = 0;
504 Printf("The global delta transformation for sector-2");
505 gdelta_2.Print();
506
507 TGeoPhysicalNode* pn3_3 = gGeoManager->MakePhysicalNode("ALIC_1/EPM3_1");
508 TGeoHMatrix gdelta_3 = g3_3->Inverse();
509 gdelta_3.MultiplyLeft(&ng_3);
510 Int_t index = 0;
511 Printf("The global delta transformation for sector-3");
512 gdelta_3.Print();
513
514 TGeoPhysicalNode* pn3_4 = gGeoManager->MakePhysicalNode("ALIC_1/EPM4_1");
515 TGeoHMatrix gdelta_4 = g3_4->Inverse();
516 gdelta_4.MultiplyLeft(&ng_4);
517 Int_t index = 0;
518 Printf("The global delta transformation for sector-4");
519 gdelta_4.Print();
520
521 new(mobj[0]) AliAlignObjMatrix("PMD/Sector1",index,gdelta_1,kTRUE);
522 new(mobj[1]) AliAlignObjMatrix("PMD/Sector2",index,gdelta_2,kTRUE);
523 new(mobj[2]) AliAlignObjMatrix("PMD/Sector3",index,gdelta_3,kTRUE);
524 new(mobj[3]) AliAlignObjMatrix("PMD/Sector4",index,gdelta_4,kTRUE);
525
526
527 TString strStorage(gSystem->Getenv("TOCDB"));
528 if(strStorage != TString("kTRUE")){
529 // save on file
530 TFile f("PMDSurvey.root","RECREATE");
531 if(!f) cerr<<"cannot open file for output\n";
532 f.cd();
533 f.WriteObject(array,"PMDSurveyObjs ","kSingleKey");
534 f.Close();
535 }else{
536 // save in CDB storage
537 AliCDBManager* cdb = AliCDBManager::Instance();
162637e4 538 AliCDBStorage* storage = cdb->GetStorage("local://$ALICE_ROOT/OCDB");
c991220d 539 AliCDBMetaData* mda = new AliCDBMetaData();
540 mda->SetResponsible(" ");
541 mda->SetComment("Alignment objects for PMD survey");
542 AliCDBId id("PMD/Align/Data",0,AliCDBRunRange::Infinity());
543 storage->Put(array,id,mda);
544 }
545
546 array->Delete();
547}