Adding misalignment with gaussian distribution
[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 // Class AliMUONGeometryMisAligner 
19 // ----------------------------
20 // Class for misalignment of geometry transformations
21 //
22 // Authors: Bruce Becker, Javier Castillo
23
24 //__________________________________________________________________
25 //
26 /////////////////////////////////////////////////////////////////////
27 //This performs the misalignment on an existing muon arm geometry
28 //  based on the standard definition of the detector elements in 
29 //  $ALICE_ROOT/MUON/data
30 //
31 //  --> User has to specify the magnitude of the alignments, in the Cartesian 
32 //  co-ordiantes (which are used to apply translation misalignments) and in the
33 //  spherical co-ordinates (which are used to apply angular displacements)
34 //  --> If the constructor is used with no arguments, user has to set 
35 //  misalignment ranges by hand using the methods : 
36 //  SetApplyMisAlig, SetMaxCartMisAlig, SetMaxAngMisAlig, SetXYAngMisAligFactor
37 //  (last method takes account of the fact that the misalingment is greatest in 
38 //  the XY plane, since the detection elements are fixed to a support structure
39 //  in this plane. Misalignments in the XZ and YZ plane will be very small 
40 //  compared to those in the XY plane, which are small already - of the order 
41 //  of microns)
42
43 //  Note : If the detection elements are allowed to be misaligned in all
44 //  directions, this has consequences for the alignment algorithm
45 //  (AliMUONAlignment), which needs to know the number of free parameters. 
46 //  Eric only allowed 3 :  x,y,theta_xy, but in principle z and the other 
47 //  two angles are alignable as well.
48
49 #include "AliMUONGeometryMisAligner.h"
50 #include "AliMUONGeometryTransformer.h"
51 #include "AliMUONGeometryModuleTransformer.h"
52 #include "AliMUONGeometryDetElement.h"
53 #include "AliMUONGeometryStore.h"
54 #include "AliMUONGeometryBuilder.h"
55
56 #include "AliLog.h"
57
58 #include <TGeoMatrix.h>
59 #include <TMath.h>
60 #include <TRandom.h>
61
62 ClassImp(AliMUONGeometryMisAligner)
63 //______________________________________________________________________________
64 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartXMisAligM, Double_t cartXMisAligW, Double_t cartYMisAligM, Double_t cartYMisAligW, Double_t angMisAligM, Double_t angMisAligW)
65 :TObject(), fDisplacementGenerator(0)
66 {
67   /// Standard constructor
68   fCartXMisAligM = cartXMisAligM;
69   fCartXMisAligW = cartXMisAligW;       // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
70   fCartYMisAligM = cartYMisAligM;
71   fCartYMisAligW = cartYMisAligW;       // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
72   fAngMisAligM = angMisAligM;
73   fAngMisAligW = angMisAligW;
74   fXYAngMisAligFactor = 0.0;
75   fZCartMisAligFactor = 0.0;
76   fDisplacementGenerator = new TRandom(0);
77   fUseUni = kFALSE;
78   fUseGaus = kTRUE;
79 }
80
81 //______________________________________________________________________________
82 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartMisAligM, Double_t cartMisAligW, Double_t angMisAligM, Double_t angMisAligW)
83 :TObject(), fDisplacementGenerator(0)
84 {
85   /// Standard constructor
86   fCartXMisAligM = cartMisAligM;
87   fCartXMisAligW = cartMisAligW;        // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
88   fCartYMisAligM = cartMisAligM;
89   fCartYMisAligW = cartMisAligW;        // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
90   fAngMisAligM = angMisAligM;
91   fAngMisAligW = angMisAligW;
92   fXYAngMisAligFactor = 0.0;
93   fZCartMisAligFactor = 0.0;
94   fDisplacementGenerator = new TRandom(0);
95   fUseUni = kFALSE;
96   fUseGaus = kTRUE;
97 }
98
99 //______________________________________________________________________________
100 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner(Double_t cartMisAlig, Double_t angMisAlig)
101 :TObject(), fDisplacementGenerator(0)
102 {
103   /// Standard constructor
104   fCartXMisAligM = 0.;
105   fCartXMisAligW = cartMisAlig; // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
106   fCartYMisAligM = 0.;
107   fCartYMisAligW = cartMisAlig; // 0.5 mm. Perhaps this should go into AliMUONConstants.h ? 
108   fAngMisAligM = 0.;
109   fAngMisAligW = angMisAlig;
110   fXYAngMisAligFactor = 0.0;
111   fZCartMisAligFactor = 0.0;
112   fDisplacementGenerator = new TRandom(0);
113   fUseUni = kTRUE;
114   fUseGaus = kFALSE;
115 }
116
117 //_____________________________________________________________________________
118 AliMUONGeometryMisAligner::AliMUONGeometryMisAligner()
119 :TObject(), fDisplacementGenerator(0)
120 {
121 /// Default constructor
122 }
123
124 //______________________________________________________________________________
125 AliMUONGeometryMisAligner::
126 AliMUONGeometryMisAligner(const AliMUONGeometryMisAligner & right):
127 TObject(right)
128 {
129   /// Copy constructor (not implemented)
130
131   AliFatal("Copy constructor not provided.");
132 }
133
134 //______________________________________________________________________________
135 AliMUONGeometryMisAligner::~AliMUONGeometryMisAligner()
136 {
137 /// Destructor
138
139   delete fDisplacementGenerator;
140 }
141
142 //______________________________________________________________________________
143 AliMUONGeometryMisAligner & AliMUONGeometryMisAligner::
144 operator=(const AliMUONGeometryMisAligner & right)
145 {
146   /// Assignement operator (not implemented)
147
148   // check assignement to self
149   if (this == &right)
150     return *this;
151
152   AliFatal("Assignement operator not provided.");
153
154   return *this;
155 }
156
157 //_________________________________________________________________________
158 void
159 AliMUONGeometryMisAligner::SetXYAngMisAligFactor(Double_t factor)
160 {
161   /// Set XY angular misalign factor 
162
163   if (TMath::Abs(factor) > 1.0 && factor > 0.)
164     fXYAngMisAligFactor = factor;
165   else
166     AliError(Form("Invalid XY angular misalign factor, %d", factor));
167 }
168
169 //_________________________________________________________________________
170 void AliMUONGeometryMisAligner::SetZCartMisAligFactor(Double_t factor) 
171 {
172   /// Set XY angular misalign factor 
173   if (TMath::Abs(factor)<1.0 && factor>0.)
174     fZCartMisAligFactor = factor;
175   else
176     AliError(Form("Invalid Z cartesian misalign factor, %d", factor));
177 }
178
179 //_________________________________________________________________________
180 void AliMUONGeometryMisAligner::GetUniMisAlign(Double_t *cartMisAlig, Double_t *angMisAlig) const
181 {
182   /// Misalign using uniform distribution
183   /*
184     misalign the centre of the local transformation
185     rotation axes : 
186     fAngMisAlig[1,2,3] = [x,y,z]
187     Assume that misalignment about the x and y axes (misalignment of z plane)
188     is much smaller, since the entire detection plane has to be moved (the 
189     detection elements are on a support structure), while rotation of the x-y
190     plane is more free.
191   */
192   cartMisAlig[0] = fDisplacementGenerator->Uniform(-fCartXMisAligW+fCartXMisAligM, fCartXMisAligM+fCartXMisAligW);
193   cartMisAlig[1] = fDisplacementGenerator->Uniform(-fCartYMisAligW+fCartYMisAligM, fCartYMisAligM+fCartYMisAligW);
194   cartMisAlig[2] = fDisplacementGenerator->Uniform(-fZCartMisAligFactor*(fCartXMisAligW+fCartXMisAligM), fZCartMisAligFactor*(fCartXMisAligM+fCartXMisAligW));  
195  
196   angMisAlig[0] = fDisplacementGenerator->Uniform(-fXYAngMisAligFactor*(fAngMisAligW+fAngMisAligM), fXYAngMisAligFactor*(fAngMisAligM+fAngMisAligW));
197   angMisAlig[1] = fDisplacementGenerator->Uniform(-fXYAngMisAligFactor*(fAngMisAligW+fAngMisAligM), fXYAngMisAligFactor*(fAngMisAligM+fAngMisAligW));
198   angMisAlig[2] = fDisplacementGenerator->Uniform(-fAngMisAligW+fAngMisAligM, fAngMisAligM+fAngMisAligW);       // degrees
199 }
200
201 //_________________________________________________________________________
202 void AliMUONGeometryMisAligner::GetGausMisAlign(Double_t cartMisAlig[3], Double_t angMisAlig[3]) const
203 {
204   /// Misalign using gaussian distribution
205   /*
206     misalign the centre of the local transformation
207     rotation axes : 
208     fAngMisAlig[1,2,3] = [x,y,z]
209     Assume that misalignment about the x and y axes (misalignment of z plane)
210     is much smaller, since the entire detection plane has to be moved (the 
211     detection elements are on a support structure), while rotation of the x-y
212     plane is more free.
213   */
214   cartMisAlig[0] = fDisplacementGenerator->Gaus(fCartXMisAligM, fCartXMisAligW);
215   cartMisAlig[1] = fDisplacementGenerator->Gaus(fCartYMisAligM, fCartYMisAligW);
216   cartMisAlig[2] = fDisplacementGenerator->Gaus(fCartXMisAligM, fZCartMisAligFactor*fCartXMisAligW);
217  
218   angMisAlig[0] = fDisplacementGenerator->Gaus(fAngMisAligM, fXYAngMisAligFactor*fAngMisAligW);
219   angMisAlig[1] = fDisplacementGenerator->Gaus(fAngMisAligM, fXYAngMisAligFactor*fAngMisAligW);
220   angMisAlig[2] = fDisplacementGenerator->Gaus(fAngMisAligM, fAngMisAligW);     // degrees
221 }
222
223 //_________________________________________________________________________
224 TGeoCombiTrans AliMUONGeometryMisAligner::MisAlign(const TGeoCombiTrans & transform) const
225 {
226   /// Misalign given transformation and return the misaligned transformation
227   
228   Double_t cartMisAlig[3] = {0,0,0};
229   Double_t angMisAlig[3] = {0,0,0};
230   const Double_t *trans = transform.GetTranslation();
231   TGeoRotation *rot;
232   // check if the rotation we obtain is not NULL
233   if (transform.GetRotation())
234     {
235       rot = transform.GetRotation();
236     }
237   else
238     {
239       rot = new TGeoRotation("rot");
240     }                   // default constructor.
241
242   if (fUseUni) { 
243     GetUniMisAlign(cartMisAlig,angMisAlig);
244   }
245   else { 
246     if (!fUseGaus) {
247       AliWarning("Neither uniform nor gausian distribution is set! Will use gausian...");
248     } 
249     GetGausMisAlign(cartMisAlig,angMisAlig);
250   }
251
252   TGeoTranslation newTrans(cartMisAlig[0] + trans[0], cartMisAlig[1] + trans[1], cartMisAlig[2] + trans[2]);
253   
254   AliInfo(Form("Rotated by %f about Z axis.", angMisAlig[2]));
255   rot->RotateX(angMisAlig[0]);
256   rot->RotateY(angMisAlig[1]);
257   rot->RotateZ(angMisAlig[2]);
258
259   return TGeoCombiTrans(newTrans, *rot);
260 }
261
262 //______________________________________________________________________
263 AliMUONGeometryTransformer *
264 AliMUONGeometryMisAligner::MisAlign(const AliMUONGeometryTransformer *
265                                  transformer, Bool_t verbose)
266 {
267   /////////////////////////////////////////////////////////////////////
268   //   Takes the internal geometry module transformers, copies them
269   // and gets the Detection Elements from them.
270   // Calculates misalignment parameters and applies these
271   // to the local transform of the Detection Element
272   // Obtains the global transform by multiplying the module transformer
273   // transformation with the local transformation 
274   // Applies the global transform to a new detection element
275   // Adds the new detection element to a new module transformer
276   // Adds the new module transformer to a new geometry transformer
277   // Returns the new geometry transformer
278
279
280   AliMUONGeometryTransformer *newGeometryTransformer =
281     new AliMUONGeometryTransformer(kTRUE);
282   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++)
283     {                           // module transformers
284       
285       const AliMUONGeometryModuleTransformer *kModuleTransformer =
286         transformer->GetModuleTransformer(iMt, true);
287       
288       AliMUONGeometryModuleTransformer *newModuleTransformer =
289         new AliMUONGeometryModuleTransformer(iMt);
290       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
291
292       TGeoCombiTrans moduleTransform =
293         TGeoCombiTrans(*kModuleTransformer->GetTransformation());
294       TGeoCombiTrans *newModuleTransform = new TGeoCombiTrans(moduleTransform); 
295               // same module transform as the previous one 
296               // no mis align object created
297       newModuleTransformer->SetTransformation(moduleTransform);
298
299       AliMUONGeometryStore *detElements =
300         kModuleTransformer->GetDetElementStore();
301
302       if (verbose)
303         AliInfo(Form
304                 ("%i DEs in old GeometryStore  %i",
305                  detElements->GetNofEntries(), iMt));
306
307       for (Int_t iDe = 0; iDe < detElements->GetNofEntries(); iDe++)
308         {                       // detection elements.
309           AliMUONGeometryDetElement *detElement =
310             (AliMUONGeometryDetElement *) detElements->GetEntry(iDe);
311           if (!detElement)
312             AliFatal("Detection element not found.");
313
314           /// make a new detection element
315           AliMUONGeometryDetElement *newDetElement =
316             new AliMUONGeometryDetElement(detElement->GetId(),
317                                           detElement->GetVolumePath());
318
319           // local transformation of this detection element.
320           TGeoCombiTrans localTransform
321             = TGeoCombiTrans(*detElement->GetLocalTransformation());
322           TGeoCombiTrans newLocalTransform = MisAlign(localTransform);
323           newDetElement->SetLocalTransformation(newLocalTransform);                                       
324
325           // global transformation
326           TGeoHMatrix newGlobalTransform =
327             AliMUONGeometryBuilder::Multiply(*newModuleTransform,
328                                               newLocalTransform);
329           newDetElement->SetGlobalTransformation(newGlobalTransform);
330           
331           // add this det element to module
332           newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
333                                                           newDetElement);
334           // Get delta transformation: 
335           // Tdelta = Tnew * Told.inverse
336           TGeoHMatrix  deltaTransform
337             = AliMUONGeometryBuilder::Multiply(
338                 newGlobalTransform, 
339                 detElement->GetGlobalTransformation()->Inverse());
340           
341           // Create mis alignment matrix
342           newGeometryTransformer
343             ->AddMisAlignDetElement(detElement->GetId(), deltaTransform);
344         }
345       if (verbose)
346         AliInfo(Form("Added module transformer %i to the transformer", iMt));
347       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
348     }
349   return newGeometryTransformer;
350 }
351
352
353
354