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 }