]>
Commit | Line | Data |
---|---|---|
7e2c89ed | 1 | void SurveyToAlignmentExample(){ |
3b583936 | 2 | // Macro to show an example of conversion of survey data into alignment |
3 | // data. The position of four fiducial marks, sticked above one surface | |
4 | // of a box is converted into the global position of the box. | |
5 | // | |
6 | gSystem->Load("libGeom"); | |
7 | TGeoManager *mgr = new TGeoManager("Geom","survey to alignment toy"); | |
8 | TGeoMedium *medium = 0; | |
9 | TGeoVolume *top = mgr->MakeBox("TOP",medium,250,250,250); | |
10 | mgr->SetTopVolume(top); | |
11 | // make shape components | |
12 | // ******** red outermost box *************** | |
13 | TGeoBBox *sbox0 = new TGeoBBox(200,200,50); | |
14 | TGeoVolume* box0 = new TGeoVolume("B0",sbox0); | |
15 | box0->SetVisDaughters(); | |
16 | box0->SetLineColor(2); //red | |
17 | top->AddNode(box0,1); | |
18 | // ******** green middle box *************** | |
19 | TGeoBBox *sbox1 = new TGeoBBox(180,180,40); | |
20 | TGeoVolume* box1 = new TGeoVolume("B1",sbox1); | |
21 | box1->SetLineColor(3);//green | |
22 | TGeoTranslation* tr = new TGeoTranslation("tr",10,0,0); | |
23 | box0->AddNode(box1,1,tr); | |
24 | // ******** bleu inner box *************** | |
25 | TGeoBBox *sbox2 = new TGeoBBox(160,160,30); | |
26 | TGeoVolume* box2 = new TGeoVolume("B2",sbox2); | |
27 | box2->SetLineColor(4);//bleu | |
28 | box1->AddNode(box2,1,tr); | |
29 | // ******** violet innermost box *************** | |
30 | Double_t zsize = 20.; | |
31 | TGeoBBox *sbox3 = new TGeoBBox(140,140,zsize); | |
32 | TGeoVolume* box3 = new TGeoVolume("B3",sbox3); | |
33 | box3->SetLineColor(6);//violet | |
34 | box2->AddNode(box3,1,tr); | |
35 | ||
36 | // Four fiducial marks on the box3, expressed in local coordinates | |
37 | // We imagine they are at 2mm above the upper surface of the volume | |
38 | // at the corners of a square of 200 cm side | |
39 | const Double_t xside = 100; | |
40 | const Double_t yside = 100; | |
41 | const Double_t zoffset = 0.2; | |
42 | const Double_t zdepth = zsize+zoffset; | |
43 | Double_t A[3]={-xside,-yside,zdepth}; | |
44 | Double_t B[3]={xside,-yside,zdepth}; | |
45 | Double_t C[3]={xside,yside,zdepth}; | |
46 | Double_t D[3]={-xside,yside,zdepth}; | |
47 | ||
48 | TGeoBBox *fmbox = new TGeoBBox(1,1,1); | |
49 | TGeoVolume* fm = new TGeoVolume("FM",fmbox); | |
50 | fm->SetLineColor(7);//color | |
51 | TGeoTranslation* Atr = new TGeoTranslation("Atr",-xside,-yside,zdepth); | |
52 | TGeoTranslation* Btr = new TGeoTranslation("Btr",xside,-yside,zdepth); | |
53 | TGeoTranslation* Ctr = new TGeoTranslation("Ctr",xside,yside,zdepth); | |
54 | TGeoTranslation* Dtr = new TGeoTranslation("Dtr",-xside,yside,zdepth); | |
55 | ||
56 | box3->AddNode(fm,1,Atr); | |
57 | box3->AddNode(fm,2,Btr); | |
58 | box3->AddNode(fm,3,Ctr); | |
59 | box3->AddNode(fm,4,Dtr); | |
60 | ||
61 | // ^ local y | |
62 | // | | |
63 | // D-------------|-------------C | |
64 | // | | | | |
65 | // | | | | |
66 | // | | | | |
67 | // | | | | |
68 | // | | | | |
69 | // | | | | |
70 | // ------------------|------------------> local x | |
71 | // | | | | |
72 | // | | | | |
73 | // | | | | |
74 | // | | | | |
75 | // | | | | |
76 | // | | | | |
77 | // A-------------|-------------B | |
78 | // | |
79 | // local z exiting the plane of the screen | |
80 | ||
81 | mgr->CloseGeometry(); | |
82 | mgr->GetTopVolume()->Draw(); | |
83 | mgr->SetVisOption(0); | |
84 | mgr->SetVisLevel(6); | |
85 | ||
86 | Int_t i; | |
87 | // ************* get ideal global matrix ******************* | |
88 | mgr->cd("TOP_1/B0_1/B1_1/B2_1/B3_1"); | |
89 | TGeoHMatrix g3 = *mgr->GetCurrentMatrix(); // !!don't declare g3 | |
90 | // as a pointer to mgr->GetCurrentMatrix(), mgr->cd("...") | |
91 | // would eventually change the content pointed by g3 behind your back | |
92 | ||
93 | // ************* get ideal local matrix ******************* | |
94 | TGeoNode* n3 = mgr->GetCurrentNode(); | |
95 | TGeoMatrix* l3 = n3->GetMatrix(); | |
96 | ||
97 | Double_t gA[3], gB[3], gC[3], gD[3]; // point coordinates in the global RS | |
98 | g3.LocalToMaster(A,gA); | |
99 | g3.LocalToMaster(B,gB); | |
100 | g3.LocalToMaster(C,gC); | |
101 | g3.LocalToMaster(D,gD); | |
102 | cout<<endl<<"Ideal fiducial marks coordinates in the global RS:\n"<< | |
103 | "A "<<gA[0]<<" "<<gA[1]<<" "<<gA[2]<<" "<<endl<< | |
104 | "B "<<gB[0]<<" "<<gB[1]<<" "<<gB[2]<<" "<<endl<< | |
105 | "C "<<gC[0]<<" "<<gC[1]<<" "<<gC[2]<<" "<<endl<< | |
106 | "D "<<gD[0]<<" "<<gD[1]<<" "<<gD[2]<<" "<<endl; | |
107 | ||
108 | // We apply a delta transformation to the surveyed vol box3 to represent | |
109 | // its real position, given below by ng3 nl3, which differs from its | |
110 | // ideal position saved above in g3 and l3 | |
111 | TGeoPhysicalNode* pn3 = mgr->MakePhysicalNode("TOP_1/B0_1/B1_1/B2_1/B3_1"); | |
112 | Double_t dphi = 3; // tilt by 3 degrees around z | |
113 | Double_t dz = 5; // shift by 5 cm along z | |
114 | TGeoRotation* rrot = new TGeoRotation("rot",dphi,0.,0.); | |
115 | TGeoCombiTrans localdelta = *(new TGeoCombiTrans(0.,0.,dz, rrot)); | |
116 | // new local matrix, representing real position | |
117 | TGeoHMatrix nlocal = *l3 * localdelta; | |
118 | TGeoHMatrix* nl3 = new TGeoHMatrix(nlocal); | |
119 | pn3->Align(nl3); | |
120 | ||
121 | //Let's get the global matrix for later comparison | |
122 | TGeoHMatrix* ng3 = pn3->GetMatrix(); //"real" global matrix, what survey sees | |
123 | printf("\n\n************ real global matrix **************\n"); | |
124 | ng3->Print(); | |
125 | Double_t ngA[3], ngB[3], ngC[3], ngD[3]; | |
126 | ng3->LocalToMaster(A,ngA); | |
127 | ng3->LocalToMaster(B,ngB); | |
128 | ng3->LocalToMaster(C,ngC); | |
129 | ng3->LocalToMaster(D,ngD); | |
130 | ||
131 | cout<<endl<<"Fiducial marks coordinates in the global RS given by survey:\n"<< | |
132 | "A "<<ngA[0]<<" "<<ngA[1]<<" "<<ngA[2]<<" "<<endl<< | |
133 | "B "<<ngB[0]<<" "<<ngB[1]<<" "<<ngB[2]<<" "<<endl<< | |
134 | "C "<<ngC[0]<<" "<<ngC[1]<<" "<<ngC[2]<<" "<<endl<< | |
135 | "D "<<ngD[0]<<" "<<ngD[1]<<" "<<ngD[2]<<" "<<endl; | |
136 | ||
137 | ||
138 | // From the new fiducial marks coordinates derive back the | |
139 | // new global position of the surveyed volume | |
140 | //*** What follows is the actual survey-to-alignment procedure which assumes, | |
141 | //*** as is the case of the present example, 4 fiducial marks | |
142 | //*** at the corners of a square lying on a plane parallel to a surface | |
143 | //*** of the surveyed box at a certain offset and with | |
144 | //*** x and y sides parallel to the box's x and y axes. | |
145 | //*** If the code below is placed in a separate class or method, it needs | |
146 | //*** as input the four points and the offset from the origin (zdepth) | |
147 | //*** The algorithm can be easily modified for different placement | |
148 | //*** and/or cardinality of the fiducial marks. | |
149 | ||
150 | Double_t ab[3], bc[3], n[3]; | |
151 | Double_t plane[4], s; | |
152 | ||
153 | // first vector on the plane of the fiducial marks | |
154 | for(i=0;i<3;i++){ | |
155 | ab[i] = ngB[i] - ngA[i]; | |
156 | } | |
157 | ||
158 | // second vector on the plane of the fiducial marks | |
159 | for(i=0;i<3;i++){ | |
160 | bc[i] = ngC[i] - ngB[i]; | |
161 | } | |
162 | ||
163 | // vector normal to the plane of the fiducial marks obtained | |
164 | // as cross product of the two vectors on the plane d0^d1 | |
165 | n[0] = ab[1] * bc[2] - ab[2] * bc[1]; | |
166 | n[1] = ab[2] * bc[0] - ab[0] * bc[2]; | |
167 | n[2] = ab[0] * bc[1] - ab[1] * bc[0]; | |
168 | ||
169 | Double_t sizen = TMath::Sqrt( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] ); | |
170 | if(sizen>1.e-8){ | |
171 | s = Double_t(1.)/sizen ; //normalization factor | |
172 | }else{ | |
173 | return 0; | |
174 | } | |
175 | ||
176 | // plane expressed in the hessian normal form, see: | |
177 | // http://mathworld.wolfram.com/HessianNormalForm.html | |
178 | // the first three are the coordinates of the orthonormal vector | |
179 | // the fourth coordinate is equal to the distance from the origin | |
180 | for(i=0;i<3;i++){ | |
181 | plane[i] = n[i] * s; | |
182 | } | |
183 | plane[3] = -( plane[0] * ngA[0] + plane[1] * ngA[1] + plane[2] * ngA[2] ); | |
184 | cout<<plane[0]<<" "<<plane[1]<<" "<<plane[2]<<" "<<plane[3]<<" "<<endl; | |
185 | ||
186 | // The center of the square with fiducial marks as corners | |
187 | // as the middle point of one diagonal - md | |
188 | // Used below to get the center - orig - of the surveyed box | |
189 | Double_t orig[3], md[3]; | |
190 | for(i=0;i<3;i++){ | |
191 | md[i] = (ngA[i] + ngC[i]) * 0.5; | |
192 | } | |
193 | ||
194 | // The center of the box | |
195 | for(i=0;i<3;i++){ | |
196 | orig[i] = md[i] - plane[i]*zdepth; | |
197 | } | |
198 | orig[1] = md[1] - plane[1]*zdepth; | |
199 | orig[2] = md[2] - plane[2]*zdepth; | |
200 | cout<<endl<<"The origin of the box: "<<orig[0]<<" "<<orig[1]<<" "<<orig[2]<<endl; | |
201 | ||
202 | // get x,y local directions needed to write the global rotation matrix | |
203 | // for the surveyed volume by normalising vectors ab and bc | |
204 | Double_t sx = TMath::Sqrt(ab[0]*ab[0] + ab[1]*ab[1] + ab[2]*ab[2]); | |
205 | if(sx>1.e-8){ | |
206 | for(i=0;i<3;i++){ | |
207 | ab[i] /= sx; | |
208 | } | |
209 | cout<<endl<<"x "<<ab[0]<<" "<<ab[1]<<" "<<ab[2]<<endl; | |
210 | } | |
211 | Double_t sy = TMath::Sqrt(bc[0]*bc[0] + bc[1]*bc[1] + bc[2]*bc[2]); | |
212 | if(sy>1.e-8){ | |
213 | for(i=0;i<3;i++){ | |
214 | bc[i] /= sy; | |
215 | } | |
216 | cout<<endl<<"y "<<bc[0]<<" "<<bc[1]<<" "<<bc[2]<<endl; | |
217 | } | |
218 | ||
219 | ||
220 | // the global matrix for the surveyed volume - ng | |
221 | Double_t rot[9] = {ab[0],bc[0],plane[0],ab[1],bc[1],plane[1],ab[2],bc[2],plane[2]}; | |
222 | TGeoHMatrix ng; | |
223 | ng.SetTranslation(orig); | |
224 | ng.SetRotation(rot); | |
225 | ||
226 | cout<<"\n********* global matrix inferred from surveyed fiducial marks ***********\n"; | |
227 | ng.Print(); | |
228 | ||
229 | // // To produce the alignment object for the given volume you would | |
230 | // // then do something like this: | |
231 | // // Calculate the global delta transformation as ng * g3-1 | |
232 | // TGeoHMatrix gdelta = g3->Inverse(); //now equal to the inverse of g3 | |
233 | // gdelta.MultiplyLeft(&ng); | |
234 | // Int_t index = 0; | |
235 | // // if the volume is in the look-up table use something like this instead: | |
ae079791 | 236 | // // AliGeomManager::LayerToVolUID(AliGeomManager::kTOF,i); |
3b583936 | 237 | // AliAlignObjMatrix* mobj = new AliAlignObjMatrix("symname",index,gdelta,kTRUE); |
238 | ||
239 | } |