ITS/AliITSAlignMille2Module.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 /// \class AliITSAlignMille2Module
  19 /// Alignment class for the ALICE ITS detector
  20 ///
  21 /// This class is used by AliITSAlignMille to build custom supermodules
  22 /// made of ITS sensitive modules. These supermodules are then aligned
  23 ///
  24 /// Custom supermodules must have VolumeID > 14335
  25 ///
  26 /// \author M. Lunardon
  27 //-----------------------------------------------------------------------------
  28
  29 #include <TGeoManager.h>
  30 #include <TGeoMatrix.h>
  31
  32 #include "AliITSAlignMille2Module.h"
  33 #include "AliITSgeomTGeo.h"
  34 #include "AliGeomManager.h"
  35 #include "AliAlignObjParams.h"
  36 #include "AliLog.h"
  37 #include "AliITSAlignMille2.h"
  38
  39 /// \cond CLASSIMP
  40 ClassImp(AliITSAlignMille2Module)
  41 /// \endcond
  42
  43 #define CORHW_
  44
  45
  46 const Float_t AliITSAlignMille2Module::fgkDummyConstraint = 1e-2;//1.E3;
  47
  48 //-------------------------------------------------------------
  49 AliITSAlignMille2Module::AliITSAlignMille2Module() :
  50   TNamed(),
  51   fNSensVol(0),
  52   fIndex(-1),
  53   fDetType(-1),
  54   fVolumeID(0),
  55   fNParTot(0),
  56   fNParFree(0),
  57   fParOffs(0),
  58   fNProcPoints(0),
  59   fParVals(0),
  60   fParErrs(0),
  61   fParCstr(0),
  62   fSensVolIndex(0),
  63   fSensVolVolumeID(0),
  64   fMatrix(NULL),
  65   fSensVolMatrix(NULL),
  66   fSensVolModifMatrix(NULL),
  67   fParent(NULL),
  68   fChildren(0)
  69 {
  70   /// void constructor
  71   fMatrix = new TGeoHMatrix;
  72   fSensVolMatrix = new TGeoHMatrix;
  73   fSensVolModifMatrix = new TGeoHMatrix;
  74   fSensVolIndex.Set(1);
  75   fSensVolVolumeID.Set(1);
  76   fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0;
  77 }
  78
  79 //-------------------------------------------------------------
  80 AliITSAlignMille2Module::AliITSAlignMille2Module(Int_t index, UShort_t volid, const char* symname,
  81                                                  const TGeoHMatrix *m, Int_t nsv, const UShort_t *volidsv) :
  82   TNamed(),
  83   fNSensVol(0),
  84   fIndex(-1),
  85   fDetType(-1),
  86   fVolumeID(0),
  87   fNParTot(0),
  88   fNParFree(0),
  89   fParOffs(kMaxParGeom),
  90   fNProcPoints(0),
  91   fParVals(0),
  92   fParErrs(0),
  93   fParCstr(0),
  94   fSensVolIndex(0),
  95   fSensVolVolumeID(0),
  96   fMatrix(NULL),
  97   fSensVolMatrix(NULL),
  98   fSensVolModifMatrix(NULL),
  99   fParent(NULL),
 100   fChildren(0)
 101 {
 102   /// void constructor
 103   fMatrix = new TGeoHMatrix;
 104   fSensVolMatrix = new TGeoHMatrix;
 105   fSensVolModifMatrix = new TGeoHMatrix;
 106   fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0;
 107   for (int i=kMaxParGeom;i--;) fParOffs[i] = -1;
 108   if (Set(index,volid,symname,m,nsv,volidsv)) {
 109     AliInfo("Error in AliITSAlignMille2Module::Set() - initializing void supermodule...");
 110   }
 111   AssignDetType();
 112 }
 113
 114 //-------------------------------------------------------------
 115 AliITSAlignMille2Module::AliITSAlignMille2Module(UShort_t volid) :
 116   TNamed(),
 117   fNSensVol(0),
 118   fIndex(-1),
 119   fDetType(-1),
 120   fVolumeID(0),
 121   fNParTot(0),
 122   fNParFree(0),
 123   fParOffs(kMaxParGeom),
 124   fNProcPoints(0),
 125   fParVals(0),
 126   fParErrs(0),
 127   fParCstr(0),
 128   fSensVolIndex(0),
 129   fSensVolVolumeID(0),
 130   fMatrix(NULL),
 131   fSensVolMatrix(NULL),
 132   fSensVolModifMatrix(NULL),
 133   fParent(NULL),
 134   fChildren(0)
 135 {
 136   /// simple constructor building a supermodule from a single sensitive volume
 137   fMatrix = new TGeoHMatrix;
 138   fSensVolMatrix = new TGeoHMatrix;
 139   fSensVolModifMatrix = new TGeoHMatrix;
 140   // temporary align object, just use the rotation...
 141   fSensVolIndex.Set(1);
 142   fSensVolVolumeID.Set(1);
 143   fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0;
 144   for (int i=kMaxParGeom;i--;) fParOffs[i] = -1;
 145   //
 146   fIndex = GetIndexFromVolumeID(volid);
 147   if (fIndex>=0 && gGeoManager) { // good sensitive module and geometry loaded
 148     SetName(AliGeomManager::SymName(volid));
 149     fVolumeID = volid;
 150     AddSensitiveVolume(volid);
 151     SetSensorsProvided(kTRUE);
 152     if (SensVolMatrix(volid, fMatrix))
 153        AliInfo("Matrix not defined");
 154   }
 155   else {
 156     AliInfo("Wrong VolumeID or Geometry not loaded - initializing void supermodule...");
 157   }
 158   AssignDetType();
 159 }
 160
 161
 162 //_____________________________________________________________________________
 163 AliITSAlignMille2Module::AliITSAlignMille2Module(const AliITSAlignMille2Module &m) :
 164   TNamed(m),
 165   fNSensVol(m.fNSensVol),
 166   fIndex(m.fIndex),
 167   fDetType(m.fDetType),
 168   fVolumeID(m.fVolumeID),
 169   fNParTot(m.fNParTot),
 170   fNParFree(m.fNParFree),
 171   fParOffs(m.fNParTot),
 172   fNProcPoints(0),
 173   fParVals(0),
 174   fParErrs(0),
 175   fParCstr(0),
 176   fSensVolIndex(m.fSensVolIndex),
 177   fSensVolVolumeID(m.fSensVolVolumeID),
 178   fMatrix(new TGeoHMatrix(*m.GetMatrix())),
 179   fSensVolMatrix(new TGeoHMatrix),
 180   fSensVolModifMatrix(new TGeoHMatrix),
 181   fParent(m.fParent),
 182   fChildren(0)
 183 {
 184   // Copy constructor
 185   fSensVolIndex = m.fSensVolIndex;
 186   fSensVolVolumeID = m.fSensVolVolumeID;
 187   for (int i=m.fNParTot;i--;) fParOffs[i] = m.fParOffs[i];
 188   for (int i=3;i--;) fSigmaFactor[i] = m.fSigmaFactor[i];
 189   if (fNParTot) {
 190     fParVals = new Float_t[fNParTot];
 191     fParErrs = new Float_t[fNParTot];
 192     fParCstr = new Float_t[fNParTot];
 193     for (int i=fNParTot;i--;) {
 194       fParVals[i] = m.fParVals[i];
 195       fParErrs[i] = m.fParErrs[i];
 196       fParCstr[i] = m.fParCstr[i];
 197     }
 198   }
 199 }
 200
 201 //_____________________________________________________________________________
 202 AliITSAlignMille2Module& AliITSAlignMille2Module::operator=(const AliITSAlignMille2Module &m)
 203 {
 204   // operator =
 205   //
 206   if(this==&m) return *this;
 207   this->TNamed::operator=(m);
 208   //
 209   fNSensVol=m.fNSensVol;
 210   fIndex=m.fIndex;
 211   fDetType = m.fDetType;
 212   fVolumeID=m.fVolumeID;
 213   fNParTot  = m.fNParTot;
 214   fNParFree = m.fNParFree;
 215   fNProcPoints = m.fNProcPoints;
 216   delete[] fParVals; fParVals = 0;
 217   delete[] fParErrs; fParErrs = 0;
 218   delete[] fParCstr; fParCstr = 0;
 219   //
 220   if (fNParTot) {
 221     fParVals = new Float_t[fNParTot];
 222     fParErrs = new Float_t[fNParTot];
 223     fParCstr = new Float_t[fNParTot];
 224     for (int i=m.GetNParTot();i--;) {
 225       fParVals[i] = m.fParVals[i];
 226       fParErrs[i] = m.fParErrs[i];
 227       fParCstr[i] = m.fParCstr[i];
 228     }
 229   }
 230   //
 231   fParOffs.Set(fNParTot);
 232   for (int i=0;i<fNParTot;i++) fParOffs[i] = m.fParOffs[i];
 233   for (int i=0;i<3;i++) fSigmaFactor[i] = m.fSigmaFactor[i];
 234   if (fMatrix) delete fMatrix;
 235   fMatrix=new TGeoHMatrix(*m.GetMatrix());
 236   if(fSensVolMatrix) delete fSensVolMatrix;
 237   fSensVolMatrix = new TGeoHMatrix(*m.fSensVolMatrix);
 238   if(fSensVolModifMatrix) delete fSensVolModifMatrix;
 239   fSensVolModifMatrix = new TGeoHMatrix(*m.fSensVolModifMatrix);
 240   fSensVolIndex = m.fSensVolIndex;
 241   fSensVolVolumeID = m.fSensVolVolumeID;
 242   fParent = m.fParent;
 243   fChildren.Clear();
 244   for (int i=0;i<m.GetNChildren();i++) fChildren.Add(m.GetChild(i));
 245   return *this;
 246 }
 247
 248
 249 //-------------------------------------------------------------
 250 AliITSAlignMille2Module::~AliITSAlignMille2Module() {
 251   /// Destructor
 252   delete fMatrix;
 253   delete fSensVolMatrix;
 254   delete fSensVolModifMatrix;
 255   delete[] fParVals;
 256   delete[] fParErrs;
 257   delete[] fParCstr;
 258   fChildren.Clear();
 259 }
 260
 261 //-------------------------------------------------------------
 262 Int_t AliITSAlignMille2Module::Set(Int_t index, UShort_t volid, const char* symname,
 263                                    const TGeoHMatrix *m, Int_t nsv, const UShort_t *volidsv)
 264 {
 265   // initialize a custom supermodule
 266   // index, volid, symname and matrix must be given
 267   // if (volidsv) add nsv sensitive volumes to the supermodules
 268   // return 0 if success
 269
 270   if (index<2198) {
 271     AliInfo("Index must be >= 2198");
 272     return -1;
 273   }
 274   if (volid<14336) {
 275     AliInfo("VolumeID must be >= 14336");
 276     return -2;
 277   }
 278
 279   if (!symname) return -3;
 280   for (Int_t i=0; i<2198; i++) {
 281     if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) {
 282       AliInfo("Symname already used by a Sensitive Volume");
 283       return -3;
 284     }
 285   }
 286
 287   if (!m) return -4;
 288
 289   // can initialize needed stuffs
 290   fIndex = index;
 291   fVolumeID = volid;
 292   SetName(symname);
 293   //
 294   (*fMatrix) = (*m);
 295   //
 296   fSensVolIndex.Set(nsv);
 297   fSensVolVolumeID.Set(nsv);
 298   // add sensitive volumes
 299   for (Int_t i=0; i<nsv; i++) AddSensitiveVolume(volidsv[i]);
 300
 301   return 0;
 302 }
 303
 304 //-------------------------------------------------------------
 305 void AliITSAlignMille2Module::SetFreeDOF(Int_t dof,Double_t cstr)
 306 {
 307   if (AliITSAlignMille2::IsZero(cstr)) fParCstr[dof] = 0;  // fixed parameter
 308   else if (cstr>0)                     fParCstr[dof] = fgkDummyConstraint+1.; // the parameter is free and unconstrained
 309   else                                 fParCstr[dof] = -cstr;                 // the parameter is free but constrained
 310 }
 311
 312 //-------------------------------------------------------------
 313 Bool_t AliITSAlignMille2Module::IsSensor(UShort_t voluid)
 314 {
 315   // Does this volid correspond to sensor ?
 316   AliGeomManager::ELayerID layId = AliGeomManager::VolUIDToLayerSafe(voluid);
 317   if (layId>0 && layId<7) {
 318     Int_t mId = Int_t(voluid & 0x7ff);
 319     if( mId>=0 && mId<AliGeomManager::LayerSize(layId)) return kTRUE;
 320   }
 321   return kFALSE;
 322 }
 323
 324 //-------------------------------------------------------------
 325 Int_t AliITSAlignMille2Module::GetIndexFromVolumeID(UShort_t voluid) {
 326   /// index from volume ID
 327   AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
 328   if (lay<1|| lay>6) return -1;
 329   Int_t idx=Int_t(voluid)-2048*lay;
 330   if (idx>=AliGeomManager::LayerSize(lay)) return -1;
 331   for (Int_t ilay=1; ilay<lay; ilay++)
 332     idx += AliGeomManager::LayerSize(ilay);
 333   return idx;
 334 }
 335
 336 //-------------------------------------------------------------
 337 void AliITSAlignMille2Module::AddSensitiveVolume(UShort_t voluid)
 338 {
 339   /// add a sensitive volume to this supermodule
 340   if (GetIndexFromVolumeID(voluid)<0) return; // bad volid
 341   //
 342   // in principle, the correct size of fSensVol... arrays was set outside but check anyway
 343   if (fSensVolVolumeID.GetSize()<fNSensVol+1) {
 344     fSensVolVolumeID.Set(fNSensVol+1);
 345     fSensVolIndex.Set(fNSensVol+1);
 346   }
 347   //
 348   fSensVolVolumeID[fNSensVol] = Short_t(voluid);
 349   fSensVolIndex[fNSensVol] = GetIndexFromVolumeID(voluid);
 350   fNSensVol++;
 351 }
 352
 353 //-------------------------------------------------------------
 354 void AliITSAlignMille2Module::DelSensitiveVolume(Int_t at)
 355 {
 356   // Suppress sensor at position "at"
 357   // in fact we are swapping with the last valid one
 358   int lastValid = --fNSensVol;
 359   int tmpv = fSensVolIndex[at];
 360   fSensVolIndex[at] = fSensVolIndex[lastValid];
 361   tmpv = fSensVolVolumeID[at];
 362   fSensVolVolumeID[at] = fSensVolVolumeID[lastValid];
 363   fSensVolVolumeID[lastValid] = tmpv;
 364   //
 365 }
 366
 367 //-------------------------------------------------------------
 368 Bool_t AliITSAlignMille2Module::IsIn(UShort_t voluid) const
 369 {
 370   /// check if voluid is defined
 371   if (!voluid) return kFALSE; // only positive voluid are accepted
 372   for (Int_t i=0; i<fNSensVol; i++) if (UShort_t(fSensVolVolumeID[i])==voluid) return kTRUE;
 373   return kFALSE;
 374 }
 375
 376 //-------------------------------------------------------------
 377 Bool_t AliITSAlignMille2Module::BelongsTo(AliITSAlignMille2Module* parent) const
 378 {
 379   /// check if parent contains the sensors of this volume
 380   if (fNSensVol<1 || fNSensVol>=parent->GetNSensitiveVolumes()) return kFALSE;
 381   return parent->IsIn( fSensVolVolumeID[0] );
 382 }
 383
 384 //-------------------------------------------------------------
 385 TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeModifiedMatrix(UShort_t voluid, const Double_t *delta,Bool_t local)
 386 {
 387   // modify the original TGeoHMatrix of the sensitive module 'voluid' according
 388   // with a delta transform. applied to the supermodule matrix
 389   // return NULL if error
 390
 391   if (!IsIn(voluid)) return NULL;
 392   if (!gGeoManager)  return NULL;
 393
 394   // prepare the TGeoHMatrix
 395   Double_t tr[3],ang[3];
 396   tr[0]=delta[0]; // in centimeter
 397   tr[1]=delta[1];
 398   tr[2]=delta[2];
 399   ang[0]=delta[3]; // psi   (X)  in deg
 400   ang[1]=delta[4]; // theta (Y)
 401   ang[2]=delta[5]; // phi   (Z)
 402   //
 403   static AliAlignObjParams tempAlignObj;
 404   tempAlignObj.SetRotation(ang[0],ang[1],ang[2]);
 405   tempAlignObj.SetTranslation(tr[0],tr[1],tr[2]);
 406   AliDebug(3,Form("Delta angles: psi=%f  theta=%f   phi=%f",ang[0],ang[1],ang[2]));
 407   TGeoHMatrix hm;
 408   tempAlignObj.GetMatrix(hm);
 409   //printf("\n0: delta matrix\n");hm.Print();
 410
 411   // 1) start setting fSensVolModif = fSensVol
 412   if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL;
 413   //
 414   if (local) {
 415     // 2) set fSensVolModif = SensVolRel
 416     fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() );
 417     // 3) multiply left by delta
 418     fSensVolModifMatrix->MultiplyLeft( &hm );
 419     // 4) multiply left by fMatrix
 420     fSensVolModifMatrix->MultiplyLeft( fMatrix );
 421   }
 422   else fSensVolModifMatrix->MultiplyLeft( &hm );
 423   //
 424   return fSensVolModifMatrix;
 425 }
 426
 427 //-------------------------------------------------------------
 428 AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeMisalignment(UShort_t voluid, const Double_t *deltalocal)
 429 {
 430   // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal'
 431   // of the mother volume. The misalignment is returned as AliAlignObjParams object
 432
 433   if (!IsIn(voluid)) return NULL;
 434   if (!gGeoManager) return NULL;
 435
 436   // prepare the TGeoHMatrix
 437   Double_t tr[3],ang[3];
 438   tr[0]=deltalocal[0]; // in centimeter
 439   tr[1]=deltalocal[1];
 440   tr[2]=deltalocal[2];
 441   ang[0]=deltalocal[3]; // psi   (X)  in deg
 442   ang[1]=deltalocal[4]; // theta (Y)
 443   ang[2]=deltalocal[5]; // phi   (Z)
 444   //
 445   static AliAlignObjParams tempAlignObj;
 446   tempAlignObj.SetRotation(ang[0],ang[1],ang[2]);
 447   tempAlignObj.SetTranslation(tr[0],tr[1],tr[2]);
 448   AliDebug(3,Form("Delta angles: psi=%f  theta=%f   phi=%f",ang[0],ang[1],ang[2]));
 449   //
 450   return GetSensitiveVolumeMisalignment(voluid,&tempAlignObj);
 451 }
 452
 453 //-------------------------------------------------------------
 454 AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeMisalignment(UShort_t voluid, const AliAlignObjParams *a)
 455 {
 456   // return the misalignment of the sens. vol. 'voluid' corresponding with
 457   // a misalignment 'a' in the mother volume
 458   // return NULL if error
 459
 460   // Gsv = Gg * Gg-1 * Gsv   -> Lsv,g = Gg-1 * Gsv
 461   // G'sv = Gg * Dg * Lsv,g === Gsv * Dsv
 462   // Gg * Dg * Gg-1 * Gsv = Gsv * Gsv-1 * Gg * Dg * Gg-1 * Gsv
 463   //
 464   // => Dsv = (Gsv-1 * Gg * Dg * Gg-1 * Gsv)
 465   //
 466
 467   if (!IsIn(voluid)) return NULL;
 468   if (!gGeoManager) return NULL;
 469
 470   //a->Print("");
 471
 472   // prepare the Delta matrix Dg
 473   TGeoHMatrix dg;
 474   a->GetMatrix(dg);
 475   //dg.Print();
 476
 477   // 1) start setting fSensVolModif = Gsv
 478   if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL;
 479   //printf("\n1: modif=orig del sensvol\n");fSensVolModifMatrix->Print();
 480
 481   // 2) set fSensVolModif = Gg-1 * Gsv
 482   fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() );
 483   //printf("\n2: modif=relative del sensvol\n");fSensVolModifMatrix->Print();
 484
 485   // 3) set fSensVolModif = Dg * Gg-1 * Gsv
 486   fSensVolModifMatrix->MultiplyLeft( &dg );
 487   //printf("\n3: modif= delta*relative\n");fSensVolModifMatrix->Print();
 488
 489   // 4) set fSensVolModif = Gg * Dg * Gg-1 * Gsv
 490   fSensVolModifMatrix->MultiplyLeft( fMatrix );
 491   //printf("\n4: modif=quasi,manca il Gsv-1...\n");fSensVolModifMatrix->Print();
 492
 493   // 5) set fSensVolModif = Gsv-1 * Gg * Dg * Gg-1 * Gsv
 494   if (SensVolMatrix(voluid, &dg)) return NULL;
 495   fSensVolModifMatrix->MultiplyLeft( &dg.Inverse() );
 496   //printf("\n5: modif=finale\n");fSensVolModifMatrix->Print();
 497   //
 498   // >> RS
 499  // 6) mo' fSensVolModif dovrebbe essere la Dsv(loc) t.c. G'sv = Gsv*Dsv(loc)
 500   // per trasformarla in Dsv(loc rispetto al Gsv0, non modificato) dovrebbe essere:
 501   // Dsv(loc) -> Dpre * Dsv(loc) * Dpre-1
 502   //TGeoHMatrix dpre; // dpre = Gsv0-1*Gsv
 503   //if (SensVolOrigGlobalMatrix(voluid, &dg)) return NULL;
 504   //if (SensVolMatrix(voluid, &dpre)) return NULL;
 505   //dpre.MultiplyLeft( &dg.Inverse() );
 506   //fSensVolModifMatrix->Multiply( &dpre.Inverse() );
 507   //fSensVolModifMatrix->MultiplyLeft( &dpre );
 508   // direi che NON FUNZIONA!!!!
 509
 510   // << RS
 511
 512   // reset align object (may not be needed...)
 513   static AliAlignObjParams tempAlignObj;
 514   tempAlignObj.SetVolUID(0);
 515   tempAlignObj.SetSymName("");
 516   tempAlignObj.SetTranslation(0,0,0);
 517   tempAlignObj.SetRotation(0,0,0);
 518   //
 519   // >> RS
 520 #ifdef CORHW_
 521   // correction for SPD y-shift
 522   if (voluid>=2048 && voluid<4256) {
 523     TGeoHMatrix deltay;
 524     double dy[3]={0.,0.0081,0.};
 525     deltay.SetTranslation(dy);
 526     fSensVolModifMatrix->MultiplyLeft( &deltay );
 527     fSensVolModifMatrix->Multiply( &deltay.Inverse() );
 528   }
 529 #endif
 530   // << RS
 531   if (!tempAlignObj.SetMatrix(*fSensVolModifMatrix)) return NULL;
 532   tempAlignObj.SetVolUID(voluid);
 533   tempAlignObj.SetSymName(AliGeomManager::SymName(voluid));
 534   //
 535   return &tempAlignObj;
 536 }
 537
 538 // >> RS
 539 //-------------------------------------------------------------
 540 AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeTotalMisalignment(UShort_t voluid, const Double_t *deltalocal)
 541 {
 542   // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal'
 543   // of the mother volume. The misalignment is returned as AliAlignObjParams object including
 544   // the (evenctual) prealignment => no merging needed
 545
 546   if (!IsIn(voluid)) return NULL;
 547   if (!gGeoManager) return NULL;
 548
 549   // prepare the TGeoHMatrix
 550   Double_t tr[3],ang[3];
 551   tr[0]=deltalocal[0]; // in centimeter
 552   tr[1]=deltalocal[1];
 553   tr[2]=deltalocal[2];
 554   ang[0]=deltalocal[3]; // psi   (X)  in deg
 555   ang[1]=deltalocal[4]; // theta (Y)
 556   ang[2]=deltalocal[5]; // phi   (Z)
 557
 558   // reset align object (may not be needed...)
 559   static AliAlignObjParams tempAlignObj;
 560   tempAlignObj.SetVolUID(0);
 561   tempAlignObj.SetSymName("");
 562   tempAlignObj.SetRotation(ang[0],ang[1],ang[2]);
 563   tempAlignObj.SetTranslation(tr[0],tr[1],tr[2]);
 564   AliDebug(3,Form("Delta angles: psi=%f  theta=%f   phi=%f",ang[0],ang[1],ang[2]));
 565
 566   // Gsv = Gg * Gg-1 * Gsv   -> Lsv,g = Gg-1 * Gsv
 567   // G'sv = Gg * Dg * Lsv,g === DGsv * Gsv
 568   //
 569   // => Dsv = (G0sv-1 * Gg * Dg * Gg-1 * GMsv)  //
 570   //
 571
 572   // prepare the Delta matrix Dg
 573   TGeoHMatrix dg;
 574   tempAlignObj.GetMatrix(dg);
 575   //dg.Print();
 576
 577   // 1) start setting fSensVolModif = Gsv
 578   if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL;
 579   //printf("\n1: modif=orig del sensvol\n");fSensVolModifMatrix->Print();
 580
 581   // 2) set fSensVolModif = Gg-1 * Gsv
 582   fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() );
 583   //printf("\n2: modif=relative del sensvol\n");fSensVolModifMatrix->Print();
 584
 585   // 3) set fSensVolModif = Dg * Gg-1 * Gsv
 586   fSensVolModifMatrix->MultiplyLeft( &dg );
 587   //printf("\n3: modif= delta*relative\n");fSensVolModifMatrix->Print();
 588
 589   // 4) set fSensVolModif = Gg * Dg * Gg-1 * Gsv
 590   fSensVolModifMatrix->MultiplyLeft( fMatrix );
 591   //printf("\n4: modif=quasi,manca il Gsv-1...\n");fSensVolModifMatrix->Print();
 592
 593   // 5) set fSensVolModif = G0sv-1 * Gg * Dg * Gg-1 * Gsv
 594   // qui usa l'orig anziche' la prealigned...
 595   if (SensVolOrigGlobalMatrix(voluid, &dg)) return NULL;
 596   fSensVolModifMatrix->MultiplyLeft( &dg.Inverse() );
 597   //printf("\n5: modif=finale\n");fSensVolModifMatrix->Print();
 598
 599   // reset align object (may not be needed...)
 600   tempAlignObj.SetVolUID(0);
 601   tempAlignObj.SetSymName("");
 602   tempAlignObj.SetTranslation(0,0,0);
 603   tempAlignObj.SetRotation(0,0,0);
 604
 605 #ifdef CORHW_
 606   // correction for SPD y-shift
 607   if (voluid>=2048 && voluid<4256) {
 608     TGeoHMatrix deltay;
 609     double dy[3]={0.,0.0081,0.};
 610     deltay.SetTranslation(dy);
 611     fSensVolModifMatrix->MultiplyLeft( &deltay );
 612     fSensVolModifMatrix->Multiply( &deltay.Inverse() );
 613   }
 614 #endif
 615   if (!tempAlignObj.SetMatrix(*fSensVolModifMatrix)) return NULL;
 616   tempAlignObj.SetVolUID(voluid);
 617   tempAlignObj.SetSymName(AliGeomManager::SymName(voluid));
 618
 619
 620   //tempAlignObj.Print("");
 621
 622   return &tempAlignObj;
 623 }
 624 //-------------------------------------------------------------
 625
 626 //-------------------------------------------------------------
 627 AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeGlobalMisalignment(UShort_t voluid, const Double_t *deltalocal)
 628 {
 629   // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal'
 630   // of the mother volume. The misalignment is returned as AliAlignObjParams object
 631
 632   if (!IsIn(voluid)) return NULL;
 633   if (!gGeoManager) return NULL;
 634
 635   // prepare the TGeoHMatrix
 636   Double_t tr[3],ang[3];
 637   tr[0]=deltalocal[0]; // in centimeter
 638   tr[1]=deltalocal[1];
 639   tr[2]=deltalocal[2];
 640   ang[0]=deltalocal[3]; // psi   (X)  in deg
 641   ang[1]=deltalocal[4]; // theta (Y)
 642   ang[2]=deltalocal[5]; // phi   (Z)
 643
 644   // reset align object (may not be needed...)
 645   static AliAlignObjParams tempAlignObj;
 646   tempAlignObj.SetTranslation(0,0,0);
 647   tempAlignObj.SetRotation(0,0,0);
 648
 649   tempAlignObj.SetRotation(ang[0],ang[1],ang[2]);
 650   tempAlignObj.SetTranslation(tr[0],tr[1],tr[2]);
 651   AliDebug(3,Form("Delta angles: psi=%f  theta=%f   phi=%f",ang[0],ang[1],ang[2]));
 652
 653   // Gsv = Gg * Gg-1 * Gsv   -> Lsv,g = Gg-1 * Gsv
 654   // G'sv = Gg * Dg * Lsv,g === DGsv * Gsv
 655   //
 656   // => DGsv = (Gg * Dg * Gg-1)
 657   //
 658
 659   // prepare the Delta matrix Dg
 660   TGeoHMatrix dg;
 661   tempAlignObj.GetMatrix(dg);
 662   //dg.Print();
 663
 664   dg.MultiplyLeft( fMatrix );
 665   dg.Multiply( &fMatrix->Inverse() );
 666
 667   // reset align object (may not be needed...)
 668   tempAlignObj.SetTranslation(0,0,0);
 669   tempAlignObj.SetRotation(0,0,0);
 670
 671   tempAlignObj.SetVolUID(voluid);
 672   tempAlignObj.SetSymName(AliGeomManager::SymName(voluid));
 673
 674   if (!tempAlignObj.SetMatrix(dg)) return NULL;
 675
 676   //tempAlignObj.Print("");
 677
 678   return &tempAlignObj;
 679 }
 680 // << RS
 681
 682 //-------------------------------------------------------------
 683 TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeMatrix(UShort_t voluid)
 684 {
 685   // return TGeoHMatrix of the sens.vol. 'voluid' in the current geometry
 686   if (SensVolMatrix(voluid,fSensVolMatrix)) return NULL;
 687   return fSensVolMatrix;
 688 }
 689
 690 //-------------------------------------------------------------
 691 TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeOrigGlobalMatrix(UShort_t voluid)
 692 {
 693   // return original ideal position (from AliGeomManager::GetOrigGlobalMatrix())
 694   if (SensVolOrigGlobalMatrix(voluid,fSensVolMatrix)) return NULL;
 695   return fSensVolMatrix;
 696 }
 697 //-------------------------------------------------------------
 698 Int_t AliITSAlignMille2Module::SensVolMatrix(UShort_t volid, TGeoHMatrix *m)
 699 {
 700   // set matrix for sensitive modules (SPD corrected)
 701   // return 0 if success
 702   Double_t rot[9];
 703   Int_t idx=GetIndexFromVolumeID(volid);
 704   if (idx<0) return -1;
 705   if (!AliITSgeomTGeo::GetRotation(idx,rot)) return -2;
 706   m->SetRotation(rot);
 707   Double_t oLoc[3]={0,0,0};
 708   Double_t oGlo[3]={0,0,0};
 709   if (!AliITSgeomTGeo::LocalToGlobal(idx,oLoc,oGlo)) return -3;
 710   m->SetTranslation(oGlo);
 711   return 0;
 712 }
 713
 714 //-------------------------------------------------------------
 715 Int_t AliITSAlignMille2Module::SensVolOrigGlobalMatrix(UShort_t volid, TGeoHMatrix *m)
 716 {
 717   // set original global matrix for sensitive modules (SPD corrected)
 718   // return 0 if success
 719   Int_t idx=GetIndexFromVolumeID(volid);
 720   if (idx<0) return -1;
 721   TGeoHMatrix mo;
 722   if (!AliGeomManager::GetOrigGlobalMatrix(AliGeomManager::SymName(volid),mo)) return -1;
 723   (*m)=mo;
 724   //
 725 #ifdef CORHW_
 726   // SPD y-shift by 81 mu
 727   if (volid<5000) {
 728     Double_t oLoc[3]={0.0,0.0081,0.0};
 729     Double_t oGlo[3]={0,0,0};
 730     m->LocalToMaster(oLoc,oGlo);
 731     m->SetTranslation(oGlo);
 732   }
 733 #endif
 734   return 0;
 735 }
 736
 737 //-------------------------------------------------------------
 738 UShort_t AliITSAlignMille2Module::GetVolumeIDFromSymname(const Char_t *symname) {
 739   /// volume ID from symname
 740   if (!symname) return 0;
 741
 742   for (UShort_t voluid=2000; voluid<13300; voluid++) {
 743     Int_t modId;
 744     AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
 745     if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
 746       if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
 747     }
 748   }
 749
 750   return 0;
 751 }
 752
 753 //-------------------------------------------------------------
 754 UShort_t AliITSAlignMille2Module::GetVolumeIDFromIndex(Int_t index) {
 755   /// volume ID from index
 756   if (index<0 || index>2197) return 0;
 757   return GetVolumeIDFromSymname(AliITSgeomTGeo::GetSymName(index));
 758 }
 759
 760 //-------------------------------------------------------------
 761 void AliITSAlignMille2Module::Print(Option_t*) const
 762 {
 763   // print data
 764   //
 765   const char* typeName[] = {"SPD","SDD","SSD"};
 766   printf("*** ITS SuperModule for AliITSAlignMille ***\n");
 767   printf("symname  : %s (type: %s)\n",GetName(),fDetType<0 ? "N/A":typeName[fDetType]);
 768   printf("parent   : %s | %d children\n",fParent ? fParent->GetName() : "N/A",GetNChildren());
 769   printf("volumeID : %4d  | index : %4d | Geom.Params are %s\n",fVolumeID,fIndex,
 770          GeomParamsGlobal() ? "Global":"Local");
 771   printf("Factors  : X=%.2f Y=%.2f Z=%.2f\n"
 772          "DOF: %cTx:%5d| %cTy:%5d| %cTz:%5d| %cPsi:%5d| %cTheta:%5d| %cPhi:%5d|",
 773          fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2],
 774          IsFreeDOF(kDOFTX) ? '+':'-',GetParOffset(kDOFTX),IsFreeDOF(kDOFTY) ? '+':'-',GetParOffset(kDOFTY),
 775          IsFreeDOF(kDOFTZ) ? '+':'-',GetParOffset(kDOFTZ),IsFreeDOF(kDOFPS) ? '+':'-',GetParOffset(kDOFPS),
 776          IsFreeDOF(kDOFTH) ? '+':'-',GetParOffset(kDOFTH),IsFreeDOF(kDOFPH) ? '+':'-',GetParOffset(kDOFPH));
 777   if (IsSDD()) {
 778     printf("%cT0:%5d| %cDVl:%5d| %cDVr:%5d|",IsFreeDOF(kDOFT0)?'+':'-',GetParOffset(kDOFT0),
 779            IsFreeDOF(kDOFDVL)?'+':'-',GetParOffset(kDOFDVL),IsFreeDOF(kDOFDVR)?'+':'-',GetParOffset(kDOFDVR));
 780     if (IsVDriftLRSame()) printf("(dVL=dVR)");
 781   }
 782   printf("\n");
 783   fMatrix->Print();
 784   printf("%4d Sensitive volumes | %6d Processed Points\n",fNSensVol,fNProcPoints);
 785   for (Int_t i=0; i<fNSensVol; i++) printf("   voluid[%d] = %d\n",i,UShort_t(fSensVolVolumeID[i]));
 786 }
 787
 788 //-------------------------------------------------------------
 789 Bool_t AliITSAlignMille2Module::IsAlignable() const
 790 {
 791   // it it alignable?
 792   TGeoManager* geoManager = AliGeomManager::GetGeometry();
 793   if (!geoManager) {
 794     AliInfo("Couldn't initialize geometry");
 795     return kFALSE;
 796   }
 797   return geoManager->GetAlignableEntry(GetName())!=0;
 798 }
 799
 800 //-------------------------------------------------------------
 801 void AliITSAlignMille2Module::GetLocalMatrix(TGeoHMatrix &mat) const
 802 {
 803   // return the local matrix for transformation to its parent
 804   mat = *fMatrix;
 805   if (fParent) mat.MultiplyLeft( &fParent->GetMatrix()->Inverse() );
 806 }
 807
 808 //-------------------------------------------------------------
 809 void AliITSAlignMille2Module::AssignDetType()
 810 {
 811   // assign the detector type
 812   TString tp = GetName();
 813   if      (tp.Contains("SPD",TString::kIgnoreCase)) fDetType = kSPD;
 814   else if (tp.Contains("SDD",TString::kIgnoreCase)) fDetType = kSDD;
 815   else if (tp.Contains("SSD",TString::kIgnoreCase)) fDetType = kSSD;
 816   else fDetType = -1;
 817   fNParTot = IsSDD() ? kMaxParTot:kMaxParGeom;
 818   fNParFree = 0;
 819   fParVals = new Float_t[fNParTot];
 820   fParErrs = new Float_t[fNParTot];
 821   fParCstr = new Float_t[fNParTot];
 822   if (fParOffs.GetSize()<fNParTot) fParOffs.Set(fNParTot);
 823   for (int i=fNParTot;i--;) {
 824     fParVals[i] = fParErrs[i] = 0.;
 825     fParCstr[i] = 0.;
 826     fParOffs[i] = -1;
 827   }
 828 }
 829
 830 //-------------------------------------------------------------
 831 void AliITSAlignMille2Module::EvaluateDOF()
 832 {
 833   // count d.o.f.
 834   fNParFree = 0;
 835   for (int i=fNParTot;i--;) if (IsFreeDOF(i)) fNParFree++;
 836 }
 837
 838 //-------------------------------------------------------------
 839 void AliITSAlignMille2Module::GetSensVolGlobalParams(UShort_t volid,Double_t *t, Double_t *r)
 840 {
 841   // return global parameters of the sensor volid
 842   for (int i=3;i--;) t[i] = r[i] = 0.;
 843   if (SensVolMatrix(volid,fSensVolMatrix)) return;
 844   AliAlignObjParams tempAlignObj;
 845   tempAlignObj.SetMatrix(*fSensVolMatrix);
 846   tempAlignObj.GetPars(t,r);
 847 }
 848
 849 //-------------------------------------------------------------
 850 void AliITSAlignMille2Module::GetSensVolLocalParams(UShort_t volid,Double_t *t, Double_t *r)
 851 {
 852   // return parameters of the sensor volid in the current module
 853   for (int i=3;i--;) t[i] = r[i] = 0.;
 854   if (SensVolMatrix(volid,fSensVolMatrix)) return;
 855   fSensVolMatrix->MultiplyLeft( &fMatrix->Inverse() );
 856   AliAlignObjParams tempAlignObj;
 857   tempAlignObj.SetMatrix(*fSensVolMatrix);
 858   tempAlignObj.GetPars(t,r);
 859 }
 860
 861 //-------------------------------------------------------------
 862 void AliITSAlignMille2Module::GetSensVolGlobalParams(UShort_t volid,const Double_t* loct, const Double_t* locr,Double_t *t, Double_t *r)
 863 {
 864   // return global parameters of the sensor volid modified by the localDelta params
 865   for (int i=3;i--;) t[i] = r[i] = 0.;
 866   if (SensVolMatrix(volid,fSensVolMatrix)) return;
 867   AliAlignObjParams tempAlignObj;
 868   tempAlignObj.SetTranslation(loct[0],loct[1],loct[2]);
 869   tempAlignObj.SetRotation(locr[0],locr[1],locr[2]);
 870   //
 871   tempAlignObj.GetMatrix(*fSensVolModifMatrix);      // obtain local delta
 872   fSensVolModifMatrix->MultiplyLeft( fSensVolMatrix ); // obtain global delta
 873   tempAlignObj.SetMatrix(*fSensVolModifMatrix);
 874   tempAlignObj.GetPars(t,r);                         // obtain global params
 875 }
 876
 877 //-------------------------------------------------------------
 878 void AliITSAlignMille2Module::GetSensVolLocalParams(UShort_t volid,const Double_t* loct,const Double_t* locr,Double_t *t, Double_t *r)
 879 {
 880   // return parameters of the sensor volid (modified by the localDelta params) in the current volume
 881   for (int i=3;i--;) t[i] = r[i] = 0.;
 882   if (SensVolMatrix(volid,fSensVolMatrix)) return;
 883   AliAlignObjParams tempAlignObj;
 884   tempAlignObj.SetTranslation(loct[0],loct[1],loct[2]);
 885   tempAlignObj.SetRotation(locr[0],locr[1],locr[2]);
 886   //
 887   tempAlignObj.GetMatrix(*fSensVolModifMatrix);      // obtain local delta
 888   fSensVolModifMatrix->MultiplyLeft( fSensVolMatrix ); // obtain global delta
 889   fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() ); // obtain delta in current volume
 890   tempAlignObj.SetMatrix(*fSensVolModifMatrix);
 891   tempAlignObj.GetPars(t,r);                         // obtain params
 892 }
 893
 894 //-------------------------------------------------------------
 895 void AliITSAlignMille2Module::SetParVals(Double_t *vl,Int_t npar)
 896 {
 897   // set parameters
 898   for (int i=TMath::Min(npar,(Int_t)fNParTot);i--;) fParVals[i] = vl[i];
 899 }
 900
 901 //-------------------------------------------------------------
 902 void AliITSAlignMille2Module::GetGeomParamsGlo(Double_t *pars)
 903 {
 904   // recompute parameters from local to global frame
 905   //
 906   // is there anything to do?
 907   if (GeomParamsGlobal()) {for (int i=kMaxParGeom;i--;) pars[i] = fParVals[i]; return;}
 908   //
 909   // IMPORTANT: It is assumed that the parents params are defined in a same way (local or global)
 910   // as for the current module. Since in the mp2 the modules are stored from parents to children,
 911   // it is safe to call this method in loop staring from the lowest level child, i.e. from the end
 912   // of the modules array.
 913   //
 914   // DeltaGlobal = (ModifParents)*DeltaLocal*(ModifParents)^-1
 915   //
 916   *fSensVolMatrix = *fMatrix;   // current global matrix
 917   AliAlignObjParams tempAlignObj;
 918   AliITSAlignMille2Module* parent = GetParent();
 919   while (parent) {
 920     if (parent->GeomParamsGlobal()) {
 921       AliError("Cannot convert params to Global when the parents are already Global\n");
 922       for (int i=kMaxParGeom;i--;) pars[i] = 0;
 923       return;
 924     }
 925     fSensVolMatrix->MultiplyLeft( &parent->GetMatrix()->Inverse() ); // Local Matrix
 926     Float_t *parpar = parent->GetParVals();
 927     tempAlignObj.SetTranslation(parpar[0],parpar[1],parpar[2]);
 928     tempAlignObj.SetRotation(parpar[3],parpar[4],parpar[5]);
 929     tempAlignObj.GetMatrix(*fSensVolModifMatrix);
 930     fSensVolMatrix->MultiplyLeft(fSensVolModifMatrix);
 931     fSensVolMatrix->MultiplyLeft(parent->GetMatrix());  // global matrix after parents modifications
 932     parent = parent->GetParent();
 933   }
 934   //
 935   tempAlignObj.SetTranslation(fParVals[0],fParVals[1],fParVals[2]);
 936   tempAlignObj.SetRotation(fParVals[3],fParVals[4],fParVals[5]);
 937   tempAlignObj.GetMatrix(*fSensVolModifMatrix);  // local delta matrix
 938   fSensVolModifMatrix->Multiply( &fSensVolMatrix->Inverse() );
 939   fSensVolModifMatrix->MultiplyLeft( fSensVolMatrix );
 940   tempAlignObj.SetMatrix( *fSensVolModifMatrix );  // global delta matrix
 941   tempAlignObj.GetPars(pars,pars+3);
 942   //
 943 }
 944
 945 //-------------------------------------------------------------
 946 void AliITSAlignMille2Module::GetGeomParamsLoc(Double_t *pars)
 947 {
 948   // recompute parameters from global to local frame
 949   //
 950   // is there anything to do?
 951   if (!GeomParamsGlobal()) {for (int i=kMaxParGeom;i--;) pars[i] = fParVals[i]; return;}
 952   //
 953   // IMPORTANT: It is assumed that the parents params are defined in a same way (local or global)
 954   // as for the current module. Since in the mp2 the modules are stored from parents to children,
 955   // it is safe to call this method in loop staring from the lowest level child, i.e. from the end
 956   // of the modules array.
 957   //
 958   //  DeltaLocal = (DeltaParents*GlobalMat)^-1*DeltaGlobal*(DeltaParents*GlobalMat)
 959   //
 960   AliITSAlignMille2Module* parent = GetParent();
 961   AliAlignObjParams tempAlignObj;
 962   tempAlignObj.SetTranslation(0.,0.,0.);
 963   tempAlignObj.SetRotation(0.,0.,0.);
 964   tempAlignObj.GetMatrix(*fSensVolMatrix); // get no-shift matrix
 965   //
 966   while (parent) { // accumulate the product of parents global modifications
 967     if (!parent->GeomParamsGlobal()) {
 968       AliError("Cannot convert params to Local when the parents are already Local\n");
 969       for (int i=kMaxParGeom;i--;) pars[i] = 0;
 970       return;
 971     }
 972     Float_t *parpar = parent->GetParVals();
 973     tempAlignObj.SetTranslation(parpar[0],parpar[1],parpar[2]);
 974     tempAlignObj.SetRotation(parpar[3],parpar[4],parpar[5]);
 975     tempAlignObj.GetMatrix(*fSensVolModifMatrix);
 976     fSensVolMatrix->Multiply(fSensVolModifMatrix);
 977     parent = parent->GetParent();
 978   }
 979   // global matrix after parents modifications
 980   fSensVolMatrix->Multiply(fMatrix);
 981   //
 982   tempAlignObj.SetTranslation(fParVals[0],fParVals[1],fParVals[2]);
 983   tempAlignObj.SetRotation(fParVals[3],fParVals[4],fParVals[5]);
 984   tempAlignObj.GetMatrix(*fSensVolModifMatrix);  // global delta matrix
 985   fSensVolModifMatrix->MultiplyLeft( &fSensVolMatrix->Inverse() );
 986   fSensVolModifMatrix->Multiply( fSensVolMatrix );
 987   tempAlignObj.SetMatrix( *fSensVolModifMatrix );  // local delta matrix
 988   tempAlignObj.GetPars(pars,pars+3);
 989   //
 990 }
 991
 992
 993 //-------------------------------------------------------------
 994 void AliITSAlignMille2Module::CalcDerivDPosDPar(Int_t sensVol,const Double_t* pl, Double_t *deriv)
 995 {
 996   // calculate jacobian of the global position vs Parameters (dPos/dParam)
 997   // for the point in the sensor sensVol
 998   const double kDel = 0.01;
 999   double pos0[3],pos1[3],pos2[3],pos3[3];
1000   double delta[kMaxParGeom];
1001   //
1002   for (int ip=kMaxParGeom;ip--;) delta[ip] = 0;
1003   //
1004   for (int ip=kMaxParGeom;ip--;) {
1005     //
1006     delta[ip] -= kDel;
1007     GetSensitiveVolumeModifiedMatrix(sensVol,delta,!GeomParamsGlobal())->LocalToMaster(pl,pos0);
1008     delta[ip] += kDel/2;
1009     GetSensitiveVolumeModifiedMatrix(sensVol,delta,!GeomParamsGlobal())->LocalToMaster(pl,pos1);
1010     delta[ip] += kDel;
1011     GetSensitiveVolumeModifiedMatrix(sensVol,delta,!GeomParamsGlobal())->LocalToMaster(pl,pos2);
1012     delta[ip] += kDel/2;
1013     GetSensitiveVolumeModifiedMatrix(sensVol,delta,!GeomParamsGlobal())->LocalToMaster(pl,pos3);
1014     //
1015     delta[ip] = 0;
1016     double *curd = deriv + ip*3;
1017     for (int i=3;i--;) curd[i] = (8.*(pos2[i]-pos1[i]) - (pos3[i]-pos0[i]))/6./kDel;
1018   }
1019   //
1020 }
1021
1022 //-------------------------------------------------------------
1023 void AliITSAlignMille2Module::CalcDerivGloLoc(Int_t idx,Double_t *deriv)
1024 {
1025   // calculate derivative of global params vs local param idx:  deriv[j] = dParGlo[j]/dParLoc[idx]
1026   Double_t lpar[kMaxParGeom];
1027   for (int i=kMaxParGeom;i--;) lpar[i] = 0.;
1028   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1029   Double_t par1[kMaxParGeom]; // f(x-h)
1030   Double_t par2[kMaxParGeom]; // f(x-h/2)
1031   Double_t par3[kMaxParGeom]; // f(x+h/2)
1032   Double_t par4[kMaxParGeom]; // f(x+h)
1033   //
1034   const Double_t dpar = 1e-3;
1035   //
1036   // first values
1037   lpar[idx] -= dpar;
1038   GetGlobalParams(lpar,lpar+3, par1,par1+3);
1039   //
1040   // second values
1041   lpar[idx] += dpar/2;
1042   GetGlobalParams(lpar,lpar+3, par2,par2+3);
1043   //
1044   // third values
1045   lpar[idx] += dpar;
1046   GetGlobalParams(lpar,lpar+3, par3,par3+3);
1047   //
1048   // fourth values
1049   lpar[idx] += dpar/2;
1050   GetGlobalParams(lpar,lpar+3, par4,par4+3);
1051   //
1052   Double_t h2 = 1./(2.*dpar);
1053   for (int i=kMaxParGeom;i--;) {
1054     Double_t d0 = par4[i]-par1[i];
1055     Double_t d2 = 2.*(par3[i]-par2[i]);
1056     deriv[i] = h2*(4*d2 - d0)/3.;
1057     if (TMath::Abs(deriv[i]) < 1.0e-9) deriv[i] = 0.0;
1058   }
1059   //
1060 }
1061
1062 //-------------------------------------------------------------
1063 void AliITSAlignMille2Module::CalcDerivLocGlo(Double_t *deriv)
1064 {
1065   // calculate derivative of local params vs global params:  deriv[i][j] = dParLoc[i]/dParGlo[j]
1066   Double_t gpar[kMaxParGeom];
1067   for (int i=kMaxParGeom;i--;) gpar[i] = 0.;
1068   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1069   Double_t par1[kMaxParGeom]; // f(x-h)
1070   Double_t par2[kMaxParGeom]; // f(x-h/2)
1071   Double_t par3[kMaxParGeom]; // f(x+h/2)
1072   Double_t par4[kMaxParGeom]; // f(x+h)
1073   //
1074   const Double_t dpar = 1e-3;
1075   //
1076   for (int ig=kMaxParGeom;ig--;) {
1077     // first values
1078     gpar[ig] -= dpar;
1079     GetLocalParams(gpar,gpar+3, par1,par1+3);
1080     //
1081     // second values
1082     gpar[ig] += dpar/2;
1083     GetLocalParams(gpar,gpar+3, par2,par2+3);
1084     //
1085     // third values
1086     gpar[ig] += dpar;
1087     GetLocalParams(gpar,gpar+3, par3,par3+3);
1088     //
1089     // fourth values
1090     gpar[ig] += dpar/2;
1091     GetLocalParams(gpar,gpar+3, par4,par4+3);
1092     //
1093     Double_t h2 = 1./(2.*dpar);
1094     for (int i=kMaxParGeom;i--;) {
1095       Double_t d0 = par4[i]-par1[i];
1096       Double_t d2 = 2.*(par3[i]-par2[i]);
1097       int idig = i*kMaxParGeom + ig;
1098       deriv[idig] = h2*(4*d2 - d0)/3.;
1099       if (TMath::Abs(deriv[idig]) < 1.0e-9) deriv[idig] = 0.0;
1100     }
1101   }
1102   //
1103 }
1104
1105 //________________________________________________________________________________________________________
1106 void AliITSAlignMille2Module::CalcDerivGloLoc(Int_t sensVol,Int_t paridx,Double_t* derivative)
1107 {
1108   /// calculate numerically the derivatives of global params vs local param paridx for sensor sensVol: dPglob/dPloc_paridx
1109   //
1110   Double_t lpar[kMaxParGeom];
1111   for (int i=kMaxParGeom;i--;) lpar[i] = 0.;
1112   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1113   Double_t par1[kMaxParGeom]; // f(x-h)
1114   Double_t par2[kMaxParGeom]; // f(x-h/2)
1115   Double_t par3[kMaxParGeom]; // f(x+h/2)
1116   Double_t par4[kMaxParGeom]; // f(x+h)
1117   //
1118   const Double_t dpar = 1e-3;
1119   //
1120   // first values
1121   lpar[paridx] -= dpar;
1122   GetSensVolGlobalParams(sensVol,lpar,lpar+3, par1,par1+3);
1123   //
1124   // second values
1125   lpar[paridx] += dpar/2;
1126   GetSensVolGlobalParams(sensVol,lpar,lpar+3, par2,par2+3);
1127   //
1128   // third values
1129   lpar[paridx] += dpar;
1130   GetSensVolGlobalParams(sensVol,lpar,lpar+3, par3,par3+3);
1131   //
1132   // fourth values
1133   lpar[paridx] += dpar/2;
1134   GetSensVolGlobalParams(sensVol,lpar,lpar+3, par4,par4+3);
1135   //
1136   Double_t h2 = 1./(2.*dpar);
1137   for (int i=kMaxParGeom;i--;) {
1138     Double_t d0 = par4[i]-par1[i];
1139     Double_t d2 = 2.*(par3[i]-par2[i]);
1140     derivative[i] = h2*(4*d2 - d0)/3.;
1141     if (TMath::Abs(derivative[i]) < 1.0e-9) derivative[i] = 0.0;
1142   }
1143   //
1144 }
1145
1146 //________________________________________________________________________________________________________
1147 void AliITSAlignMille2Module::CalcDerivCurLoc(Int_t sensVol,Int_t paridx,Double_t* derivative)
1148 {
1149   /// calculate numerically the derivatives of sensor params in the current volume vs sensor local param paridx
1150   //
1151   Double_t lpar[kMaxParGeom];
1152   for (int i=kMaxParGeom;i--;) lpar[i] = 0.;
1153   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1154   Double_t par1[kMaxParGeom]; // f(x-h)
1155   Double_t par2[kMaxParGeom]; // f(x-h/2)
1156   Double_t par3[kMaxParGeom]; // f(x+h/2)
1157   Double_t par4[kMaxParGeom]; // f(x+h)
1158   //
1159   const Double_t dpar = 1e-3;
1160   //
1161   // first values
1162   lpar[paridx] -= dpar;
1163   GetSensVolLocalParams(sensVol,lpar,lpar+3, par1,par1+3);
1164   //
1165   // second values
1166   lpar[paridx] += dpar/2;
1167   GetSensVolLocalParams(sensVol,lpar,lpar+3, par2,par2+3);
1168   //
1169   // third values
1170   lpar[paridx] += dpar;
1171   GetSensVolLocalParams(sensVol,lpar,lpar+3, par3,par3+3);
1172   //
1173   // fourth values
1174   lpar[paridx] += dpar/2;
1175   GetSensVolLocalParams(sensVol,lpar,lpar+3, par4,par4+3);
1176   //
1177   Double_t h2 = 1./(2.*dpar);
1178   for (int i=kMaxParGeom;i--;) {
1179     Double_t d0 = par4[i]-par1[i];
1180     Double_t d2 = 2.*(par3[i]-par2[i]);
1181     derivative[i] = h2*(4*d2 - d0)/3.;
1182     if (TMath::Abs(derivative[i]) < 1.0e-9) derivative[i] = 0.0;
1183   }
1184   //
1185 }
1186
1187
1188 //-------------------------------------------------------------
1189 void AliITSAlignMille2Module::GetGlobalParams(Double_t *t, Double_t *r) const
1190 {
1191   // global parameters of the module
1192   AliAlignObjParams tempAlignObj;
1193   tempAlignObj.SetMatrix( *fMatrix );
1194   tempAlignObj.GetPars(t,r);
1195 }
1196
1197 //-------------------------------------------------------------
1198 void AliITSAlignMille2Module::GetGlobalParams(const Double_t* loct, const Double_t* locr, Double_t *t, Double_t *r)
1199 {
1200   // global parameters of the module after the modification by local loct,locr
1201   AliAlignObjParams tempAlignObj;
1202   tempAlignObj.SetTranslation(loct[0],loct[1],loct[2]);
1203   tempAlignObj.SetRotation(locr[0],locr[1],locr[2]);
1204   tempAlignObj.GetMatrix(*fSensVolModifMatrix);
1205   *fSensVolMatrix = *fMatrix;
1206   fSensVolMatrix->Multiply(fSensVolModifMatrix);
1207   tempAlignObj.SetMatrix(*fSensVolMatrix);
1208   tempAlignObj.GetPars(t,r);
1209 }
1210
1211 //-------------------------------------------------------------
1212 void AliITSAlignMille2Module::GetLocalParams(const Double_t* glot, const Double_t* glor, Double_t *t, Double_t *r)
1213 {
1214   // obtain local delta parameters from global delta params
1215   AliAlignObjParams tempAlignObj;
1216   tempAlignObj.SetTranslation(glot[0],glot[1],glot[2]);
1217   tempAlignObj.SetRotation(glor[0],glor[1],glor[2]);
1218   tempAlignObj.GetMatrix(*fSensVolMatrix);
1219   fSensVolMatrix->Multiply( fMatrix );
1220   fSensVolMatrix->MultiplyLeft( &fMatrix->Inverse() );
1221   tempAlignObj.SetMatrix(*fSensVolMatrix);
1222   tempAlignObj.GetPars(t,r);
1223 }