1 void PMDSurveyToAlignment_v1(){
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
8 if(!gGeoManager) AliGeomManager::LoadGeometry("geometry.root");
10 TClonesArray *array = new TClonesArray("AliAlignObjMatrix",10);
11 TClonesArray &mobj = *array;
13 Double_t l_vect[3]={0.,0.,0.}; // a local vector (the origin)
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
20 //**************** get global matrix *******************
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
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;
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);
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;
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);
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;
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);
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;
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.
68 // Now specify the path of the module to be misaligned
69 // as followed for the PMD geometry in Geant.
74 _____________ D(3)|_____________C(2)
80 | 4 | | |_____________|_______
81 |________|____| A(0) B(1) X
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.
89 // All the units are in milimeter in survey data file.
91 // Retrieval of real survey data from ALICE Survey Data Depot :
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.
97 so->FillFromLocalFile("PMDGenSurveyPoints_v1.txt");//This file can be generated by the user.
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());
111 TObjArray *colNames = so->GetColumnNames();
112 TObjArray *points = so->GetData();
114 Printf("Relevant points to be used for alignment procedure (in mm):");
115 printf("%d \n",points->GetEntries());
117 Float_t surveyFidX[28];
118 Float_t surveyFidY[28];
119 Float_t surveyFidZ[28];
121 for (Int_t i = 0; i < points->GetEntries(); ++i)
124 surveyFidX[i] = (((AliSurveyPoint *) points->At(i))->GetX());
125 surveyFidY[i] = (((AliSurveyPoint *) points->At(i))->GetY());
126 surveyFidZ[i] = (((AliSurveyPoint *) points->At(i))->GetZ());
129 for (Int_t i=0; i < 28; i++)
131 // printf("%f %f %f\n",surveyFidX[i],surveyFidY[i],surveyFidZ[i]);
135 const Int_t kNDIM1 = 3;
136 const Int_t kNDIM2 = 4;
137 const Double_t kIdealOrig[kNDIM1] = {0., 0., 3645.}; // Geant values
142 Double_t ab1[kNDIM1], bc1[kNDIM1], nn1[kNDIM1], plane1[kNDIM1];
143 Double_t nga1[kNDIM2][kNDIM1];
146 Double_t ab2[kNDIM1], bc2[kNDIM1], nn2[kNDIM1], plane2[kNDIM1];
147 Double_t nga2[kNDIM2][kNDIM1];
150 Double_t ab3[kNDIM1], bc3[kNDIM1], nn3[kNDIM1], plane3[kNDIM1];
151 Double_t nga3[kNDIM2][kNDIM1];
154 Double_t ab4[kNDIM1], bc4[kNDIM1], nn4[kNDIM1], plane4[kNDIM1];
155 Double_t nga4[kNDIM2][kNDIM1];
157 // These are the sequence of fiducial marks for sector-1,2,3 & 4 - to be
158 // given by the user.
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};
165 for (i = 0; i < 28; i++)
167 for (Int_t j = 0; j < 4; j++)
169 if (surseq1[j] == i+1)
171 nga1[j][0] = surveyFidX[i];
172 nga1[j][1] = surveyFidY[i];
173 nga1[j][2] = surveyFidZ[i];
175 else if (surseq2[j] == i+1)
177 nga2[j][0] = surveyFidX[i];
178 nga2[j][1] = surveyFidY[i];
179 nga2[j][2] = surveyFidZ[i];
181 else if (surseq3[j] == i+1)
183 nga3[j][0] = surveyFidX[i];
184 nga3[j][1] = surveyFidY[i];
185 nga3[j][2] = surveyFidZ[i];
187 else if (surseq4[j] == i+1)
189 nga4[j][0] = surveyFidX[i];
190 nga4[j][1] = surveyFidY[i];
191 nga4[j][2] = surveyFidZ[i];
197 // First vector on the plane of the fiducial marks.
199 for(i = 0; i < kNDIM1; i++)
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];
207 // Second vector on the plane of the fiducial marks.
209 for(i = 0; i < kNDIM1; i++)
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];
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.
220 // Vector normal to the sector-1.
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];
226 // Vector normal to the sector-2.
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];
232 // Vector normal to the sector-3.
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];
238 // Vector normal to the sector-4.
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];
244 Double_t sizen1 = TMath::Sqrt( nn1[0]*nn1[0] + nn1[1]*nn1[1] +
246 Double_t sizen2 = TMath::Sqrt( nn2[0]*nn2[0] + nn2[1]*nn2[1] +
248 Double_t sizen3 = TMath::Sqrt( nn3[0]*nn3[0] + nn3[1]*nn3[1] +
250 Double_t sizen4 = TMath::Sqrt( nn4[0]*nn4[0] + nn4[1]*nn4[1] +
255 s1 = Double_t(1.)/sizen1 ; //normalization factor
264 s2 = Double_t(1.)/sizen2 ; //normalization factor
273 s3 = Double_t(1.)/sizen3 ; //normalization factor
282 s4 = Double_t(1.)/sizen4 ; //normalization factor
289 // Plane expressed in the hessian normal form, see:
290 // http://mathworld.wolfram.com/HessianNormalForm.html
292 for(i = 0; i < kNDIM1; i++)
294 plane1[i] = nn1[i] * s1;
295 plane2[i] = nn2[i] * s2;
296 plane3[i] = nn3[i] * s3;
297 plane4[i] = nn4[i] * s4;
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;
309 // The center of the sector with fiducial marks as corners,
310 // is the middle point of one diagonal.
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];
317 for(i = 0; i < kNDIM1; i++)
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;
325 const Double_t kZdepth = 23.;//(Zdepth+zoffset) is the distance
326 //from surface to centre of the PMD.
328 for(i = 0; i < kNDIM1; i++)
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;
336 Double_t PMDorig1[kNDIM1];
337 Double_t PMDorig2[kNDIM1];
338 for(i = 0; i < kNDIM1; i++)
340 PMDorig1[i] = (orig1[i] + orig2[i]) * 0.5;
341 PMDorig2[i] = (orig3[i] + orig4[i]) * 0.5;
344 // Get x,y local directions needed to write the global rotation matrix
345 // for the surveyed volume by normalising vectors ab1 and bc1
347 Double_t sx1 = TMath::Sqrt(ab1[0]*ab1[0] + ab1[1]*ab1[1] +
351 for(i = 0; i < 3; i++)
357 Double_t sy1 = TMath::Sqrt(bc1[0]*bc1[0] + bc1[1]*bc1[1] +
361 for(i = 0; i < kNDIM1; i++)
367 Double_t sx2 = TMath::Sqrt(ab2[0]*ab2[0] + ab2[1]*ab2[1] +
371 for(i = 0; i < kNDIM1; i++)
377 Double_t sy2 = TMath::Sqrt(bc2[0]*bc2[0] + bc2[1]*bc2[1] +
381 for(i = 0; i < kNDIM1; i++)
387 Double_t sx3 = TMath::Sqrt(ab3[0]*ab3[0] + ab3[1]*ab3[1] +
391 for(i = 0; i < kNDIM1; i++)
397 Double_t sy3 = TMath::Sqrt(bc3[0]*bc3[0] + bc3[1]*bc3[1] +
401 for(i = 0; i < kNDIM1; i++)
407 Double_t sx4 = TMath::Sqrt(ab4[0]*ab4[0] + ab4[1]*ab4[1] +
411 for(i = 0; i < kNDIM1; i++)
417 Double_t sy4 = TMath::Sqrt(bc4[0]*bc4[0] + bc4[1]*bc4[1] +
421 for(i = 0; i < kNDIM1; i++)
427 // The global matrix for the surveyed volume - ng_1
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]};
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.
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);
442 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-1 **\n";
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]};
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);
456 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-2 **\n";
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]};
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);
470 // cout<<"\n** Global Matrix from surveyed fiducial marks for sector-3 **\n";
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]};
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);
485 //cout<<"\n** Global Matrix from surveyed fiducial marks for sector-4 **\n";
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
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);
497 Printf("The global delta transformation for sector-1");
500 TGeoPhysicalNode* pn3_2 = gGeoManager->MakePhysicalNode("ALIC_1/EPM2_1");
501 TGeoHMatrix gdelta_2 = g3_2->Inverse();
502 gdelta_2.MultiplyLeft(&ng_2);
504 Printf("The global delta transformation for sector-2");
507 TGeoPhysicalNode* pn3_3 = gGeoManager->MakePhysicalNode("ALIC_1/EPM3_1");
508 TGeoHMatrix gdelta_3 = g3_3->Inverse();
509 gdelta_3.MultiplyLeft(&ng_3);
511 Printf("The global delta transformation for sector-3");
514 TGeoPhysicalNode* pn3_4 = gGeoManager->MakePhysicalNode("ALIC_1/EPM4_1");
515 TGeoHMatrix gdelta_4 = g3_4->Inverse();
516 gdelta_4.MultiplyLeft(&ng_4);
518 Printf("The global delta transformation for sector-4");
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);
527 TString strStorage(gSystem->Getenv("TOCDB"));
528 if(strStorage != TString("kTRUE")){
530 TFile f("PMDSurvey.root","RECREATE");
531 if(!f) cerr<<"cannot open file for output\n";
533 f.WriteObject(array,"PMDSurveyObjs ","kSingleKey");
536 // save in CDB storage
537 AliCDBManager* cdb = AliCDBManager::Instance();
538 AliCDBStorage* storage = cdb->GetStorage("local://$ALICE_ROOT/OCDB");
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);