Updated macros for PHOS alignment calculation
[u/mrichter/AliRoot.git] / PHOS / macros / AlignmentDB / CheckPHOSAlignment.C
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
10 TGeoHMatrix *mPHOS[6] ;
11 Int_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
98 Double_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
104 Double_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 }
109 void 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
148 void 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