Fix coverity defect
[u/mrichter/AliRoot.git] / MUON / AliMUONGeometryMisAligner.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *      SigmaEffect_thetadegrees                                          * 
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 purpeateose. It is      *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // $Id$
17 //
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONGeometryMisAligner
20 ///
21 /// This performs the misalignment on an existing muon arm geometry
22 /// based on the standard definition of the detector elements in 
23 /// $ALICE_ROOT/MUON/data
24 ///
25 /// --> User has to specify the magnitude of the alignments, in the Cartesian 
26 /// co-ordiantes (which are used to apply translation misalignments) and in the
27 /// spherical co-ordinates (which are used to apply angular displacements)
28 ///
29 /// --> If the constructor is used with no arguments, user has to set 
30 /// misalignment ranges by hand using the methods : 
31 /// SetApplyMisAlig, SetMaxCartMisAlig, SetMaxAngMisAlig, SetXYAngMisAligFactor
32 /// (last method takes account of the fact that the misalingment is greatest in 
33 /// the XY plane, since the detection elements are fixed to a support structure
34 /// in this plane. Misalignments in the XZ and YZ plane will be very small 
35 /// compared to those in the XY plane, which are small already - of the order 
36 /// of microns)
37 ///
38 /// Note : If the detection elements are allowed to be misaligned in all
39 /// directions, this has consequences for the alignment algorithm
40 /// (AliMUONAlignment), which needs to know the number of free parameters. 
41 /// Eric only allowed 3 :  x,y,theta_xy, but in principle z and the other 
42 /// two angles are alignable as well.
43 ///
44 /// \author Bruce Becker, Javier Castillo
45 //-----------------------------------------------------------------------------
46
47 #include "AliMUONGeometryMisAligner.h"
48 #include "AliMUONGeometryTransformer.h"
49 #include "AliMUONGeometryModuleTransformer.h"
50 #include "AliMUONGeometryDetElement.h"
51 #include "AliMUONGeometryBuilder.h"
52 #include "AliMpExMap.h"
53 #include "AliMpExMapIterator.h"
54
55 #include "AliAlignObjMatrix.h"
56 #include "AliMathBase.h"
57 #include "AliLog.h"
58
59 #include <TClonesArray.h>
60 #include <TGeoMatrix.h>
61 #include <TMatrixDSym.h>
62 #include <TMath.h>
63 #include <TRandom.h>
64 #include <Riostream.h>
65
66 /// \cond CLASSIMP
67 ClassImp(AliMUONGeometryMisAligner)
68 /// \endcond
69
70 //______________________________________________________________________________
71 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartXMisAligM, Double_t cartXMisAligW, Double_t cartYMisAligM, Double_t cartYMisAligW, Double_t angMisAligM, Double_t angMisAligW)
72   : TObject(),
73     fUseUni(kFALSE),
74     fUseGaus(kTRUE),
75     fXYAngMisAligFactor(0.0),
76     fZCartMisAligFactor(0.0)
77 {
78   /// Standard constructor
79   for (Int_t i=0; i<6; i++){
80     for (Int_t j=0; j<2; j++){
81       fDetElemMisAlig[i][j] = 0.0;
82       fModuleMisAlig[i][j] = 0.0;
83     }
84   }
85   fDetElemMisAlig[0][0] = cartXMisAligM; 
86   fDetElemMisAlig[0][1] = cartXMisAligW; 
87   fDetElemMisAlig[1][0] = cartYMisAligM; 
88   fDetElemMisAlig[1][1] = cartYMisAligW; 
89   fDetElemMisAlig[5][0] = angMisAligM; 
90   fDetElemMisAlig[5][1] = angMisAligW;
91
92 }
93
94 //______________________________________________________________________________
95 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartMisAligM, Double_t cartMisAligW, Double_t angMisAligM, Double_t angMisAligW)
96   : TObject(), 
97     fUseUni(kFALSE),
98     fUseGaus(kTRUE),
99     fXYAngMisAligFactor(0.0),
100     fZCartMisAligFactor(0.0)
101 {
102   /// Standard constructor
103   for (Int_t i=0; i<6; i++){
104     for (Int_t j=0; j<2; j++){
105       fDetElemMisAlig[i][j] = 0.0;
106       fModuleMisAlig[i][j] = 0.0;
107     }
108   }
109   fDetElemMisAlig[0][0] = cartMisAligM; 
110   fDetElemMisAlig[0][1] = cartMisAligW; 
111   fDetElemMisAlig[1][0] = cartMisAligM; 
112   fDetElemMisAlig[1][1] = cartMisAligW; 
113   fDetElemMisAlig[5][0] = angMisAligM; 
114   fDetElemMisAlig[5][1] = angMisAligW;
115
116 }
117
118 //______________________________________________________________________________
119 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartMisAlig, Double_t angMisAlig)
120   : TObject(), 
121     fUseUni(kTRUE),
122     fUseGaus(kFALSE),
123     fXYAngMisAligFactor(0.0),
124     fZCartMisAligFactor(0.0)
125 {
126   /// Standard constructor
127   for (Int_t i=0; i<6; i++){
128     for (Int_t j=0; j<2; j++){
129       fDetElemMisAlig[i][j] = 0.0;
130       fModuleMisAlig[i][j] = 0.0;
131     }
132   }
133   fDetElemMisAlig[0][1] = cartMisAlig; 
134   fDetElemMisAlig[1][1] = cartMisAlig; 
135   fDetElemMisAlig[5][1] = angMisAlig;
136
137 }
138
139 //_____________________________________________________________________________
140 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner()
141   : TObject(), 
142     fUseUni(kTRUE),
143     fUseGaus(kFALSE),
144     fXYAngMisAligFactor(0.0),
145     fZCartMisAligFactor(0.0)
146 {
147   /// Default constructor
148   for (Int_t i=0; i<6; i++){
149     for (Int_t j=0; j<2; j++){
150       fDetElemMisAlig[i][j] = 0.0;
151       fModuleMisAlig[i][j] = 0.0;
152     }
153   }
154 }
155
156 //______________________________________________________________________________
157 AliMUONGeometryMisAligner::~AliMUONGeometryMisAligner()
158 {
159 /// Destructor
160
161 }
162
163 //_________________________________________________________________________
164 void
165 AliMUONGeometryMisAligner::SetXYAngMisAligFactor(Double_t factor)
166 {
167   /// Set XY angular misalign factor 
168
169   if (TMath::Abs(factor) > 1.0 && factor > 0.){
170     fXYAngMisAligFactor = factor;
171     fDetElemMisAlig[3][0] = fDetElemMisAlig[5][0]*factor; // These lines were 
172     fDetElemMisAlig[3][1] = fDetElemMisAlig[5][1]*factor; // added to keep
173     fDetElemMisAlig[4][0] = fDetElemMisAlig[5][0]*factor; // backward 
174     fDetElemMisAlig[4][1] = fDetElemMisAlig[5][1]*factor; // compatibility 
175   }
176   else
177     AliError(Form("Invalid XY angular misalign factor, %f", factor));
178 }
179
180 //_________________________________________________________________________
181 void AliMUONGeometryMisAligner::SetZCartMisAligFactor(Double_t factor) 
182 {
183   /// Set XY angular misalign factor 
184   if (TMath::Abs(factor)<1.0 && factor>0.) {
185     fZCartMisAligFactor = factor;
186     fDetElemMisAlig[2][0] = fDetElemMisAlig[0][0];        // These lines were added to 
187     fDetElemMisAlig[2][1] = fDetElemMisAlig[0][1]*factor; // keep backward compatibility
188   }
189   else
190     AliError(Form("Invalid Z cartesian misalign factor, %f", factor));
191 }
192
193 //_________________________________________________________________________
194 void AliMUONGeometryMisAligner::GetUniMisAlign(Double_t cartMisAlig[3], Double_t angMisAlig[3], const Double_t lParMisAlig[6][2]) const
195 {
196   /// Misalign using uniform distribution
197   /**
198     misalign the centre of the local transformation
199     rotation axes : 
200     fAngMisAlig[1,2,3] = [x,y,z]
201     Assume that misalignment about the x and y axes (misalignment of z plane)
202     is much smaller, since the entire detection plane has to be moved (the 
203     detection elements are on a support structure), while rotation of the x-y
204     plane is more free.
205   */
206   cartMisAlig[0] = gRandom->Uniform(-lParMisAlig[0][1]+lParMisAlig[0][0], lParMisAlig[0][0]+lParMisAlig[0][1]);
207   cartMisAlig[1] = gRandom->Uniform(-lParMisAlig[1][1]+lParMisAlig[1][0], lParMisAlig[1][0]+lParMisAlig[1][1]);
208   cartMisAlig[2] = gRandom->Uniform(-lParMisAlig[2][1]+lParMisAlig[2][0], lParMisAlig[2][0]+lParMisAlig[2][1]);  
209  
210   angMisAlig[0] = gRandom->Uniform(-lParMisAlig[3][1]+lParMisAlig[3][0], lParMisAlig[3][0]+lParMisAlig[3][1]);
211   angMisAlig[1] = gRandom->Uniform(-lParMisAlig[4][1]+lParMisAlig[4][0], lParMisAlig[4][0]+lParMisAlig[4][1]);
212   angMisAlig[2] = gRandom->Uniform(-lParMisAlig[5][1]+lParMisAlig[5][0], lParMisAlig[5][0]+lParMisAlig[5][1]);  // degrees
213 }
214
215 //_________________________________________________________________________
216 void AliMUONGeometryMisAligner::GetGausMisAlign(Double_t cartMisAlig[3], Double_t angMisAlig[3], const Double_t lParMisAlig[6][2]) const
217 {
218   /// Misalign using gaussian distribution
219   /**
220     misalign the centre of the local transformation
221     rotation axes : 
222     fAngMisAlig[1,2,3] = [x,y,z]
223     Assume that misalignment about the x and y axes (misalignment of z plane)
224     is much smaller, since the entire detection plane has to be moved (the 
225     detection elements are on a support structure), while rotation of the x-y
226     plane is more free.
227   */
228   cartMisAlig[0] = AliMathBase::TruncatedGaus(lParMisAlig[0][0], lParMisAlig[0][1], 3.*lParMisAlig[0][1]);
229   cartMisAlig[1] = AliMathBase::TruncatedGaus(lParMisAlig[1][0], lParMisAlig[1][1], 3.*lParMisAlig[1][1]);
230   cartMisAlig[2] = AliMathBase::TruncatedGaus(lParMisAlig[2][0], lParMisAlig[2][1], 3.*lParMisAlig[2][1]);
231  
232   angMisAlig[0] = AliMathBase::TruncatedGaus(lParMisAlig[3][0], lParMisAlig[3][1], 3.*lParMisAlig[3][1]);
233   angMisAlig[1] = AliMathBase::TruncatedGaus(lParMisAlig[4][0], lParMisAlig[4][1], 3.*lParMisAlig[4][1]);
234   angMisAlig[2] = AliMathBase::TruncatedGaus(lParMisAlig[5][0], lParMisAlig[5][1], 3.*lParMisAlig[5][1]);       // degrees
235 }
236
237 //_________________________________________________________________________
238 TGeoCombiTrans AliMUONGeometryMisAligner::MisAlignDetElem(const TGeoCombiTrans & transform) const
239 {
240   /// Misalign given transformation and return the misaligned transformation. 
241   /// Use misalignment parameters for detection elements.
242   /// Note that applied misalignments are small deltas with respect to the detection 
243   /// element own ideal local reference frame. Thus deltaTransf represents 
244   /// the transformation to go from the misaligned d.e. local coordinates to the 
245   /// ideal d.e. local coordinates. 
246   /// Also note that this -is not- what is in the ALICE alignment framework known 
247   /// as local nor global (see AliMUONGeometryMisAligner::MisAlign) 
248   
249   Double_t cartMisAlig[3] = {0,0,0};
250   Double_t angMisAlig[3] = {0,0,0};
251
252   if (fUseUni) { 
253     GetUniMisAlign(cartMisAlig,angMisAlig,fDetElemMisAlig);
254   }
255   else { 
256     if (!fUseGaus) {
257       AliWarning("Neither uniform nor gausian distribution is set! Will use gausian...");
258     } 
259     GetGausMisAlign(cartMisAlig,angMisAlig,fDetElemMisAlig);
260   }
261
262   TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
263   TGeoRotation deltaRot;
264   deltaRot.RotateX(angMisAlig[0]);
265   deltaRot.RotateY(angMisAlig[1]);
266   deltaRot.RotateZ(angMisAlig[2]);
267
268   TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
269   TGeoHMatrix newTransfMat = transform * deltaTransf;
270     
271   AliInfo(Form("Rotated DE by %f about Z axis.", angMisAlig[2]));
272
273   return TGeoCombiTrans(newTransfMat);
274 }
275
276 //_________________________________________________________________________
277 TGeoCombiTrans AliMUONGeometryMisAligner::MisAlignModule(const TGeoCombiTrans & transform) const
278 {
279   /// Misalign given transformation and return the misaligned transformation. 
280   /// Use misalignment parameters for modules.
281   /// Note that applied misalignments are small deltas with respect to the module 
282   /// own ideal local reference frame. Thus deltaTransf represents 
283   /// the transformation to go from the misaligned module local coordinates to the 
284   /// ideal module local coordinates. 
285   /// Also note that this -is not- what is in the ALICE alignment framework known 
286   /// as local nor global (see AliMUONGeometryMisAligner::MisAlign) 
287   
288   Double_t cartMisAlig[3] = {0,0,0};
289   Double_t angMisAlig[3] = {0,0,0};
290
291   if (fUseUni) { 
292     GetUniMisAlign(cartMisAlig,angMisAlig,fModuleMisAlig);
293   }
294   else { 
295     if (!fUseGaus) {
296       AliWarning("Neither uniform nor gausian distribution is set! Will use gausian...");
297     } 
298     GetGausMisAlign(cartMisAlig,angMisAlig,fModuleMisAlig);
299   }
300
301   TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
302   TGeoRotation deltaRot;
303   deltaRot.RotateX(angMisAlig[0]);
304   deltaRot.RotateY(angMisAlig[1]);
305   deltaRot.RotateZ(angMisAlig[2]);
306
307   TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
308   TGeoHMatrix newTransfMat = transform * deltaTransf;
309
310   AliInfo(Form("Rotated Module by %f about Z axis.", angMisAlig[2]));
311
312   return TGeoCombiTrans(newTransfMat);
313 }
314
315 //______________________________________________________________________
316 AliMUONGeometryTransformer *
317 AliMUONGeometryMisAligner::MisAlign(const AliMUONGeometryTransformer *
318                                     transformer, Bool_t verbose)
319 {
320   /// Takes the internal geometry module transformers, copies them to
321   /// new geometry module transformers. 
322   /// Calculates  module misalignment parameters and applies these
323   /// to the new module transformer.
324   /// Calculates the module misalignment delta transformation in the 
325   /// Alice Alignment Framework newTransf = delta * oldTransf.
326   /// Add a module misalignment to the new geometry transformer.
327   /// Gets the Detection Elements from the module transformer. 
328   /// Calculates misalignment parameters and applies these
329   /// to the local transformation of the Detection Element.
330   /// Obtains the new global transformation by multiplying the new 
331   /// module transformer transformation with the new local transformation. 
332   /// Applies the new global transform to a new detection element.
333   /// Adds the new detection element to a new module transformer.
334   /// Calculates the d.e. misalignment delta transformation in the 
335   /// Alice Alignment Framework (newGlobalTransf = delta * oldGlobalTransf).
336   /// Add a d.e. misalignment to the new geometry transformer.
337   /// Adds the new module transformer to a new geometry transformer.
338   /// Returns the new geometry transformer.
339
340
341   AliMUONGeometryTransformer *newGeometryTransformer =
342     new AliMUONGeometryTransformer();
343   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++)
344     {                           // module transformers
345       const AliMUONGeometryModuleTransformer *kModuleTransformer =
346         transformer->GetModuleTransformer(iMt, true);
347       
348       AliMUONGeometryModuleTransformer *newModuleTransformer =
349         new AliMUONGeometryModuleTransformer(iMt);
350       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
351
352       TGeoCombiTrans moduleTransform =
353         TGeoCombiTrans(*kModuleTransformer->GetTransformation());
354       // New module transformation
355       TGeoCombiTrans newModuleTransform = MisAlignModule(moduleTransform);
356       newModuleTransformer->SetTransformation(newModuleTransform);
357
358       // Get delta transformation: 
359       // Tdelta = Tnew * Told.inverse
360       TGeoHMatrix deltaModuleTransform = 
361         AliMUONGeometryBuilder::Multiply(
362           newModuleTransform, 
363           kModuleTransformer->GetTransformation()->Inverse());
364
365       // Create module mis alignment matrix
366       newGeometryTransformer
367         ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
368
369       AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
370
371       if (verbose)
372         AliInfo(Form("%i DEs in old GeometryStore  %i",detElements->GetSize(), iMt));
373
374       TIter next(detElements->CreateIterator());
375       AliMUONGeometryDetElement *detElement;
376       
377       while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
378       {
379           /// make a new detection element
380           AliMUONGeometryDetElement *newDetElement =
381             new AliMUONGeometryDetElement(detElement->GetId(),
382                                           detElement->GetVolumePath());
383
384           // local transformation of this detection element.
385           TGeoCombiTrans localTransform
386             = TGeoCombiTrans(*detElement->GetLocalTransformation());
387           TGeoCombiTrans newLocalTransform = MisAlignDetElem(localTransform);
388           newDetElement->SetLocalTransformation(newLocalTransform);
389
390
391           // global transformation
392           TGeoHMatrix newGlobalTransform =
393             AliMUONGeometryBuilder::Multiply(newModuleTransform,
394                                              newLocalTransform);
395           newDetElement->SetGlobalTransformation(newGlobalTransform);
396           
397           // add this det element to module
398           newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
399                                                           newDetElement);
400
401           // In the Alice Alignment Framework misalignment objects store
402           // global delta transformation
403           // Get detection "intermediate" global transformation
404           TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform;
405           // Get detection element global delta transformation: 
406           // Tdelta = Tnew * Told.inverse
407           TGeoHMatrix  deltaGlobalTransform
408             = AliMUONGeometryBuilder::Multiply(
409                 newGlobalTransform, 
410                 newOldGlobalTransform.Inverse());
411           
412           // Create mis alignment matrix
413           newGeometryTransformer
414             ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
415         }
416
417       
418       if (verbose)
419         AliInfo(Form("Added module transformer %i to the transformer", iMt));
420       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
421     }
422   return newGeometryTransformer;
423 }
424
425
426 void AliMUONGeometryMisAligner::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY){
427
428   Int_t chIdMin = (rChId<0)? 0 : rChId;
429   Int_t chIdMax = (rChId<0)? 9 : rChId;
430   Double_t chResX = (rChResX<0)? fModuleMisAlig[0][1] : rChResX;
431   Double_t chResY = (rChResY<0)? fModuleMisAlig[1][1] : rChResY;
432   Double_t deResX = (rDeResX<0)? fDetElemMisAlig[0][1] : rDeResX;
433   Double_t deResY = (rDeResY<0)? fDetElemMisAlig[1][1] : rDeResY;
434
435   TMatrixDSym mChCorrMatrix(6);
436   mChCorrMatrix[0][0]=chResX*chResX;
437   mChCorrMatrix[1][1]=chResY*chResY;
438   //  mChCorrMatrix.Print();
439
440   TMatrixDSym mDECorrMatrix(6);
441   mDECorrMatrix[0][0]=deResX*deResX;
442   mDECorrMatrix[1][1]=deResY*deResY;
443   //  mDECorrMatrix.Print();
444
445   AliAlignObjMatrix *alignMat = 0x0;
446
447   for(Int_t chId=chIdMin; chId<=chIdMax; chId++) {
448     TString chName1;
449     TString chName2;
450     if (chId<4){
451       chName1 = Form("GM%d",chId);
452       chName2 = Form("GM%d",chId);
453     } else {
454       chName1 = Form("GM%d",4+(chId-4)*2);
455       chName2 = Form("GM%d",4+(chId-4)*2+1);
456     }
457     
458     for (int i=0; i<misAlignArray->GetEntries(); i++) {
459       alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
460       TString volName(alignMat->GetSymName());
461       if((volName.Contains(chName1)&&
462           ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
463            (volName.Length()==volName.Index(chName1)+chName1.Length())))||
464          (volName.Contains(chName2)&&
465           ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
466            (volName.Length()==volName.Index(chName2)+chName2.Length())))){
467         volName.Remove(0,volName.Last('/')+1);
468         if (volName.Contains("GM")) {
469           //    alignMat->Print("NULL");
470           alignMat->SetCorrMatrix(mChCorrMatrix);
471         } else if (volName.Contains("DE")) {
472           //    alignMat->Print("NULL");
473           alignMat->SetCorrMatrix(mDECorrMatrix);
474         }
475       }
476     }
477   }
478 }
479
480
481