Update master to aliroot
[u/mrichter/AliRoot.git] / PHOS / macros / AlignmentDB / CheckPHOSAlignment.C
CommitLineData
b79eb35f 1//Macros to calculate PHOS alignment matrixes from survay
2//To run select module number (mod=1,2,3,4) and uncomment correcponding
3//pairs of matrixes: survay points in crystal plane and in ALICE global system
4//These data one takes from syrvay results.
5// Note that for Mod.4 there was now survay repers in crystal surface frame.
6// So we choose module with closest ditances between repers (mod2==mod4)
7// and assume that mod4 has similr reper coordinates wrt crytal surface.
8//
9
10TGeoHMatrix *mPHOS[6] ;
11Int_t mod=1 ; //Which module minimize
12
13
14//Module 1, points 500x survey
15 const Float_t srv5000[3]={-217.61,-386.44,-77.98} ;
16 const Float_t srv5001[3]={ -79.32,-439.05,-77.54} ;
17 const Float_t srv5002[3]={ -79.40,-439.29, 80.80} ;
18 const Float_t srv5003[3]={-217.64,-386.49, 80.04} ;
19
20//Module 1, points 500x wrt crystals (copied from Mod3)
21 const Float_t p5000[3]={ -3.08,28.97,-16.43} ;
22 const Float_t p5001[3]={145.94,28.90,-16.03} ;
23 const Float_t p5002[3]={145.58,29.01,142.17} ;
24 const Float_t p5003[3]={ -3.58,28.88,141.88} ;
25
26
27
28
29
30/*
31//Module 2, points 500x survey
32 const Float_t srv5000[3]={-74.67,-429.97,-90.95} ;
33 const Float_t srv5001[3]={ 72.39,-430.96,-91.25} ;
34 const Float_t srv5002[3]={ 64.43,-431.06, 92.01} ;
35 const Float_t srv5003[3]={-70.56,-430.27, 92.43} ;
36*/
37/*
38//Module 3, points 500x survey
39 const Float_t srv5000[3]={ 77.95,-430.77,-78.85} ;
40 const Float_t srv5001[3]={217.78,-380.47,-78.76} ;
41 const Float_t srv5002[3]={217.86,-380.25, 79.28} ;
42 const Float_t srv5003[3]={ 77.57,-430.93, 79.58} ;
43*/
44/*
45//Module 4, points 500x survey
46 const Float_t srv5000[3]={222.58,-377.81,-78.96} ;
47 const Float_t srv5001[3]={333.99,-284.63,-78.21} ;
48 const Float_t srv5002[3]={335.38,-283.62, 80.39} ;
49 const Float_t srv5003[3]={220.94,-379.38, 79.89} ;
50*/
51
52
53/*
54 //Module 2 reference points on crystal surface (in cm)
55 const Float_t p5000[3]={0.0000,0.0000, 0.00} ;
56 const Float_t p5001[3]={0.0000,0.0000,122.62} ;
57 const Float_t p5002[3]={142.42,0.0000, 1.95} ;
58 const Float_t p5003[3]={142.49,0.0001,124.69} ;
59
60 //Module 2 survey: position of these points in global ALICE system
61 const Float_t srv5000[3]={-71.06,-459.00,-61.82} ;
62 const Float_t srv5001[3]={-71.16,-459.15, 61.06} ;
63 const Float_t srv5002[3]={ 71.65,-459.99,-59.76} ;
64 const Float_t srv5003[3]={ 71.63,-460.13, 63.24} ;
65*/
66/*
67 //Module 3 reference points on crystal surface
68 const Float_t p5000[3]={0.0000,0.0000,0.0000 } ;
69 const Float_t p5001[3]={0.0000,0.0000,122.63 } ;
70 const Float_t p5002[3]={142.53,0.0000, 2.37 } ;
71 const Float_t p5003[3]={142.28, -0.11,125.06 } ;
72
73 //Module 3 survey
74 const Float_t srv5000[3]={ 90.51,-457.01,-62.35 } ;
75 const Float_t srv5001[3]={ 90.77,-456.86, 60.19 } ;
76 const Float_t srv5002[3]={224.51,-408.73,-60.33 } ;
77 const Float_t srv5003[3]={224.57,-408.77, 62.27 } ;
78
79*/
80/*
81 //Module 4 reference points on crystal surface
82 const Float_t p5000[3]={0.0000,0.0000,0.0000 } ;
83 const Float_t p5001[3]={0.0000,0.0000,122.98 } ;
84 const Float_t p5002[3]={142.39,0.0000, 2.05 } ;
85 const Float_t p5003[3]={142.41,0.0019,124.82 } ;
86
87 //Module 4 survey
88 const Float_t srv5000[3]={244.57,-399.74,-62.67 } ;
89 const Float_t srv5001[3]={244.10,-400.79, 60.38 } ;
90 const Float_t srv5002[3]={353.95,-308.47,-59.42 } ;
91 const Float_t srv5003[3]={353.38,-309.36, 63.42 } ;
92*/
93
94
95
96
97
98Double_t DistV(TVector3 &a, Float_t * b){
99
100 return sqrt((a[0]-b[0])*(a[0]-b[0])+ (a[1]-b[1])*(a[1]-b[1])+ (a[2]-b[2])*(a[2]-b[2]) ) ;
101
102}
103
104Double_t Dist(Float_t *a, Float_t *b){
105
106 return sqrt((a[0]-b[0])*(a[0]-b[0])+ (a[1]-b[1])*(a[1]-b[1])+ (a[2]-b[2])*(a[2]-b[2]) ) ;
107
108}
109void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
110{
111 //Function used in Minuit minimization:
112 //Minimize distances between measred reper coordinates
113 //and calculate after PHOS rotaion/shifts.
114 //Uses dedicated function AliPHOSGeoUtils::TestSurvey()
115
116 Double_t chisq = 0;
117
118 //Create current rotation matrix
119 TGeoRotation r1;
120 r1.SetAngles(par[0],par[1],par[2]) ;
121 mPHOS[mod]->SetRotation(r1.GetRotationMatrix()) ;
122 mPHOS[mod]->SetDx(par[3]) ;
123 mPHOS[mod]->SetDy(par[4]) ;
124 mPHOS[mod]->SetDz(par[5]) ;
125
126 AliPHOSGeoUtils * geom = new AliPHOSGeoUtils("PHOS") ;
127 geom->SetMisalMatrix(mPHOS[mod],mod-1) ;
128
129 TVector3 globaPos ;
130 //Compare ditances for all 4 points
131 geom->TestSurvey(mod,p5000,globaPos) ;
132 chisq+=DistV(globaPos,srv5000) ;
133
134 geom->TestSurvey(mod,p5001,globaPos) ;
135 chisq+=DistV(globaPos,srv5001) ;
136
137 geom->TestSurvey(mod,p5002,globaPos) ;
138 chisq+=DistV(globaPos,srv5002) ;
139
140 geom->TestSurvey(mod,p5003,globaPos) ;
141 chisq+=DistV(globaPos,srv5003) ;
142
143 delete geom ;
144 f = chisq;
145}
146
147
148void CheckPHOSAlignment(){
149 //Checks if PHOS alignment agree with results of photogrammetry
150
151
152 mPHOS[mod] = new TGeoHMatrix() ;
153
154 //Initial conditions for minimization, close to ideal PHOS positions
155 TGeoRotation r1;
156 r1.SetAngles(0,90+(mod-2)*20,0) ;
157 r1.Print() ;
158
159 mPHOS[mod]->SetRotation(r1.GetRotationMatrix()) ;
160 mPHOS[mod]->SetDx(-480.*TMath::Cos((210+mod*20)*TMath::DegToRad())) ;
161 mPHOS[mod]->SetDy(-480.*TMath::Sin((210+mod*20)*TMath::DegToRad())) ;
162 mPHOS[mod]->SetDz(0.) ;
163
164
165
166 //Each module has 4 reference points
167 //Define them w.r.t. center of top right crystal (see more details in EDH comenents
168 // 992344, 1002387, 1012391.
169
170 //Print distances between repers in current module
171//printf("%6.2f,\t %6.2f,\t %6.2f,\t %6.2f\n",Dist(srv5000,srv5001),Dist(srv5000,srv5002),Dist(srv5000,srv5003),Dist(srv5001,srv5003)) ;
172//return ;
173
174
175
176 //Module 3
177 const Float_t m3p5000[3]={-30.76, 191.15, -164.32} ;
178 const Float_t m3p5001[3]={1459.44, 190.47, -160.25} ;
179 const Float_t m3p5002[3]={1455.82, 191.57, 1421.70} ;
180 const Float_t m3p5003[3]={-35.83, 190.32, 1418.78} ;
181
182 //Module 3 survay (in meters!)
183 const Float_t m3srv5000[3]={0.7795,-4.3077,-0.7885} ;
184 const Float_t m3srv5001[3]={2.1778,-3.8047,-0.7876} ;
185 const Float_t m3srv5002[3]={2.1786,-3.8025,0.7928} ;
186 const Float_t m3srv5003[3]={0.7757,-4.3093,0.7958} ;
187
188
189 //Module 4
190 const Float_t m4p5000[3]={-34.59, 210.66, -165.23} ;
191 const Float_t m4p5001[3]={1456.73, 210.71, -164.91} ;
192 const Float_t m4p5002[3]={1454.76, 213.62, 1418.01} ;
193 const Float_t m4p5003[3]={-36.60, 210.14, 1419.84} ;
194
195
196 const Float_t m4srv5000[3]={2.2258,-3.7781,-0.7896} ;
197 const Float_t m4srv5001[3]={3.3399,-2.8463,-0.7821} ;
198 const Float_t m4srv5002[3]={3.3538,-2.8362, 0.8039} ;
199 const Float_t m4srv5003[3]={2.2094,-3.7938, 0.7989} ;
200
201
202
203
204 // prepare minuit
205 const Double_t nPar = 6; // number of fit parameters
206 TMinuit* gMinuit = new TMinuit(nPar);
207 gMinuit->SetFCN(fcn);
208 gMinuit->SetPrintLevel(1); // -1 quiet, 0 normal, 1 verbose
209
210 Double_t arglist[10];
211 Int_t ierflg = 0;
212
213 // start minimization
214 Double_t chi2 = 0;
215
216 ierflg = 0;
217 arglist[0] = 1;
218 gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
219 gMinuit->SetMaxIterations(5000);
220
221 // parameters: parameter no., name, start value, step size, parameter range min., parameter range max., status
222 gMinuit->mnparm(0,"A1",0., 1., 0, 0, ierflg);
223 gMinuit->mnparm(1,"A2",90., 1., 0, 0, ierflg);
224 gMinuit->mnparm(2,"A3",0., 1., 0, 0, ierflg);
225 gMinuit->mnparm(3,"dx",0., 1., 0, 0, ierflg);
226 gMinuit->mnparm(4,"dy",-460., 1., 0, 0, ierflg);
227 gMinuit->mnparm(5,"dz",0., 1., 0, 0, ierflg);
228
229 // Now ready for minimization step
230 arglist[0] = 5000;
231 arglist[1] = 1.;
232 gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
233
234 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("PHOS_mod1234") ;
235 geom->SetMisalMatrix(mPHOS[mod],mod-1) ;
236
237
238
239 TVector3 globaPos ;
240
241
242 geom->TestSurvey(mod,p5000,globaPos) ;
243 printf("P5000: (%f, %f, %f) \n",globaPos.X(),globaPos.Y(),globaPos.Z()) ;
244 printf(" srv: (%f, %f, %f) \n",srv5000[0],srv5000[1],srv5000[2]) ;
245 geom->TestSurvey(mod,p5001,globaPos) ;
246 printf("P5001: (%f, %f, %f) \n",globaPos.X(),globaPos.Y(),globaPos.Z()) ;
247 printf(" srv: (%f, %f, %f) \n",srv5001[0],srv5001[1],srv5001[2]) ;
248 geom->TestSurvey(mod,p5002,globaPos) ;
249 printf("P5002: (%f, %f, %f) \n",globaPos.X(),globaPos.Y(),globaPos.Z()) ;
250 printf(" srv: (%f, %f, %f) \n",srv5002[0],srv5002[1],srv5002[2]) ;
251 geom->TestSurvey(mod,p5003,globaPos) ;
252 printf("P5003: (%f, %f, %f) \n",globaPos.X(),globaPos.Y(),globaPos.Z()) ;
253 printf(" srv: (%f, %f, %f) \n",srv5003[0],srv5003[1],srv5003[2]) ;
254
255 mPHOS[mod]->Print() ;
256}
257
258
259