end-of-line normalization
[u/mrichter/AliRoot.git] / 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 }