Resolving all symbols in the library
[u/mrichter/AliRoot.git] / MUON / MUONSurveyUtil.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /// \ingroup macros
19 /// \file MUONSurveyUtil.C
20 /// \brief Utility macro for survey data to alignment transformation.
21 ///  
22 /// Macro contains various functions to calculate misalignement parameters
23 /// from survey data and designed positions of survey targets.
24 /// Macro also includes a method to get the new AliMUONGeometryTransformer.
25 /// It is intended to be loaded by chamber specific macros.
26 /// 
27 /// \author Javier Castillo
28
29 #if !defined(__CINT__) || defined(__MAKECINT__)
30
31 #include "AliMUONGeometryTransformer.h"
32 #include "AliMUONGeometryMisAligner.h"
33 #include "AliMUONGeometryModuleTransformer.h"
34 #include "AliMUONGeometryDetElement.h"
35 #include "AliMUONGeometryBuilder.h"
36 #include "AliMpExMap.h"
37 #include "AliMpExMapIterator.h"
38 #include "AliGeomManager.h"
39 #include "AliCDBManager.h"
40 #include "AliCDBMetaData.h"
41 #include "AliCDBId.h"
42
43 #include <TGeoManager.h>
44 #include <TClonesArray.h>
45 #include <TMath.h>
46 #include <TString.h>
47 #include <Riostream.h>
48
49 #include <fstream>
50
51 #endif
52
53 static int fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
54
55 Bool_t MatrixToAngles(const Double_t *rot, Double_t *angles)
56 {
57   // Calculates the Euler angles in "x y z" notation
58   // using the rotation matrix
59   // Returns false in case the rotation angles can not be
60
61   // extracted from the matrix
62   //
63   if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
64     printf("Failed to extract roll-pitch-yall angles!");
65     return kFALSE;
66   }
67   //   Double_t raddeg = TMath::RadToDeg();
68   angles[0]=TMath::ATan2(-rot[5],rot[8]);
69   angles[1]=TMath::ASin(rot[2]);
70   angles[2]=TMath::ATan2(-rot[1],rot[0]);
71   return kTRUE;
72 }
73
74 Double_t eqPlane(Double_t *x, Double_t *par){
75   return (-par[0]*x[0] -par[1]*x[1] -par[2]); 
76 }
77
78 Double_t xpCenter(Double_t *x, Double_t *par){
79
80   Double_t lCos2Tht = TMath::Cos(2*par[6]);
81   Double_t lSinTht = TMath::Sin(par[6]);
82
83   Double_t inSqrt = TMath::Abs((par[0] - par[3])*(par[0] - par[3]) 
84                                -2*(x[0] -x[1])*(x[0] -x[1]) 
85                                +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
86                                +((par[0] - par[3])*(par[0] - par[3]) 
87                                  +(par[1] - par[4])*(par[1] - par[4]) 
88                                  +(par[2] - par[5])*(par[2] - par[5]))*lCos2Tht 
89                                +4*(x[0] - x[1])*(par[2] - par[5])*lSinTht);
90
91   Double_t xD = ((2*(par[0]*par[0]*x[1] 
92                      -par[0]*par[3]*(x[0] + x[1]) 
93                      +x[1]*par[1]*(par[1] - par[4]) 
94                      +x[0]*(par[3]*par[3] - par[1]*par[4] + par[4]*par[4])) 
95                   -2*(par[3]*par[3]*par[2] 
96                       +par[0]*par[0]*par[5] 
97                       -par[0]*par[3]*(par[2] + par[5]) 
98                       +(par[1] - par[4])*(-par[4]*par[2] +par[1]*par[5]))*lSinTht 
99                   +TMath::Sqrt(2)*(-par[3]*par[1] + par[0]*par[4])
100                   *TMath::Sqrt(inSqrt))
101                  /(2*((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4]))));
102
103   return xD;
104 }
105
106 Double_t xnCenter(Double_t *x, Double_t *par){
107
108   Double_t lCos2Tht = TMath::Cos(2*par[6]);
109   Double_t lSinTht = TMath::Sin(par[6]);
110
111   Double_t inSqrt = TMath::Abs((par[0] - par[3])*(par[0] - par[3]) 
112                                -2*(x[0] - x[1])*(x[0] - x[1]) 
113                                +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
114                                +((par[0] - par[3])*(par[0] - par[3]) 
115                                  +(par[1] - par[4])*(par[1] - par[4]) 
116                                  +(par[2] - par[5])*(par[2] - par[5]))*lCos2Tht
117                                +4*(x[0] - x[1])*(par[2] - par[5])*lSinTht);
118
119   Double_t xD = ((2*(par[0]*par[0]*x[1] 
120                      -par[0]*par[3]*(x[0] + x[1]) 
121                      +x[1]*par[1]*(par[1] - par[4]) 
122                      +x[0]*(par[3]*par[3] - par[1]*par[4] + par[4]*par[4])) 
123                   -2*(par[3]*par[3]*par[2] + par[0]*par[0]*par[5] 
124                       -par[0]*par[3]*(par[2] + par[5]) 
125                       +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lSinTht 
126                   +TMath::Sqrt(2)*(par[3]*par[1] - par[0]*par[4])
127                   *TMath::Sqrt(inSqrt))
128                  /(2*((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4]))));
129
130   return xD;
131 }
132
133 Double_t phixpn(Double_t *x, Double_t *par){
134
135   Double_t inSqrt = TMath::Abs(((par[0] - par[3])*(par[0] - par[3]) 
136                                 -2*(x[0] - x[1])*(x[0] - x[1]) 
137                                 +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
138                                 +(+(par[0] - par[3])*(par[0] - par[3]) 
139                                   +(par[1] - par[4])*(par[1] - par[4]) 
140                                   +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6]) 
141                                 +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6])));
142   
143   Double_t phix = ((+2*(par[0] - par[3])*(x[0] - x[1]) 
144                     -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6]) 
145                     +TMath::Sqrt(2)*(par[1] - par[4])
146                     *TMath::Sqrt(inSqrt))
147                    /(2*(+(par[0] - par[3])*(par[0] - par[3]) 
148                         +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
149
150   phix = -TMath::ACos(phix);
151
152   return phix;
153 }
154
155 Double_t phixpp(Double_t *x, Double_t *par){
156
157   Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3]) 
158                                -2*(x[0] - x[1])*(x[0] - x[1]) 
159                                +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
160                                +(+(par[0] - par[3])*(par[0] - par[3]) 
161                                  +(par[1] - par[4])*(par[1] - par[4]) 
162                                  +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6]) 
163                                +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
164
165   Double_t phix = ((+2*(par[0] - par[3])*(x[0] - x[1]) 
166                     -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6]) 
167                     +TMath::Sqrt(2)*(par[1] - par[4])
168                     *TMath::Sqrt(inSqrt))
169                    /(2*(+(par[0] - par[3])*(par[0] - par[3]) 
170                         +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
171
172   phix = TMath::ACos(phix);
173
174   return phix;
175 }
176
177 Double_t phixnn(Double_t *x, Double_t *par){
178
179   Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3]) 
180                                -2*(x[0] - x[1])*(x[0] - x[1]) 
181                                +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
182                                +(+(par[0] - par[3])*(par[0] - par[3]) 
183                                  +(par[1] - par[4])*(par[1] - par[4]) 
184                                  +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6]) 
185                                + 4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
186   
187   Double_t phix = (+(+2*(par[0] - par[3])*(x[0] - x[1]) 
188                      -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6]) 
189                      +TMath::Sqrt(2)*(-par[1] + par[4])
190                      *TMath::Sqrt(inSqrt))
191                    /(2*(+(par[0] - par[3])*(par[0] - par[3]) 
192                         +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
193
194   phix = -TMath::ACos(phix);
195
196   return phix;
197 }
198
199 Double_t phixnp(Double_t *x, Double_t *par){
200
201   Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3]) 
202                                -2*(x[0] - x[1])*(x[0] - x[1]) 
203                                +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5]) 
204                                +(+(par[0] - par[3])*(par[0] - par[3]) 
205                                  +(par[1] - par[4])*(par[1] - par[4]) 
206                                  +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6]) 
207                                +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
208
209   Double_t phix = (+(+2*(par[0] - par[3])*(x[0] - x[1]) 
210                      -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6]) 
211                      +TMath::Sqrt(2)*(-par[1] + par[4])
212                      *TMath::Sqrt(inSqrt))
213                    /(2*(+(par[0] - par[3])*(par[0] - par[3]) 
214                         +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
215
216   phix = TMath::ACos(phix);
217
218   return phix;
219 }
220
221 Double_t ypCenter(Double_t *x, Double_t *par){
222   // par : x1l, y1l, z1l, x2l, y2l, z2, lpsi, tht,
223
224   Double_t lCosPsi = TMath::Cos(par[6]);
225   Double_t lSinPsi = TMath::Sin(par[6]);
226   Double_t lCosTht = TMath::Cos(par[7]);
227   Double_t lSinTht = TMath::Sin(par[7]);
228
229   Double_t yD = ((1./((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4])))
230                  *(+par[3]*par[3]*x[0] 
231                    +par[0]*par[0]*x[1] 
232                    -par[0]*par[3]*(x[0] + x[1]) 
233                    +(par[1] - par[4])*(-x[0]*par[4] + par[1]*x[1]) 
234                    +(par[3]*par[3]*par[2] 
235                      +par[0]*par[0]*par[5] 
236                      -par[0]*par[3]*(par[2] + par[5]) 
237                      +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lCosTht*lSinPsi 
238                    +(-par[3]*par[1] + par[0]*par[4])
239                    *TMath::Sqrt(-(x[0] - x[1] 
240                                   +(par[2] - par[5])*lCosTht*lSinPsi)
241                                 *(x[0] - x[1] 
242                                   +(par[2] - par[5])*lCosTht*lSinPsi) 
243                                 + ((par[0] - par[3])*(par[0] - par[3]) 
244                                    +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
245                                                                           +lSinPsi*lSinPsi*lSinTht*lSinTht))));
246
247   return yD;  
248 }
249
250 Double_t phiypn(Double_t *x, Double_t *par){
251
252   Double_t lCosPsi = TMath::Cos(par[6]);
253   Double_t lSinPsi = TMath::Sin(par[6]);
254   Double_t lCosTht = TMath::Cos(par[7]);
255   Double_t lSinTht = TMath::Sin(par[7]);
256
257   Double_t phiy = ((lCosPsi*((par[1] - par[4])*(x[0] - x[1]) 
258                              +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi 
259                              +(-par[0] + par[3])
260                              *TMath::Sqrt(-(x[0] - x[1] 
261                                             +(par[2] - par[5])*lCosTht*lSinPsi)
262                                           *(x[0] - x[1] 
263                                             +(par[2] - par[5])*lCosTht*lSinPsi) 
264                                           +(+(par[0] - par[3])*(par[0] - par[3]) 
265                                             +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
266                                                                                    +lSinPsi*lSinPsi*lSinTht*lSinTht))) 
267                     +lSinPsi*lSinTht*((par[0] - par[3])*(x[0] - x[1]) 
268                                       +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi 
269                                       +(par[1] - par[4])
270                                       *TMath::Sqrt(-(x[0] - x[1] 
271                                                      +(par[2] - par[5])*lCosTht*lSinPsi)
272                                                    *(x[0] - x[1] 
273                                                      +(par[2] - par[5])*lCosTht*lSinPsi) 
274                                                    + ((par[0] - par[3])*(par[0] - par[3]) 
275                                                       +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
276                                                                                              +lSinPsi*lSinPsi*lSinTht*lSinTht))))
277                    /((+(par[0] - par[3])*(par[0] - par[3]) 
278                       +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
279                                                              +lSinPsi*lSinPsi*lSinTht*lSinTht)));
280   
281   phiy = -TMath::ACos(phiy);
282
283
284   return phiy;
285 }
286
287 Double_t phiypp(Double_t *x, Double_t *par){
288
289   Double_t lCosPsi = TMath::Cos(par[6]);
290   Double_t lSinPsi = TMath::Sin(par[6]);
291   Double_t lCosTht = TMath::Cos(par[7]);
292   Double_t lSinTht = TMath::Sin(par[7]);
293
294   Double_t phiy = ((lCosPsi*((par[1] - par[4])*(x[0] - x[1]) 
295                              +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi 
296                              +(-par[0] + par[3])
297                              *TMath::Sqrt(-(x[0] - x[1] 
298                                             +(par[2] - par[5])*lCosTht*lSinPsi)
299                                           *(x[0] - x[1] 
300                                             +(par[2] - par[5])*lCosTht*lSinPsi) 
301                                           +((par[0] - par[3])*(par[0] - par[3]) 
302                                             +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
303                                                                                    +lSinPsi*lSinPsi*lSinTht*lSinTht))) 
304                     +lSinPsi*lSinTht*((par[0] - par[3])*(x[0] - x[1]) 
305                                       +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi 
306                                       +(par[1] - par[4])*TMath::Sqrt(-(x[0] - x[1] 
307                                                                        +(par[2] - par[5])*lCosTht*lSinPsi)
308                                                                      *(x[0] - x[1] 
309                                                                        +(par[2] - par[5])*lCosTht*lSinPsi) 
310                                                                      +((par[0] - par[3])*(par[0] - par[3])
311                                                                        +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
312                                                                                                               +lSinPsi*lSinPsi*lSinTht*lSinTht))))
313                    /(((par[0] - par[3])*(par[0] - par[3]) 
314                       +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi 
315                                                              +lSinPsi*lSinPsi*lSinTht*lSinTht)));
316   
317   phiy = TMath::ACos(phiy);
318
319   return phiy;
320 }
321  
322 Double_t ynCenter(Double_t *x, Double_t *par){
323
324   Double_t lCosPsi = TMath::Cos(par[6]);
325   Double_t lSinPsi = TMath::Sin(par[6]);
326   Double_t lCosTht = TMath::Cos(par[7]);
327   Double_t lSinTht = TMath::Sin(par[7]);
328
329   Double_t yD = ((1./(+(par[0] - par[3])*(par[0] - par[3]) 
330                       +(par[1] - par[4])*(par[1] - par[4])))
331                  *(+par[3]*par[3]*x[0] 
332                    +par[0]*par[0]*x[1] 
333                    -par[0]*par[3]*(x[0] + x[1]) 
334                    +(par[1] - par[4])*(-x[0]*par[4] + par[1]*x[1]) 
335                    +(+par[3]*par[3]*par[2] 
336                      +par[0]*par[0]*par[5] 
337                      -par[0]*par[3]*(par[2] + par[5]) 
338                      +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lCosTht*lSinPsi 
339                    +(par[3]*par[1] - par[0]*par[4])
340                    *TMath::Sqrt(-(+x[0] - x[1] 
341                                   +(par[2] - par[5])*lCosTht*lSinPsi)
342                                 *(x[0] - x[1] 
343                                   +(par[2] - par[5])*lCosTht*lSinPsi) 
344                                 +((par[0] - par[3])*(par[0] - par[3]) 
345                                   +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
346                                                                          +lSinPsi*lSinPsi*lSinTht*lSinTht))));
347   
348   return yD;  
349 }
350  
351
352 Double_t phiynn(Double_t *x, Double_t *par){
353
354   Double_t lCosPsi = TMath::Cos(par[6]);
355   Double_t lSinPsi = TMath::Sin(par[6]);
356   Double_t lCosTht = TMath::Cos(par[7]);
357   Double_t lSinTht = TMath::Sin(par[7]);
358
359   Double_t phiy = ((lCosPsi*(+(par[1] - par[4])*(x[0] - x[1]) 
360                              +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi 
361                              +(par[0] - par[3])
362                              *TMath::Sqrt(-(x[0] - x[1] 
363                                             +(par[2] - par[5])*lCosTht*lSinPsi)
364                                           *(x[0] - x[1] 
365                                             +(par[2] - par[5])*lCosTht*lSinPsi) 
366                                           +(+(par[0] - par[3])*(par[0] - par[3]) 
367                                             +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
368                                                                                    +lSinPsi*lSinPsi*lSinTht*lSinTht))) 
369                     +lSinPsi*lSinTht*(+(par[0] - par[3])*(x[0] - x[1]) 
370                                       +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi 
371                                       +(-par[1] + par[4])
372                                       *TMath::Sqrt(-(x[0] - x[1] 
373                                                      +(par[2] - par[5])*lCosTht*lSinPsi)
374                                                    *(x[0] - x[1] 
375                                                      +(par[2] - par[5])*lCosTht*lSinPsi) 
376                                                    +(+(par[0] - par[3])*(par[0] - par[3]) 
377                                                      +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
378                                                                                             +lSinPsi*lSinPsi*lSinTht*lSinTht))))
379                    /((+(par[0] - par[3])*(par[0] - par[3]) 
380                       +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
381                                                              +lSinPsi*lSinPsi*lSinTht*lSinTht)));
382   
383   phiy = -TMath::ACos(phiy);
384   
385   return phiy;
386 }
387
388
389 Double_t phiynp(Double_t *x, Double_t *par){
390
391   Double_t lCosPsi = TMath::Cos(par[6]);
392   Double_t lSinPsi = TMath::Sin(par[6]);
393   Double_t lCosTht = TMath::Cos(par[7]);
394   Double_t lSinTht = TMath::Sin(par[7]);
395
396   Double_t phiy = ((lCosPsi*(+(par[1] - par[4])*(x[0] - x[1]) 
397                              +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi 
398                              +(par[0] - par[3])
399                              *TMath::Sqrt(-(x[0] - x[1] 
400                                             +(par[2] - par[5])*lCosTht*lSinPsi)
401                                           *(x[0] - x[1] 
402                                             +(par[2] - par[5])*lCosTht*lSinPsi) 
403                                           +((par[0] - par[3])*(par[0] - par[3]) 
404                                             +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
405                                                                                    +lSinPsi*lSinPsi*lSinTht*lSinTht))) 
406                     +lSinPsi*lSinTht*(+(par[0] - par[3])*(x[0] - x[1]) 
407                                       +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi 
408                                       +(-par[1] + par[4])
409                                       *TMath::Sqrt(-(x[0] - x[1] 
410                                                      +(par[2] - par[5])*lCosTht*lSinPsi)
411                                                    *(x[0] - x[1] 
412                                                      +(par[2] - par[5])*lCosTht*lSinPsi) 
413                                                    +((par[0] - par[3])*(par[0] - par[3]) 
414                                                      +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
415                                                                                             +lSinPsi*lSinPsi*lSinTht*lSinTht))))
416                    /((+(par[0] - par[3])*(par[0] - par[3])
417                       +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi 
418                                                              +lSinPsi*lSinPsi*lSinTht*lSinTht)));
419   
420   phiy = TMath::ACos(phiy);
421   
422   return phiy;
423 }
424
425 Double_t znCenter(Double_t *x, Double_t *par){
426   // par :  x1l, y1l, z1l, x2l, y2l, z2l, psi, tht
427
428   Double_t lCosPsi = TMath::Cos(par[6]);
429   Double_t lSinPsi = TMath::Sin(par[6]);
430   Double_t lCosTht = TMath::Cos(par[7]);
431   Double_t lSinTht = TMath::Sin(par[7]);
432
433   Double_t inSqrt = ((par[3]*par[1] - par[0]*par[4])*(par[3]*par[1] - par[0]*par[4])
434                      *((-(x[0] - x[1])*(x[0] - x[1]))
435                        +(((par[0] - par[3])*(par[0] - par[3])
436                           +(par[1] - par[4])*(par[1] - par[4])))*lSinPsi*lSinPsi
437                        +lCosPsi*((-(par[2] - par[5]))
438                                  *lCosTht*(-2*x[0]+2*x[1]
439                                            +(par[2] - par[5])*lCosPsi*lCosTht)
440                                  +((par[0] - par[3])*(par[0] - par[3])
441                                    +(par[1] - par[4])*(par[1] - par[4]))*lCosPsi*lSinTht*lSinTht)));
442
443   Double_t zD = ((1./((par[0] - par[3])*(par[0] - par[3]) 
444                       +(par[1] - par[4])*(par[1] - par[4])))
445                  *(-par[1]*par[4]*x[0] 
446                    +par[4]*par[4]*x[0] 
447                    +par[0]*par[0]*x[1] 
448                    +par[1]*par[1]*x[1] 
449                    -par[1]*par[4]*x[1] 
450                    -par[0]*par[3]*(x[0] + x[1]) 
451                    +par[3]*par[3]*x[0]
452                    +(+par[1]*par[4]*par[2] 
453                      -par[4]*par[4]*par[2] 
454                      -par[0]*par[0]*par[5] 
455                      -par[1]*par[1]*par[5] 
456                      +par[1]*par[4]*par[5] 
457                      +par[0]*par[3]*(par[2] + par[5])
458                      -par[3]*par[3]*par[2])*lCosPsi*lCosTht
459                    -TMath::Sqrt(inSqrt)));
460   
461   return zD;
462 }
463
464 Double_t zpCenter(Double_t *x, Double_t *par){
465   // par :  x1l, y1l, z1l, x2l, y2l, z2l, psi, tht
466
467   Double_t lCosPsi = TMath::Cos(par[6]);
468   Double_t lSinPsi = TMath::Sin(par[6]);
469   Double_t lCosTht = TMath::Cos(par[7]);
470   Double_t lSinTht = TMath::Sin(par[7]);
471
472   Double_t inSqrt = ((par[3]*par[1] - par[0]*par[4])*(par[3]*par[1] - par[0]*par[4])
473                      *((-(x[0] - x[1])*(x[0] - x[1]))
474                        +(((par[0] - par[3])*(par[0] - par[3])
475                           +(par[1] - par[4])*(par[1] - par[4])))*lSinPsi*lSinPsi
476                        +lCosPsi*((-(par[2] - par[5]))
477                                  *lCosTht*(-2*x[0]+2*x[1]
478                                            +(par[2] - par[5])*lCosPsi*lCosTht)
479                                  +((par[0] - par[3])*(par[0] - par[3])
480                                    +(par[1] - par[4])*(par[1] - par[4]))*lCosPsi*lSinTht*lSinTht)));
481   
482   Double_t zD = ((1./((par[0] - par[3])*(par[0] - par[3]) 
483                       +(par[1] - par[4])*(par[1] - par[4])))
484                  *(-par[1]*par[4]*x[0] 
485                    +par[4]*par[4]*x[0] 
486                    +par[0]*par[0]*x[1] 
487                    +par[1]*par[1]*x[1] 
488                    -par[1]*par[4]*x[1] 
489                    -par[0]*par[3]*(x[0] + x[1]) 
490                    +par[3]*par[3]*x[0]
491                    +(+par[1]*par[4]*par[2] 
492                      -par[4]*par[4]*par[2] 
493                      -par[0]*par[0]*par[5] 
494                      -par[1]*par[1]*par[5] 
495                      +par[1]*par[4]*par[5] 
496                      +par[0]*par[3]*(par[2] + par[5])
497                      -par[3]*par[3]*par[2])*lCosPsi*lCosTht
498                    +TMath::Sqrt(inSqrt)));
499
500   return zD;
501 }
502
503 //______________________________________________________________________
504 AliMUONGeometryTransformer *ReAlign(const AliMUONGeometryTransformer * transformer, 
505                                     int rMod, TGeoCombiTrans deltaDetElemTransf[], Bool_t verbose)
506 {
507   /////////////////////////////////////////////////////////////////////
508   //   Takes the internal geometry module transformers, copies them
509   // and gets the Detection Elements from them.
510   // Takes misalignment parameters and applies these
511   // to the local transform of the Detection Element
512   // Obtains the global transform by multiplying the module transformer
513   // transformation with the local transformation 
514   // Applies the global transform to a new detection element
515   // Adds the new detection element to a new module transformer
516   // Adds the new module transformer to a new geometry transformer
517   // Returns the new geometry transformer
518
519
520   Int_t iDetElemId = 0;
521   Int_t iDetElemNumber = 0;
522   Int_t iDetElemIndex = 0;
523   Int_t iCh = 0;
524
525   AliMUONGeometryTransformer *newGeometryTransformer =
526     new AliMUONGeometryTransformer();
527   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
528     // module transformers    
529     const AliMUONGeometryModuleTransformer *kModuleTransformer =
530       transformer->GetModuleTransformer(iMt, true);
531       
532     AliMUONGeometryModuleTransformer *newModuleTransformer =
533       new AliMUONGeometryModuleTransformer(iMt);
534     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
535     
536     TGeoCombiTrans moduleTransform =
537       TGeoCombiTrans(*kModuleTransformer->GetTransformation());
538     // New module transformation
539     TGeoCombiTrans *newModuleTransform;
540     if (iMt==rMod) {
541       newModuleTransform = new TGeoCombiTrans(moduleTransform);
542     } else {
543       newModuleTransform = new TGeoCombiTrans(moduleTransform);
544     }
545     newModuleTransformer->SetTransformation(*newModuleTransform);
546     
547     // For the selected chamber add misalign module
548     if (iMt==rMod) {
549       // Get delta transformation: 
550       // Tdelta = Tnew * Told.inverse
551       TGeoHMatrix deltaModuleTransform = 
552         AliMUONGeometryBuilder::Multiply(*newModuleTransform, 
553                                          kModuleTransformer->GetTransformation()->Inverse());    
554       // Create module mis alignment matrix
555       newGeometryTransformer
556         ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
557     }
558       
559     AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
560     
561     if (verbose)
562       printf("%i DEs in old GeometryStore  %i\n",detElements->GetSize(), iMt);
563     TGeoCombiTrans *deltaLocalTransform;
564     TIter next(detElements->CreateIterator());
565     AliMUONGeometryDetElement *detElement;
566     while ((detElement = static_cast<AliMUONGeometryDetElement*>(next()))){
567       /// make a new detection element
568       AliMUONGeometryDetElement *newDetElement =
569         new AliMUONGeometryDetElement(detElement->GetId(),
570                                       detElement->GetVolumePath());
571       TString lDetElemName(detElement->GetDEName());
572       lDetElemName.ReplaceAll("DE","");
573       iDetElemId = lDetElemName.Atoi();
574       iDetElemNumber = iDetElemId%100;
575       iCh = iDetElemId/100 -1;
576       if(iMt==rMod){
577         if (iCh<4) {
578           iDetElemIndex = iDetElemId;
579         } else {
580           if ((iDetElemNumber > (fgNDetElemCh[iCh]-2)/4) &&
581               (iDetElemNumber < fgNDetElemCh[iCh]-(fgNDetElemCh[iCh]-2)/4)) {
582             iDetElemIndex = (+fgNDetElemCh[iCh] 
583                              -(1+(fgNDetElemCh[iCh]-2)/4) 
584                              -iDetElemNumber);
585           } else {
586             iDetElemIndex = (+fgNDetElemCh[iCh] 
587                              -fgNDetElemCh[iCh]/2
588                              -((1+(fgNDetElemCh[iCh]-2)/4) 
589                                -TMath::Min(iDetElemNumber,
590                                            TMath::Abs(iDetElemNumber-fgNDetElemCh[iCh]))));
591           }
592         }
593         deltaLocalTransform = new TGeoCombiTrans(deltaDetElemTransf[iDetElemIndex]);       
594       } else {
595         deltaLocalTransform = new TGeoCombiTrans(*gGeoIdentity);
596       }
597
598       // local transformation of this detection element.
599       TGeoCombiTrans localTransform
600         = TGeoCombiTrans(*detElement->GetLocalTransformation());
601       //      TGeoHMatrix newLocalMatrix = localTransform * (*deltaLocalTransform);
602       TGeoCombiTrans newLocalTransform 
603         = TGeoCombiTrans(localTransform * (*deltaLocalTransform));
604       newDetElement->SetLocalTransformation(newLocalTransform);   
605       // global transformation
606       TGeoHMatrix newGlobalTransform =
607         AliMUONGeometryBuilder::Multiply(*newModuleTransform,
608                                          newLocalTransform);
609       newDetElement->SetGlobalTransformation(newGlobalTransform);
610       
611       // add this det element to module
612       newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
613                                                       newDetElement);
614       
615       // In the Alice Alignment Framework misalignment objects store
616       // global delta transformation
617       // Get detection "intermediate" global transformation
618       TGeoHMatrix newOldGlobalTransform = (*newModuleTransform) * localTransform;
619       // Get detection element global delta transformation: 
620       // Tdelta = Tnew * Told.inverse
621       TGeoHMatrix  deltaGlobalTransform
622         = AliMUONGeometryBuilder::Multiply(newGlobalTransform, 
623                                            newOldGlobalTransform.Inverse());
624       
625       // Create mis alignment matrix
626       newGeometryTransformer
627         ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
628     }
629     
630     if (verbose)
631       printf("Added module transformer %i to the transformer\n", iMt);
632     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
633   }
634   return newGeometryTransformer;
635 }